Algoritma SqDRIFT untuk estimasi keadaan dasar

Perkiraan penggunaan: 180 detik pada prosesor Heron r3 (CATATAN: Ini hanya perkiraan. Waktu jalanmu mungkin berbeda.)

Hasil pembelajaran

• Pelajari cara membuat Circuit dengan kedalaman lebih kecil dibandingkan Trotterisasi

• Telusuri workflow end-to-end untuk estimasi keadaan dasar memakai qDRIFT dan SQD

• Pelajari cara memakai qiskit-fermions bersama addon Qiskit lainnya untuk mengimplementasikan workflow seperti itu

Tutorial ini disajikan sebagai notebook Python untuk tujuan pengajaran.

Prasyarat

• Baca ikhtisar Sample-based quantum diagonalization (SQD)

• Baca pelajaran Sample-based Krylov Quantum Diagonalization (SKQD)

Latar belakang

SqDRIFT adalah varian SKQD yang menggantikan keharusan memilih ansatz untuk men-sampel bitstring dengan ensemble Circuit evolusi-waktu yang dibangun langsung dari Hamiltonian target. Ini dicapai dengan men-subsampel operator evolusi-waktu yang lebih kecil dari Hamiltonian berdasarkan koefisiennya, yang dikenal sebagai metode Trotterisasi qDRIFT.

Tutorial ini memakai Qiskit Fermions untuk membuat Circuit fermionik yang lebih natural untuk algoritma qDRIFT, lalu memakai pass layout dan sintesis fermionik sebelum memasukkan Circuit ke pipeline Qiskit tradisional untuk eksekusi di hardware.

Misalkan Hamiltonian berbentuk:

$$H = \sum_{i=1}^{N} c_i h_i$$

di mana, tanpa mengurangi keumuman, kita mensyaratkan $c_i > 0$ dan eigenvalue terbesar dari $h_i$ sama, dalam nilai absolut, dengan $1$. Prafaktor bertanda atau kompleks apa pun diserap ke dalam $h_i$, sehingga koefisien $c_i$ adalah bobot yang benar-benar positif sementara $h_i$ membawa arah tiap suku. Di sini $N$ adalah jumlah suku (atau, setelah pengelompokan, jumlah kelompok) dalam Hamiltonian; ia adalah properti Hamiltonian dan berbeda dari jumlah operator yang disampel ke dalam satu Circuit, ditulis $n$ di bawah.

Algoritma qDRIFT kemudian merealisasikan, untuk waktu target $t$, suatu operator $V_k$, di mana $k$ berjalan dari $1 \cdots K$ dan menandakan Circuit SqDRIFT ke-$k_{th}$, didefinisikan sebagai:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Di sini $n$ adalah jumlah operator yang disampel per Circuit dan $K$ adalah jumlah Circuit dalam ensemble. Perkaliannya berjalan atas $n$ pengundian, bukan atas semua $N$ suku Hamiltonian, dan karena suku diundi dengan pengembalian, $h_i$ yang sama bisa muncul lebih dari sekali dalam satu $V_k$.

Besaran:

$$\lambda = \sum_{i=1}^{N} c_i$$

adalah norma $L_1$ dari koefisien, sehingga masing-masing dari $n$ langkah berevolusi selama durasi yang sama $\lambda t / n$ terlepas dari suku mana yang diundi. Keseragaman sudut langkah adalah ciri khas qDRIFT: sebuah koefisien memengaruhi hasil lewat seberapa sering sukunya diundi, bukan lewat seberapa jauh suku itu dirotasi. Indeks-indeksnya disampel dari distribusi:

$$P[k_i] = \frac{c_i}{\lambda}$$

jadi deret $(k_1, \ldots, k_n)$ adalah urutan acak indeks suku yang diundi dari distribusi ini. Karena $c_i$ positif dan berjumlah $\lambda$, ini adalah distribusi probabilitas ternormalisasi, dan ekspektasi kanal yang dihasilkan atas undian acak mendekati evolusi di bawah $H$, dengan galat yang menurun seiring $n$ bertambah. Perhatikan bahwa galat aproksimasi bergantung pada $\lambda$, bukan pada jumlah suku $N$.

(Paper SqDRIFT menulis jumlah suku sebagai $\mathcal{N}$ dan panjang urutan sebagai $N$; kita memakai $N$ dan $n$ di sini agar keduanya tetap jelas berbeda.)

Tutorial ini menunjukkan cara membuat ensemble Circuit acak seperti itu. Setelah Circuit-Circuit ini dibuat, mirip cara kita membuat subruang Krylov untuk operator yang berbeda, kita men-sampel bitstring dari beberapa operator seperti itu dengan parameter waktu yang berbeda. Ini memastikan overlap yang lebih tinggi antara vektor keadaan dasar dan bitstring yang disampel.

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu sudah menginstal

• Virtual environment Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Perhatikan bahwa namanya jamak)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Kamu bisa menginstal semua paket yang dibutuhkan dengan:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Setup

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

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "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")

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Contoh simulator

Step 1: Petakan input klasik ke masalah kuantum

Membaca dan Menyiapkan FCIDump

Untuk tutorial ini, kita akan memuat Hamiltonian struktur elektronik untuk nitrogen (N2). Ada cara lain juga untuk membuat operator fermionik. Lihat dokumentasi di qiskit_fermions.operators.library.

Tentang FCIDump ini. File N2_sto_3g menggambarkan molekul nitrogen ($N_2$) dalam basis STO-3G minimal pada jarak antaratom 1.09 $\AA$, panjang ikatan kesetimbangan eksperimental. Headernya menyatakan NORB=10, NELEC=14, dan MS2=0: 10 orbital spasial (jadi 20 orbital spin, dan 20 Qubit di bawah Jordan-Wigner), 14 elektron dalam singlet spin, jadi tujuh elektron $\alpha$ dan tujuh $\beta$. Semua orbital diberi label simetri 1, artinya tidak ada simetri grup-titik yang dimanfaatkan. Karena ini dump STO-3G ruang penuh, tidak ada orbital yang dibekukan dan ruang korelasinya cukup kecil sehingga energi referensi FCI eksak bisa dihitung secara klasik sebagai pembanding, seperti ditunjukkan di sel berikutnya.

File yang setara bisa dibuat ulang dengan PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Karena integralnya bergantung pada orbital SCF yang konvergen, file hasil buat ulang bisa berbeda dari yang disertakan dalam fase atau urutan orbital; energi totalnya tidak terpengaruh.

Mendapatkan file. Temukan FCIDump di repositori GitHub ini. Kamu bisa menjalankan sel di bawah untuk mengambilnya ke lokasi yang diharapkan oleh bagian tutorial lainnya.

Pertama kita pakai cisolver dari pyscf untuk mendapatkan energi referensi. Ini adalah energi keadaan dasar sebenarnya dari molekul yang kita pakai. Untuk itu kita deklarasikan dulu norb dan nelec, yaitu jumlah orbital dan jumlah elektron. Lalu kita deklarasikan h1e dan h2e, yaitu integral satu- dan dua-elektron. Semuanya nanti juga akan dipakai untuk SQD.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Memuat Hamiltonian

Dengan data yang dibutuhkan sudah siap, kita baca Hamiltonian dari file FCI dalam format yang kompatibel dengan qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Workflow fermionik dengan qiskit-fermions

Pertama kita petakan Hamiltonian ke model Circuit fermionik memakai qiskit-fermions, yang menyediakan pass Transpiler dan Gate khusus untuk Circuit fermionik. Ini nanti dipakai sebelum pass Transpiler tradisional Qiskit untuk workflow ini.

Pengelompokan suku

Untuk memastikan reproduksibilitas hasil, pertama kita pakai canonical_order untuk mengurutkan suku hanya berdasarkan strukturnya. Urutan operator dalam daftar canon dengan demikian tetap. Ini menjamin reproduksibilitas operator yang dibuat karena pass QDriftTrotterization yang akan kita pakai nanti men-sampel indeks acak untuk membuat operator qDRIFT.

Pada langkah ini, kita memanfaatkan banyak simetri yang ada dalam Hamiltonian struktur elektronik dengan mengelompokkan suku terkait yang berkoefisien identik. Meski ini mengubah distribusi koefisien operator yang disampel protokol qDRIFT, hal ini tidak memengaruhi jaminan konvergensinya. Yang penting, mengelompokkan suku yang terkait simetri menghasilkan pembatalan suku Pauli yang menguntungkan dan kedalaman Circuit keseluruhan yang lebih pendek saat mengevolusi-waktu sebuah keadaan di bawah aksinya.

qiskit-fermions menyediakan fungsi group_terms_by_electronic_structure yang melakukan pengelompokan ini untuk kita.

Perhatikan bahwa group_terms_by_electronic_structure mengasumsikan suku yang normal ordering.

Menyaring suku diagonal

Kita menghapus suku diagonal dari Hamiltonian yang dipakai untuk membuat Circuit, supaya $n$ slot sampling qDRIFT dipakai untuk suku yang memindahkan populasi antar konfigurasi. Suku seperti itu sebaiknya disaring dari Hamiltonian pada titik ini, sebelum Gate Evolution dibangun di langkah berikutnya.

Suku yang dimaksud adalah yang diagonal dalam basis bilangan-okupansi, yaitu perkalian operator bilangan $a^\dagger_i a_i$. Tiga jenis suku masuk dalam deskripsi ini:

• offset energi konstan, perkalian nol operator bilangan, yang evolusi waktunya hanya menyumbang fase global;

• operator bilangan individual $n_i$, yang evolusi waktunya menyusut menjadi rotasi $Z$ Qubit-tunggal;

• perkalian orde lebih tinggi seperti $n_i n_j$.

Dengan sendirinya, tidak satu pun dari ini memindahkan populasi antar konfigurasi bilangan-okupansi; mereka hanya bekerja pada fase konfigurasi yang sudah ada. Namun mereka tidak inert: fase relatif itu masuk ke interferensi yang dihasilkan oleh suku eksitasi di bagian Circuit berikutnya, jadi menyaringnya mengubah evolusi yang sebenarnya dihasilkan dan bisa mengubah distribusi sampling. Ini adalah aproksimasi yang disengaja dalam langkah pembuatan Circuit, dibuat untuk memfokuskan sampling pada suku eksitasi, bukan langkah yang membiarkan distribusi tersampel tidak tersentuh. Berbeda dengan pengelompokan simetri di atas, yang menjaga jaminan konvergensi qDRIFT tetap utuh, penyaring ini mengubah operator yang dievolusikan. Circuit dengan demikian tidak lagi mengaproksimasi evolusi di bawah Hamiltonian penuh, dan batas galat qDRIFT berlaku untuk operator yang telah disaring, bukan yang asli. Ini dapat diterima di sini karena Circuit hanyalah heuristik sampling yang dipakai untuk mengusulkan konfigurasi: tidak ada suku yang hilang dari estimasi energi itu sendiri, karena penyaring hanya berlaku untuk Hamiltonian yang dipakai membangun Circuit, sedangkan diagonalisasi klasik nanti memakai Hamiltonian penuh, termasuk suku diagonal. Akurasi SQD bergantung pada langkah klasik itu, yang tetap variasional dalam subruang yang disampel terlepas dari bagaimana konfigurasi diusulkan.

Fungsi filter_diagonal_terms() menghapus suku seperti itu dari sebuah operator di tempat. Ia mengidentifikasinya dari struktur normal-ordered-nya — multiset mode kreasi yang cocok dengan multiset mode anihilasi — jadi hanya valid pada operator yang sudah normal-ordered. Asumsi ini tidak dicek saat runtime.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Setelah mengelompokkan suku dalam Hamiltonian, kita akan menentukan parameter berikut untuk membuat ensemble Circuit:

• Jumlah Circuit yang dibuat: num_circuits
• Panjang setiap Circuit dalam kelompok eksitasi: num_exc
• Faktor untuk waktu evolusi yang berbeda: times

Membuat Circuit fermionik

Sekarang kita buat Circuit fermionik untuk setiap langkah waktu. Setiap Circuit terdiri dari satu Gate evolusi, dengan waktu evolusi yang sudah kita deklarasikan tadi. Operator evolusinya adalah Hamiltonian. Nanti kita jalankan pass Transpiler pada Circuit ini untuk membuat Circuit qDRIFT.

Persiapan ansatz

Kita menyiapkan keadaan Hartree-Fock memakai kelas InitializeModes. Untuk nitrogen, prosesnya sederhana: menerapkan Gate X pada num_elec_a Qubit pertama lalu pada num_elec_b Qubit, keduanya sama dengan tujuh untuk nitrogen. Keadaan ini merepresentasikan tujuh elektron $\alpha$ dan tujuh $\beta$ dari nitrogen.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Step 2: Optimalkan masalah untuk eksekusi di hardware kuantum

Setelah Circuit kita siap, pertama kita pakai pass yang tersedia di qiskit-fermions untuk melakukan optimasi tingkat fermionik, lalu men-transpile Circuit kita untuk Backend pilihan. Karena ini eksperimen simulator, pertama kita lakukan untuk AerSimulator. Perhitungan bobot untuk setiap kelompok

Di langkah ini, kita melakukan sampling qDRIFT atas suku-suku secara stokastik dengan probabilitas sebanding dengan koefisiennya di Hamiltonian. Pass Transpiler qDRIFT yang mengerjakannya buat kita. Sekarang kita bisa membuat circuit yang lebih dangkal dan bisa dijalankan di hardware dengan lebih efisien meski konektivitas qubit terbatas, bahkan ketika Hamiltonian punya kopling jarak jauh dan suku berorde lebih tinggi dari kuadratik. Setelah pengelompokan suku, operator di-sampling berdasarkan bobotnya. Untuk setiap operator $h_i$, bobot $W_{h_i}$ didefinisikan sebagai berikut:

$$W_{h_i} = |c_i| / \lambda$$

Karena suku-suku sudah dikelompokkan di Langkah 1, setiap $h_i$ di sini adalah satu grup utuh: $c_i$ adalah rata-rata koefisien absolut dari suku-suku di grup $i$, dan setiap suku di grup itu dievolusikan dengan koefisien yang direduksi menjadi tandanya saja.

Optimasi fermionik dan optimasi native hardware

Fungsi generate_preset_jw_pass_manager() mengembalikan MultiStagePassManager yang menerima FermionicCircuit dan menghasilkan circuit akhir yang sudah dioptimalkan, yang bisa kita transpile untuk dijalankan di hardware kita. Kita mengganti tahap optimasi bawaannya dengan FermionicPassManager yang berisi pass QDriftTrotterization kita:

• Pass QDriftTrotterization memakai penghitungan bobot dan sampling secara internal untuk menghasilkan circuit yang akan kita pakai untuk sampling

• Pass RelabelModes adalah pass optimasi lain yang bisa dipakai untuk mempermutasi mode fermionik guna mengoptimalkan konektivitas antar qubit dan mengurangi kedalaman gate; baca selengkapnya di referensi API

Tahap-tahap sisa dari MultiStagePassManager berjalan otomatis dan menangani seluruh pemetaan fermion ke qubit:

• F2QLayout: Preset pass manager menerapkan pass TrivialF2QLayout, yang memetakan $n$ bit fermionik ke $n$ qubit secara trivial.

• F2QSynth: Pass transpilasi untuk memetakan instruksi circuit berbasis fermion ke instruksi berbasis qubit.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Sekarang optimasi tingkat fermionik sudah selesai, kita bisa mentranspile circuit untuk dijalankan di simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Langkah 3: Jalankan dengan primitif Qiskit

Sekarang circuit sudah siap, kita bisa menjalankannya dengan primitif Qiskit di AerSimulator. Kita akan menggabungkan semua hasil count dari circuit yang berbeda. Kita ubah jadi vektor boolean sebelum akhirnya diproses lanjut dengan SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

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

Memakai bitstring untuk SQD

Sekarang kita bisa menjalankan skema diagonalisasi pada bitstring terpilih untuk menemukan eigenvalue terendah, yang sesuai dengan energi keadaan dasar molekul. Kita buat fungsi callback, deklarasikan okupansi awal, dan atur parameternya sebelum akhirnya menjalankan skema diagonalisasi. Fungsi callback dipakai untuk mencetak iterasi saat ini dan estimasi eigenvalue saat ini di setiap iterasi.

Terakhir, untuk mendapatkan estimasi keadaan dasar, kita tambahkan nuclear_repulsion_energy ke energi hasilnya.

Catatan: Dimensi subruang tidak tetap antar iterasi, bahkan di simulator tanpa noise — setiap subsampel menarik konfigurasi yang berbeda, dan langkah pemulihan membentuk ulang pool di antara iterasi, jadi dimensi yang dilaporkan bervariasi dari satu subsampel ke subsampel berikutnya. Sampling tanpa noise saja tidak otomatis mengunci dimensi subruang terpilih. Namun, run di hardware cenderung menghasilkan subruang yang secara sistematis lebih besar, karena shot yang noisy melanggar simetri jumlah partikel dan pemulihan konfigurasi mengubahnya menjadi vektor basis tambahan. Karena itu, di bagian hardware kita juga akan memperkenalkan satu langkah lagi untuk memangkas bitstring.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Contoh hardware

Contoh ini memakai 20 qubit (10 orbital spasial). Pilihan itu cuma demi kemudahan untuk tutorial yang harus berjalan cepat, bukan batas keras bagi metodenya.

Biaya langkah klasik tidak ditentukan langsung oleh jumlah qubit. SQD mendiagonalisasi Hamiltonian yang diproyeksikan ke subruang yang direntang oleh konfigurasi hasil sampling, jadi yang menentukan biaya klasik adalah dimensi subruang terpilih itu — di sini diatur oleh samples_per_batch, num_batches, dan berapa banyak konfigurasi berbeda yang benar-benar dihasilkan circuit — bersama aljabar linear sparse yang dibutuhkan untuk menerapkan Hamiltonian terproyeksi. Ruang CI penuh tumbuh secara kombinatorial seiring orbital dan elektron, tetapi subruang terpilih hanyalah irisan kecil yang bisa diatur darinya, dan ukurannya kita kendalikan langsung. Akibatnya, jumlah qubit dan kesulitan klasik bisa divariasikan agak independen: ruang orbital yang lebih lebar yang di-sampling ke subruang sederhana bisa lebih murah daripada sistem yang lebih kecil yang didiagonalisasi di subruang yang sangat besar.

Jadi dalam praktiknya, ukuran sistem yang layak bergantung pada dimensi subruang yang kamu butuhkan untuk akurasi yang diinginkan, serta memori dan core yang tersedia untuk eigensolver. Ruang orbital yang lebih besar biasanya memang butuh subruang lebih besar untuk mencapai akurasi kimia, dan itulah yang akhirnya mendorong penggunaan sumber daya terdistribusi — lihat qiskit-addon-sqd-hpc untuk menskalakan langkah ini. Daripada berasumsi ada batas tetap, pendekatan praktisnya adalah memantau dimensi subruang yang dilaporkan dan konvergensi energi antar iterasi, lalu memperbesar ukuran subruang sampai energi berhenti membaik atau memori yang tersedia habis.

Catatan: Karena galat sampling akibat noise di hardware, subruang yang dibuat untuk diagonalisasi pada run hardware akan lebih besar daripada yang kita dapat saat memakai simulator. Meski ini menambah dimensi subruang yang ingin kita diagonalisasi, workflow tetap memberi jawaban akurat berkat ketahanan SQD terhadap noise.

Pemangkasan string palsu

Di sini kita bisa memilih untuk melakukan satu langkah tambahan. Setelah punya semua bitstring dari eksekusi circuit, kita bisa menyaring bitstring yang tidak valid sebelum menjalankan SQD, atau lanjut tanpa pemangkasan. Melewatkan pemangkasan umumnya lebih baik untuk run hardware, karena shot yang simetrinya rusak tetap tersedia untuk pemulihan konfigurasi, yang bisa memperbaikinya menjadi konfigurasi valid dan dengan begitu memperlebar subruang, alih-alih langsung membuang shot tersebut.

Karena nitrogen hanya bisa punya tujuh elektron $\alpha$ dan tujuh elektron $\beta$, bitstring apa pun yang punya lebih atau kurang dari tujuh angka 1 di paruh pertama dan kedua output bisa dibuang. Kita definisikan fungsi yang memeriksa apakah bitstring valid, dan kalau tidak, membuangnya. Setelah bitstring palsu disaring, sisanya dikirim ke skema diagonalisasi. Pakai flag PRUNE di bawah untuk berpindah antara dua perilaku itu.

Perlu diingat bahwa pemangkasan hanyalah satu dari beberapa pilihan yang membentuk subruang akhir, selain jumlah circuit, himpunan waktu evolusi, dan penyaringan suku diagonal. Membandingkan run dengan pemangkasan dan tanpa pemangkasan baru informatif kalau semua hal lain dibuat tetap; versi C++ dari tutorial ini membahasnya lebih detail, karena ia melakukan postselection (bukan pemulihan) dan juga berbeda di parameter-parameter lain tersebut.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

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
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Langkah berikutnya

Rekomendasi

Kalau karya ini menarik buatmu, kamu mungkin tertarik dengan materi berikut:

• Sample-based Krylov quantum diagonalization of a fermionic lattice model - tutorial terkait yang memakai circuit evolusi waktu alih-alih ansatz variasional.
• Sample-based quantum diagonalization of a chemistry Hamiltonian - tutorial tentang cara membangun circuit local unitary cluster Jastrow (LUCJ) untuk simulasi kimia kuantum.
• Paper SqDRIFT - literatur yang jadi dasar tutorial ini. (Perhatikan bahwa sebagian optimasi yang dibahas di paper ini saat ini masih dalam pengerjaan, dan tutorial ini bisa berubah di masa depan mengikuti perkembangan library yang dipakai.)