Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia

Estimasi penggunaan: kurang dari satu menit pada prosesor Heron r2 (CATATAN: Ini hanya estimasi. Waktu aktual bisa berbeda.)

Hasil pembelajaran

• Cara menggunakan SQD Qiskit addon untuk mengaproksimasi energi ground state sistem molekuler menggunakan bitstring yang disampling dari unit pemrosesan kuantum (QPU).

• Cara menggunakan ffsim untuk membangun Circuit local unitary cluster Jastrow (LUCJ) untuk simulasi kimia kuantum.

Prasyarat

• Kimia kuantum dan second quantization

• Menggunakan Sampler primitive untuk mengambil sampel dari Circuit kuantum

Latar Belakang

Dalam tutorial ini, kita akan menunjukkan cara memproses sampel kuantum yang berisik untuk mengaproksimasi ground state molekul nitrogen $\text{N}_2$ pada panjang ikatan kesetimbangan, dengan menggunakan SQD Qiskit addon untuk mengimplementasikan algoritma diagonalisasi kuantum berbasis sampel (SQD). Detail lebih lanjut tentang perangkat lunaknya dapat ditemukan di dokumentasi yang sesuai, termasuk contoh sederhana untuk memulai.

Kimia nitrogen menghubungkan contoh ini dengan isu keberlanjutan melalui produksi amonia untuk pupuk. Memahami bagaimana katalis mengaktifkan nitrogen adalah bagian dari penelitian tentang metode produksi yang lebih hemat energi. Menghitung ground state dari molekul nitrogen yang terisolasi memperkenalkan metode struktur elektronik yang relevan dengan tantangan yang lebih luas tersebut.

Tutorial ini direkomendasikan bagi pengguna yang sudah familiar dengan kimia kuantum: khususnya, keakraban dalam menemukan energi ground state suatu molekul. Untuk penjelasan mendetail tentang alur kerja, lihat kursus algoritma diagonalisasi kuantum.

SQD adalah teknik untuk menemukan nilai eigen dan vektor eigen dari operator kuantum, seperti Hamiltonian sistem kuantum, dengan menggunakan komputasi kuantum dan komputasi klasik terdistribusi bersama-sama. Komputasi klasik terdistribusi digunakan untuk memproses sampel yang diperoleh dari prosesor kuantum, serta untuk memproyeksikan dan mendiagonalisasi Hamiltonian target dalam subruang yang dibentangkan oleh sampel tersebut. Alur kerja berbasis SQD memiliki langkah-langkah berikut:

1. Pilih ansatz sirkuit dan terapkan pada komputer kuantum ke sebuah state referensi (dalam kasus ini, state Hartree-Fock).

2. Sampel bitstring dari state kuantum yang dihasilkan.

3. Jalankan prosedur self-consistent configuration recovery pada bitstring untuk mendapatkan aproksimasi ke ground state.

SQD diketahui bekerja dengan baik ketika eigenstate target bersifat jarang: fungsi gelombang didukung dalam sekumpulan state basis $\mathcal{S} = \{|x\rangle \}$ yang ukurannya tidak bertambah secara eksponensial seiring dengan ukuran masalah.

Kimia kuantum

Hamiltonian sistem molekuler dapat ditulis sebagai

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},$$

di mana $h_{pr}$ dan $h_{prqs}$ adalah bilangan kompleks yang disebut integral molekuler yang dapat dihitung dari spesifikasi molekul menggunakan sebuah program komputer. Dalam tutorial ini, kita menghitung integral menggunakan paket perangkat lunak PySCF.

Untuk detail tentang cara Hamiltonian molekuler diturunkan, konsultasikan buku teks kimia kuantum (misalnya, Modern Quantum Chemistry oleh Szabo dan Ostlund). Untuk penjelasan tingkat tinggi tentang bagaimana masalah kimia kuantum dipetakan ke komputer kuantum, lihat kuliah Mapping Problems to qubits dari Qiskit Global Summer School 2024.

Ansatz local unitary cluster Jastrow (LUCJ)

SQD memerlukan ansatz sirkuit kuantum untuk mengambil sampel. Dalam tutorial ini, kita akan menggunakan ansatz local unitary cluster Jastrow (LUCJ) karena kombinasinya antara motivasi fisik dan keramahan terhadap perangkat keras. Kita akan menggunakan ffsim untuk membangun Circuit ansatz.

Ansatz LUCJ beradaptasi dengan QPU yang memiliki konektivitas qubit terbatas. Orbital-orbital spin dipetakan ke qubit sedemikian rupa sehingga ansatz tidak memerlukan routing dengan gate SWAP. Perangkat keras IBM® memiliki topologi qubit kisi persegi atau kisi heavy-hex. Untuk kisi heavy-hex, kita dapat mengadopsi pola "zig-zag" yang digambarkan di bawah. Dalam pola ini, orbital dengan spin yang sama dipetakan ke qubit dengan topologi garis (lingkaran merah dan biru), dan terdapat koneksi antara orbital dengan spin berbeda pada setiap orbital spasial ke-4, dengan koneksi yang difasilitasi oleh qubit ancilla (lingkaran ungu).

Pemulihan konfigurasi yang konsisten sendiri

Prosedur self-consistent configuration recovery dirancang untuk mengekstrak sinyal sebanyak mungkin dari sampel kuantum yang berisik. Karena Hamiltonian molekuler mengkonservasi jumlah partikel dan spin Z, masuk akal untuk memilih ansatz sirkuit yang juga mengkonservasi simetri-simetri ini. Ketika diterapkan ke state Hartree-Fock, state yang dihasilkan memiliki jumlah partikel dan spin Z yang tetap dalam pengaturan tanpa noise. Oleh karena itu, bagian spin-$\alpha$ dan spin-$\beta$ dari setiap bitstring yang disampling dari state ini seharusnya memiliki Hamming weight yang sama seperti pada state Hartree-Fock. Karena adanya noise pada prosesor kuantum saat ini, beberapa bitstring yang diukur akan melanggar sifat ini. Bentuk sederhana dari postseleksi akan membuang bitstring-bitstring ini, tetapi hal itu boros karena bitstring tersebut mungkin masih mengandung sinyal. Prosedur self-consistent recovery berupaya memulihkan sebagian sinyal tersebut dalam post-processing. Prosedur ini bersifat iteratif dan memerlukan estimasi awal rata-rata penghunian setiap orbital dalam ground state sebagai input, yang pertama kali dihitung dari sampel mentah. Prosedur ini dijalankan dalam sebuah loop, dan setiap iterasi memiliki langkah-langkah berikut:

1. Untuk setiap bitstring yang melanggar simetri yang ditentukan, balik bit-nya dengan prosedur probabilistik yang dirancang untuk membawa bitstring lebih dekat ke estimasi saat ini dari rata-rata penghunian orbital, untuk mendapatkan bitstring baru.

2. Kumpulkan semua bitstring lama dan baru yang memenuhi simetri, lalu bagi subsampel dari bitstring tersebut menjadi batch terpisah dengan ukuran tetap yang dipilih sebelumnya.

3. Untuk setiap batch bitstring, proyeksikan Hamiltonian ke subruang yang dibentangkan oleh vektor basis yang sesuai (lihat bagian sebelumnya untuk deskripsi vektor basis ini), dan hitung estimasi ground state dari Hamiltonian yang diproyeksikan pada komputer klasik. Setiap bitstring terdiri dari bagian spin-$\alpha$ dan bagian spin-$\beta$, yang masing-masing kita sebut CI string. Subruang untuk sebuah batch dibentangkan oleh vektor basis yang terbentuk dengan memasangkan setiap CI string spin-$\alpha$ dalam batch dengan setiap CI string spin-$\beta$ dalam batch.

4. Pangkas (trim) setiap batch, dengan hanya mempertahankan CI string yang memiliki bobot terbesar dalam estimasi ground state-nya. Gabungkan CI string yang bertahan dari semua batch, lalu diagonalisasi Hamiltonian sekali lagi pada subruang gabungan tersebut.

5. Perbarui estimasi rata-rata penghunian orbital menggunakan estimasi ground state dari subruang gabungan. CI string dengan bobot tertinggi darinya juga dibawa ke subruang pada iterasi berikutnya.

Pemangkasan dan penggabungan batch pada langkah 4 adalah varian trim SQD dari algoritma ini. Karena mendiagonalisasi sebuah batch untuk menyaring CI string-nya lebih murah daripada mendiagonalisasi subruang yang cukup besar untuk menampung semua CI string kandidat, pemangkasan memungkinkan SQD memilih dari kumpulan CI string hasil sampling yang lebih besar daripada yang bisa didiagonalisasi sekaligus.

Diagram alur kerja SQD

Alur kerja SQD digambarkan dalam diagram berikut (pemangkasan dan penggabungan batch yang dijelaskan pada langkah 4 tidak ditampilkan):

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu telah menginstal yang 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.14 atau lebih baru (pip install qiskit-addon-sqd)

• ffsim v0.0.75 atau lebih baru (pip install ffsim)

Pengaturan

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"ffsim": "ffsim", "matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Contoh simulator skala kecil

Dalam tutorial ini, kita akan mencari aproksimasi terhadap ground state molekul nitrogen mendekati jarak ikatan kesetimbangannya. Kita pertama menggunakan basis set STO-6G yang kecil agar bisa mensimulasikan eksperimen dan memastikannya berfungsi.

Langkah 1: Petakan input klasik ke masalah kuantum

Pertama, kita tentukan molekul dan properti-propertinya.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="sto-6g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

converged SCF energy = -108.464957764795

CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000

norb = 8
nelec = (5, 5)

Sebelum membangun Circuit ansatz LUCJ, kita terlebih dahulu melakukan perhitungan CCSD pada sel kode berikut. Amplitudo $t_1$ dan $t_2$ dari perhitungan ini akan digunakan untuk menginisialisasi parameter ansatz.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -108.5933309085007  E_corr = -0.1283731437052356

Sekarang, kita gunakan ffsim untuk membuat Circuit ansatz. Karena molekul kita memiliki state Hartree-Fock closed-shell, kita gunakan varian spin-balanced dari ansatz UCJ, yaitu UCJOpSpinBalanced. Kita set optimize=True pada metode from_t_amplitudes untuk mengaktifkan double-factorization "terkompresi" dari amplitudo $t_2$ (lihat The local unitary cluster Jastrow (LUCJ) ansatz di dokumentasi ffsim untuk detailnya).

Karena ansatz LUCJ beradaptasi dengan konektivitas yang tersedia pada QPU, kita harus menginisialisasi backend QPU sebelum membuat ansatz. Untuk saat ini, kita akan membuat backend generik dengan coupling map kisi persegi dan gate set yang secara alami menjadi hasil dekomposisi ansatz LUCJ. Kemudian, kita akan menggunakan ffsim.qiskit.generate_lucj_pass_manager untuk membuat pass manager yang dikhususkan untuk melakukan transpile ansatz LUCJ ke backend yang diberikan sesuai dengan layout "ladder" (lihat penjelasan di dokumentasi ffsim). Fungsi ini menggunakan heuristik penilaian untuk meminimalkan error yang terkait dengan layout yang dipilih, yang penting jika backend kamu adalah QPU nyata atau simulator dengan model noise. Selain mengembalikan pass manager, fungsi ini juga mengembalikan pasangan interaksi alpha-beta yang dapat diimplementasikan pada hardware. Jika tidak semua pasangan dapat diimplementasikan, fungsi ini mengeluarkan peringatan.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 4
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
rng = np.random.default_rng(12345)
coupling_map = CouplingMap.from_grid(norb, 2)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
    seed=rng,
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
    seed_transpiler=int(rng.integers(2**31)),
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

Langkah 2: Optimalkan untuk eksekusi hardware kuantum

Selanjutnya, kita optimalkan Circuit untuk hardware target. Biasanya, langkah ini melibatkan inisialisasi backend hardware dan pass manager untuk backend tersebut. Namun, karena ansatz LUCJ sudah disesuaikan dengan konektivitas hardware, kita sudah melakukan tindakan-tindakan tersebut di langkah sebelumnya. Yang perlu dilakukan hanyalah menjalankan pass manager pada Circuit untuk melakukan transpile menjadi ISA circuit yang dapat langsung dieksekusi di QPU.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})

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. Karena kita hanya memiliki satu sirkuit, kita akan menggunakan mode eksekusi Job dari IBM Quantum Compute Service dan mengeksekusi Circuit kita.

sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

primitive_result = job.result()
pub_result = primitive_result[0]

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang diinginkan

Metrik yang berguna untuk menilai kualitas output QPU adalah jumlah konfigurasi valid yang dikembalikan. Konfigurasi valid memiliki jumlah partikel dan spin Z yang benar, yang berarti separuh kanan dari bitstring memiliki Hamming weight yang sama dengan jumlah elektron spin-up, dan separuh kiri memiliki Hamming weight yang sama dengan jumlah elektron spin-down. Sel berikut menghitung fraksi konfigurasi yang disampling yang valid.

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(
    f"Fraction of sampled configurations that are valid (simulation): {valid_fraction}"
)

Fraction of sampled configurations that are valid (simulation): 1.0

Semua bitstring valid karena kita mengambil sampel Circuit pada simulator tanpa noise. Ketika menjalankan pada QPU yang berisik, fraksi tersebut akan kurang dari satu, tetapi diharapkan lebih besar dari fraksi yang akan didapat jika bitstring disampling secara acak seragam, yang dihitung pada sel berikut.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Sekarang, kita estimasi energi ground state dari Hamiltonian menggunakan fungsi diagonalize_fermionic_hamiltonian. Fungsi ini melakukan prosedur configuration recovery yang self-consistent untuk secara iteratif menyempurnakan sampel kuantum yang berisik guna meningkatkan estimasi energi. Kita berikan fungsi callback agar bisa menyimpan hasil antara untuk analisis lebih lanjut. Lihat dokumentasi API untuk penjelasan argumen-argumen diagonalize_fermionic_hamiltonian.

Di sini, kita menggunakan argumen initial_occupancies pada diagonalize_fermionic_hamiltonian untuk menentukan konfigurasi Hartree-Fock sebagai tebakan awal untuk penghunian orbital dalam ground state. Pendekatan ini masuk akal untuk sistem di mana ground state memiliki dukungan signifikan pada konfigurasi Hartree-Fock, tetapi mungkin tidak cocok dalam situasi lain, meskipun metode komputasi yang lebih canggih mungkin menghasilkan tebakan awal yang lebih baik dalam kasus-kasus tersebut. Menentukan initial_occupancies juga memungkinkan configuration recovery berjalan bahkan jika tidak ada konfigurasi valid yang disampling, seperti yang mungkin terjadi saat mengambil sampel dari Circuit besar pada QPU yang berisik. Tanpa argumen ini, configuration recovery akan gagal dan menampilkan error jika tidak ada konfigurasi valid yang diberikan.

Kita mengaktifkan trim SQD dengan meneruskan TrimPolicy sebagai argumen policy. Setiap iterasi kemudian melakukan num_batches + 1 diagonalisasi: satu untuk setiap batch, dan satu untuk subruang gabungan. Kita menggunakan dua batch, sehingga setiap iterasi memerlukan tiga diagonalisasi. Callback hanya menerima hasil gabungan, sehingga setiap iterasi mencetak satu energi. Opsi TrimPolicy adalah:

• trim_ratio: fraksi CI string setiap batch, per sektor spin, yang bertahan setelah pemangkasan dan masuk ke subruang gabungan. Subruang gabungan berisi hingga sekitar num_batches * trim_ratio kali lipat jumlah CI string sebuah batch tunggal, dikurangi CI string yang dimiliki bersama oleh batch-batch tersebut. Argumen max_dim membatasi setiap subruang, termasuk subruang gabungan.
• carryover_threshold dan carryover_ratio: CI string mana dari subruang gabungan yang dibawa ke iterasi berikutnya. String harus muncul dalam konfigurasi yang amplitudonya melebihi carryover_threshold, dan dari string tersebut, fraksi carryover_ratio dengan bobot terbesar dipertahankan.

Sistem ini cukup kecil sehingga satu batch saja sudah memuat sebagian besar CI string yang mungkin, sehingga batch-batch saling tumpang tindih secara signifikan. Karena itu, kita menggunakan trim_ratio yang tinggi agar subruang gabungan tidak menjadi jauh lebih kecil daripada batch. Lihat panduan trim SQD untuk detail lebih lanjut tentang cara memilih opsi-opsi ini.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)
from qiskit_addon_sqd.trim import TrimPolicy

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5
samples_per_batch = 1000
num_batches = 2
max_dim = 50
policy = TrimPolicy(
    trim_ratio=0.8,
    carryover_threshold=1e-4,
    carryover_ratio=0.5,
)
symmetrize_spin = True

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)

# 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"\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    policy=policy,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy (simulation): {final_energy}")
print(f"Final energy error (simulation): {energy_error}")

Iteration 1
	Energy: -108.59585971421083
	Subspace dimension: 1681

Iteration 2
	Energy: -108.59585971421083
	Subspace dimension: 1681
Final energy (simulation): -108.59585971421083
Final energy error (simulation): 0.00012763677520410965

Visualisasikan hasilnya

Plot pertama menunjukkan error energi dari subruang tergabung di setiap iterasi. Dalam simulasi ini, kita sudah berada dalam rentang 1 mH dari jawaban eksak setelah iterasi pertama (akurasi kimia umumnya diterima sebagai 1 kcal/mol $\approx$ 1.6 mH). Namun ini adalah sistem kecil, dan karena sampelnya tanpa noise, configuration recovery tidak diperlukan. Pada sistem yang lebih besar yang dijalankan pada QPU yang berisik, beberapa iterasi configuration recovery mungkin diperlukan, dan akurasi akhir mungkin lebih buruk. Secara umum, energi dapat ditingkatkan dengan mengizinkan lebih banyak iterasi configuration recovery atau dengan menambah jumlah sampel per batch.

Plot kedua menunjukkan rata-rata occupancy setiap orbital spasial setelah iterasi terakhir. Kita bisa melihat bahwa elektron spin-up maupun spin-down menempati lima orbital pertama dengan probabilitas tinggi dalam solusi kita.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# 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, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-5, 1e-2)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", 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})

plt.tight_layout()
plt.show()

Contoh hardware skala besar

Sekarang, kita menjalankan contoh yang lebih besar pada hardware kuantum nyata. Di sini, kita akan menurunkan active space untuk molekul nitrogen dari basis set cc-pVDZ. Kita juga akan membuat perubahan berikut:

• Turunkan n_reps menjadi 1 untuk mengurangi noise dari error gate
• Atur konektivitas ansatz berdasarkan keluarga prosesor backend (heavy-hex untuk Heron, persegi untuk Nighthawk)
• Tingkatkan max_dim menjadi 1000 untuk memungkinkan dimensi subruang yang lebih besar.
• Turunkan trim_ratio menjadi 0.5. Sistem ini memiliki jauh lebih banyak CI string yang mungkin daripada yang muat dalam satu subruang, jadi dengan dua batch, rasio 0.5 menjaga subruang gabungan tidak lebih besar daripada satu batch.

Langkah 1-4

Di sini kita merangkum semua langkah menjadi satu alur kerja tunggal dalam skala yang lebih besar, yang kemudian dijalankan pada hardware kuantum nyata.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="cc-pvdz",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

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

# Create pass manager
connectivity_map = {"Heron": "heavy-hex", "Nighthawk": "square"}
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity_map[
        backend.configuration().processor_type["family"]
    ],
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(
    f"Fraction of sampled configurations that are valid (hardware): {valid_fraction}"
)
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5
samples_per_batch = 1000
num_batches = 2
max_dim = 1000
policy = TrimPolicy(
    trim_ratio=0.5,
    carryover_threshold=1e-4,
    carryover_ratio=0.5,
)
symmetrize_spin = True

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the
# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)

# List to capture intermediate results
result_history = []

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    policy=policy,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy (hardware): {final_energy}")
print(f"Final energy error (hardware): {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# 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, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-3, 1e-1)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", 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})

plt.tight_layout()
plt.show()

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450045
Using backend ibm_kingston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7039, 'rz': 6977, 'cz': 1858, 'measure': 52, 'x': 45, 'barrier': 1})
Fraction of sampled configurations that are valid (hardware): 0.00248
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Energy: -109.19178247530117
	Subspace dimension: 644809
Iteration 2
	Energy: -109.20055015645298
	Subspace dimension: 339889
Iteration 3
	Energy: -109.20621939005301
	Subspace dimension: 446224
Iteration 4
	Energy: -109.21228961696536
	Subspace dimension: 393129
Iteration 5
	Energy: -109.21328434918014
	Subspace dimension: 403225
Final energy (hardware): -109.21328434918014
Final energy error (hardware): 0.014744867477020307

Langkah selanjutnya

Rekomendasi

Jika kamu merasa topik ini menarik, kamu mungkin tertarik dengan materi berikut:

• Diagonalisasi kuantum Krylov berbasis sampel untuk model kisi fermionic - tutorial terkait yang menggunakan Circuit evolusi waktu alih-alih ansatz variasional
• Skalakan alur kerja kimia SQD dengan pemecah Dice - halaman yang menunjukkan cara menggunakan perangkat lunak Dice yang lebih efisien untuk diagonalisasi
• Dokumentasi API addon SQD - referensi untuk fungsi diagonalize_fermionic_hamiltonian
• Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer - makalah yang menjadi dasar tutorial ini