Penyepelasan kuantum Krylov berasaskan sampel bagi model kekisi fermionic

Anggaran penggunaan: Sembilan saat pada pemproses Heron r2 (NOTA: Ini adalah anggaran sahaja. Masa jalan anda mungkin berbeza.)

Hasil pembelajaran

Selepas melalui tutorial ini, pengguna seharusnya memahami:

• Cara menggunakan SQD Qiskit addon untuk menganggarkan tenaga keadaan dasar model kekisi menggunakan rentetan bit yang disampling daripada unit pemprosesan kuantum (QPU).
• Cara menggunakan ffsim untuk membina litar evolusi masa bagi simulasi fermionic.
• Cara menggabungkan sampel daripada pelbagai litar untuk pasca-pemprosesan dengan algoritma penyepelasan Krylov berasaskan sampel (SKQD).

Prasyarat

Kami mencadangkan pengguna biasa dengan topik-topik berikut sebelum melalui tutorial ini:

• Penyepelasan kuantum berasaskan sampel bagi Hamiltonian kimia
• Penyepelasan kuantum Krylov bagi Hamiltonian kekisi
• Primitif Qiskit

Latar belakang

Tutorial ini menunjukkan cara menggunakan penyepelasan kuantum berasaskan sampel (SQD) untuk menganggarkan tenaga keadaan dasar bagi model kekisi fermionic. Khususnya, kami mengkaji model Anderson impuriti tunggal satu dimensi (SIAM), yang digunakan untuk menerangkan impuriti magnetik yang tertanam dalam logam.

Tutorial ini mengikuti aliran kerja yang serupa dengan tutorial berkaitan Penyepelasan kuantum berasaskan sampel bagi Hamiltonian kimia. Namun, perbezaan utama terletak pada cara litar kuantum dibina. Tutorial yang lain menggunakan ansatz variasi heuristik, yang menarik untuk Hamiltonian kimia dengan berpotensi jutaan sebutan interaksi. Sebaliknya, tutorial ini menggunakan litar yang menghampiri evolusi masa oleh Hamiltonian. Litar sedemikian boleh menjadi sangat dalam, menjadikan pendekatan ini lebih sesuai untuk aplikasi pada model kekisi. Vektor keadaan yang disediakan oleh litar-litar ini membentuk asas bagi subruang Krylov, dan hasilnya, algoritma terbukti dan cekap menumpu ke keadaan dasar, di bawah andaian yang sesuai.

Pendekatan yang digunakan dalam tutorial ini boleh dilihat sebagai gabungan teknik yang digunakan dalam SQD dan penyepelasan kuantum Krylov (KQD). Pendekatan gabungan ini kadang-kadang dirujuk sebagai penyepelasan kuantum Krylov berasaskan sampel (SQKD). Lihat Penyepelasan kuantum Krylov bagi Hamiltonian kekisi untuk tutorial tentang kaedah KQD.

Tutorial ini berdasarkan kerja "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", yang boleh dirujuk untuk maklumat lanjut.

Model Anderson impuriti tunggal (SIAM)

Hamiltonian SIAM satu dimensi ialah jumlah tiga sebutan:

$$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}$ ialah operator penciptaan/pemusnahan fermionic bagi tapak mandi ke-$\mathbf{j}^{\textrm{th}}$ dengan spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ ialah operator penciptaan/pemusnahan bagi mod impuriti, dan $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$, dan $V$ ialah nombor nyata yang menerangkan interaksi lompatan, tapak, hibridasi, dan $\varepsilon$ ialah nombor nyata yang menentukan potensi kimia.

Perhatikan bahawa Hamiltonian ini ialah contoh khusus bagi 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 daripada sebutan satu badan, yang bersifat kuadratik dalam operator penciptaan dan pemusnahan fermionic, dan $H_2$ terdiri daripada sebutan dua badan, yang bersifat kuartik. Untuk SIAM,

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

dan $H_1$ mengandungi selebihnya sebutan dalam Hamiltonian. Bagi mewakili Hamiltonian secara aturcara, kami menyimpan matriks $h_{pq}$ dan tensor $h_{pqrs}$.

Asas kedudukan dan momentum

Disebabkan simetri translasi hampir dalam $H_\textrm{bath}$, kami tidak menjangkakan keadaan dasar akan jarang dalam asas kedudukan (asas orbital di mana Hamiltonian ditentukan di atas). Prestasi SQD dijamin hanya jika keadaan dasar adalah jarang, iaitu ia mempunyai berat ketara hanya pada sebilangan kecil keadaan asas pengiraan. Untuk meningkatkan kejarangan keadaan dasar, kami menjalankan simulasi dalam asas orbital di mana $H_\textrm{bath}$ adalah pepenjuru. Kami menyebut asas ini sebagai asas momentum. Kerana $H_\textrm{bath}$ ialah Hamiltonian fermionic kuadratik, ia boleh dipesongjuru dengan cekap melalui putaran orbital.

Evolusi masa hampiran oleh Hamiltonian

Untuk menghampiri evolusi masa oleh Hamiltonian, kami menggunakan penguraian Trotter-Suzuki peringkat 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 masa oleh $H_2$ bersamaan dengan satu gate CPhase antara orbital spin-atas dan spin-bawah pada tapak impuriti. Kerana $H_1$ ialah Hamiltonian fermionic kuadratik, evolusi masa oleh $H_1$ bersamaan dengan putaran orbital.

Keadaan asas Krylov $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, di mana $D$ ialah dimensi subruang Krylov, dibentuk melalui penggunaan berulang satu langkah Trotter, jadi

$$|\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 aliran kerja berasaskan SQD berikut, kami akan mengambil sampel daripada set litar ini dan memproses gabungan set rentetan bit dengan SQD. Pendekatan ini berbeza dengan yang digunakan dalam tutorial berkaitan Penyepelasan kuantum berasaskan sampel bagi Hamiltonian kimia, di mana sampel diambil daripada satu litar variasi heuristik.

Keperluan

Sebelum memulakan tutorial ini, pastikan anda telah memasang yang berikut:

• Qiskit SDK v1.0 atau lebih baru, dengan sokongan 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 berskala kecil

Langkah 1: Petakan masalah ke Circuit kuantum

Pertama, kami menjana Hamiltonian SIAM dalam asas kedudukan. Hamiltonian diwakili oleh matriks $h_{pq}$ dan tensor $h_{pqrs}$. Kemudian, kami memutarnya ke asas momentum. Dalam asas kedudukan, kami meletakkan impuriti di tapak pertama. Namun, apabila kami berputar ke asas momentum, kami memindahkan impuriti ke tapak 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())

Seterusnya, kami menjana litar untuk menghasilkan keadaan asas Krylov. Bagi setiap spesies spin, keadaan awal $\ket{\psi_0}$ diberikan oleh superposisi semua kemungkinan pengujaan tiga elektron yang paling hampir dengan paras Fermi ke dalam 4 mod kosong yang paling hampir bermula dari keadaan $|00\cdots 0011 \cdots 11\rangle$, dan direalisasikan melalui penggunaan tujuh XXPlusYYGate. Keadaan yang telah berevolusi masa dihasilkan melalui penggunaan berturut-turut langkah Trotter peringkat kedua.

Untuk penerangan yang lebih terperinci tentang model ini dan cara litar direka bentuk, rujuk "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)

Langkah 2: Optimumkan masalah untuk pelaksanaan kuantum

Seterusnya, kami mengoptimumkan litar untuk perkakasan sasaran. Buat masa ini, kami akan mencipta Backend generik dengan bilangan qubit yang ditentukan dan set gate yang boleh didekomposkan secara semula jadi oleh litar evolusi masa.

from qiskit.providers.fake_provider import GenericBackendV2

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

Kini, kami menggunakan Qiskit untuk menukar litar kepada Backend sasaran.

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)

Langkah 3: Jalankan menggunakan primitif Qiskit

Selepas mengoptimumkan litar untuk pelaksanaan perkakasan, kami sedia untuk menjalankannya pada perkakasan sasaran dan mengumpul sampel untuk anggaran tenaga keadaan dasar. Selepas menggunakan primitif Sampler untuk mengambil sampel rentetan bit daripada setiap litar, kami menggabungkan semua keputusan ke dalam satu kamus kiraan dan memplot 20 rentetan bit yang paling kerap disampling.

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

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

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 4: Pasca-proses dan kembalikan keputusan ke format klasik yang dikehendaki

Kini, kami menjalankan algoritma SQD menggunakan fungsi diagonalize_fermionic_hamiltonian. Lihat dokumentasi API untuk penjelasan argumen fungsi ini.

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

Sel kod berikut memplot keputusan. Plot pertama menunjukkan tenaga yang dikira sebagai fungsi bilangan iterasi pemulihan konfigurasi, dan plot kedua menunjukkan purata penghunian setiap orbital ruang selepas iterasi akhir. Memandangkan ini adalah masalah yang sangat kecil, iterasi pertama sudah membawa kita sangat hampir dengan tenaga tepat (perhatikan skala paksi y).

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

Sahkan tenaga

Tenaga yang dikembalikan oleh SQD dijamin sebagai had atas kepada tenaga keadaan dasar sebenar. Nilai tenaga boleh disahkan kerana SQD juga mengembalikan pekali vektor keadaan yang menghampiri keadaan dasar. Anda boleh mengira tenaga daripada vektor keadaan menggunakan matriks ketumpatan tereduksi satu- dan dua-zarah, seperti yang ditunjukkan dalam sel kod berikut.

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

Contoh perkakasan berskala besar

Kini, kami menjalankan contoh yang lebih besar pada QPU sebenar. Untuk tenaga rujukan, kami menggunakan keputusan pengiraan DMRG yang dilakukan secara berasingan.

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 seterusnya

Cadangan

Jika anda mendapati karya ini menarik, anda mungkin berminat dengan bahan berikut:

• Penyepelasan kuantum berasaskan sampel bagi Hamiltonian kimia - tutorial berkaitan yang menggunakan ansatz variasi heuristik berbanding litar Trotter
• Penyepelasan kuantum Krylov bagi Hamiltonian kekisi - tutorial tentang kaedah KQD
• Dokumentasi API addon SQD - rujukan untuk fungsi diagonalize_fermionic_hamiltonian
• Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization - kertas yang tutorial ini berdasarkan padanya