Algoritma SqDRIFT untuk estimasi ground state

Estimasi penggunaan: 180 detik pada prosesor Heron r3 (CATATAN: Ini hanya perkiraan. Runtime Anda mungkin berbeda.)

Mencari versi C++?

Tutorial ini menggunakan Python. Untuk implementasi C++, termasuk source code dan instruksi build, lihat tutorial SqDRIFT C++.

Hasil pembelajaran

• Pelajari cara membuat circuit dengan kedalaman lebih kecil dibandingkan Trotterization

• Telusuri workflow end-to-end untuk estimasi ground state menggunakan qDRIFT dan SQD

• Pelajari cara menggunakan qiskit-fermions bersama dengan addon Qiskit lainnya untuk mengimplementasikan workflow semacam itu

Tutorial ini disajikan sebagai notebook Python untuk tujuan pengajaran.

Prasyarat

• Baca ringkasan Sample-based quantum diagonalization (SQD)

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

Latar Belakang

SqDRIFT adalah varian dari SKQD yang menggantikan kebutuhan untuk memilih ansatz tempat men-sampling bitstring dengan sebuah ensemble circuit evolusi-waktu yang dibangun langsung dari Hamiltonian target. Ini dicapai dengan men-subsample operator evolusi-waktu yang lebih kecil dari Hamiltonian berdasarkan koefisiennya, yang dikenal sebagai metode Trotterization qDRIFT.

Tutorial ini menggunakan Qiskit Fermions untuk membuat circuit fermionik yang lebih alami untuk algoritma qDRIFT, diikuti dengan penggunaan pass layout dan sintesis fermionik sebelum memasukkan circuit ke dalam pipeline Qiskit tradisional untuk eksekusi hardware.

Misalkan Hamiltonian berbentuk:

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

dengan, tanpa mengurangi keumuman, kita mensyaratkan $c_i > 0$ dan eigenvalue terbesar dari $h_i$ sama dengan $1$ dalam nilai absolut. Prefaktor 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 dari setiap suku. Di sini $N$ adalah jumlah suku (atau, setelah pengelompokan, jumlah grup) dalam Hamiltonian; ini adalah properti dari Hamiltonian dan berbeda dari jumlah operator yang di-sampling ke dalam satu circuit, yang ditulis $n$ di bawah ini.

Algoritma qDRIFT kemudian merealisasikan, untuk waktu target $t$, sebuah operator $V_k$, di mana $k$ berjalan dari $1 \cdots K$ dan menandakan circuit SqDRIFT ke-$k$, didefinisikan sebagai:

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

Di sini $n$ adalah jumlah operator yang di-sampling per circuit dan $K$ adalah jumlah circuit dalam ensemble. Produk tersebut berjalan pada $n$ pengambilan, bukan pada seluruh $N$ suku Hamiltonian, dan karena suku-suku diambil dengan pengembalian (with replacement), $h_i$ yang sama dapat muncul lebih dari sekali dalam satu $V_k$.

Besaran:

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

adalah norma $L_1$ dari koefisien-koefisien tersebut, sehingga setiap $n$ langkah berevolusi selama durasi yang sama $\lambda t / n$ terlepas dari suku mana yang diambil. Keseragaman sudut langkah adalah fitur khas dari qDRIFT: sebuah koefisien memengaruhi hasil melalui seberapa sering sukunya diambil, bukan melalui seberapa jauh suku tersebut dirotasi. Indeks-indeks tersebut di-sampling dari distribusi:

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

sehingga deret $(k_1, \ldots, k_n)$ adalah barisan acak indeks suku yang diambil dari distribusi ini. Karena $c_i$ positif dan berjumlah $\lambda$, ini adalah distribusi probabilitas yang ternormalisasi, dan ekspektasi channel yang dihasilkan pada pengambilan acak ini mengaproksimasi evolusi di bawah $H$, dengan error yang menurun seiring bertambahnya $n$. Perhatikan bahwa error aproksimasi bergantung pada $\lambda$, bukan pada jumlah suku $N$.

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

Tutorial ini menunjukkan cara menghasilkan ensemble dari circuit yang diacak semacam itu. Setelah kita membuat circuit-circuit ini, mirip dengan cara kita membuat subspace Krylov untuk operator yang berbeda-beda, kita men-sampling bitstring dari beberapa operator semacam itu dengan parameter waktu yang berbeda-beda. Ini memastikan overlap yang lebih tinggi antara vektor ground state dan bitstring yang di-sampling.

Persyaratan

Sebelum memulai tutorial ini, pastikan Anda sudah menginstal

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

Anda dapat menginstal semua paket yang diperlukan dengan:

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

Penyiapan

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

# 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

Langkah 1: Memetakan input klasik ke masalah kuantum

Membaca dan Menyiapkan FCIDump

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

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

File yang setara dapat 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 integral bergantung pada orbital SCF yang telah konvergen, file yang dibuat ulang mungkin berbeda dari file bawaan dalam hal fase atau urutan orbital; total energinya tidak terpengaruh.

Mendapatkan file. Temukan FCIDump di repositori GitHub ini. Anda dapat menjalankan sel di bawah untuk mengambilnya ke lokasi yang diharapkan oleh sisa tutorial ini.

Pertama kita menggunakan cisolver yang disediakan oleh pyscf untuk mendapatkan energi referensi. Ini adalah energi keadaan dasar sebenarnya dari molekul yang kita kerjakan. Untuk ini kita akan mendeklarasikan norb dan nelec terlebih dahulu, yang masing-masing adalah jumlah orbital dan jumlah elektron. Kemudian kita mendeklarasikan h1e dan h2e, yang masing-masing adalah integral satu- dan dua-elektron. Semua ini nantinya juga akan digunakan 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 = "assets/sqdrift/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 = "assets/sqdrift/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 diperlukan sudah siap, kita membaca 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

Alur kerja fermionik dengan qiskit-fermions

Pertama kita akan memetakan Hamiltonian ke dalam model sirkuit fermionik menggunakan qiskit-fermions, yang menyediakan pass transpiler dan gate khusus untuk sirkuit fermionik. Pass ini nantinya akan digunakan sebelum pass transpiler Qiskit yang biasa untuk alur kerja ini.

Pengelompokan suku

Untuk memastikan reproduktibilitas hasil, pertama kita menggunakan canonical_order untuk mengurutkan suku hanya berdasarkan strukturnya. Urutan operator dalam daftar canon dengan demikian bersifat tetap. Ini memastikan reproduktibilitas operator yang dibuat karena pass QDriftTrotterization yang akan kita gunakan nanti mengambil sampel indeks acak untuk membuat operator qDRIFT.

Pada langkah ini, kita memanfaatkan banyak simetri yang ada dalam Hamiltonian struktur elektronik dengan mengelompokkan suku-suku terkait yang memiliki koefisien identik. Meskipun melakukan hal ini mengubah distribusi koefisien operator yang menjadi sumber sampling protokol qDRIFT, hal ini tidak memengaruhi jaminan konvergensinya. Yang penting, pengelompokan suku yang terkait melalui simetri menghasilkan pembatalan suku Pauli yang menguntungkan dan kedalaman sirkuit yang lebih pendek secara keseluruhan saat mengevolusikan 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 dengan urutan normal.

Penyaringan suku diagonal

Kita menghapus suku diagonal dari Hamiltonian yang digunakan untuk membuat sirkuit, sehingga $n$ slot sampling qDRIFT digunakan untuk suku yang memindahkan populasi antar konfigurasi. Suku semacam ini paling baik disaring dari Hamiltonian pada titik ini, sebelum gate Evolution dibuat pada langkah berikutnya.

Suku yang dimaksud adalah suku yang diagonal dalam basis bilangan okupasi, yaitu hasil kali operator bilangan $a^\dagger_i a_i$. Ada tiga jenis suku yang termasuk dalam deskripsi ini:

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

• operator bilangan individual $n_i$, yang evolusi waktunya tereduksi menjadi rotasi $Z$ satu-qubit;

• hasil kali orde lebih tinggi seperti $n_i n_j$.

Dengan sendirinya, tidak satu pun dari suku ini memindahkan populasi antar konfigurasi bilangan okupasi; mereka hanya bekerja pada fase konfigurasi yang sudah ada. Namun, mereka bukan tidak berpengaruh: fase relatif tersebut memberi masukan pada interferensi yang dihasilkan oleh suku eksitasi belakangan dalam sirkuit, sehingga menyaringnya mengubah evolusi yang sebenarnya dihasilkan dan dapat mengubah distribusi sampling. Ini adalah aproksimasi yang disengaja dalam langkah pembuatan sirkuit, dibuat untuk memfokuskan sampling pada suku eksitasi, bukan langkah yang membiarkan distribusi yang disampel tidak berubah. Berbeda dengan pengelompokan simetri di atas, yang membiarkan jaminan konvergensi qDRIFT tetap utuh, filter ini mengubah operator yang dievolusikan. Sirkuit oleh karena itu tidak lagi mengaproksimasi evolusi di bawah Hamiltonian penuh, dan batas kesalahan qDRIFT berlaku untuk operator yang telah difilter, bukan operator aslinya. Hal ini dapat diterima di sini karena sirkuit hanyalah heuristik sampling yang digunakan untuk mengusulkan konfigurasi: tidak ada suku yang hilang dari estimasi energi itu sendiri, karena filter hanya berlaku pada Hamiltonian yang digunakan untuk membangun sirkuit, sementara diagonalisasi klasik nantinya menggunakan Hamiltonian penuh, termasuk suku diagonal. Akurasi SQD bergantung pada langkah klasik tersebut, yang tetap variasional dalam subruang yang disampel terlepas dari bagaimana konfigurasi diusulkan.

Fungsi filter_diagonal_terms() menghapus suku semacam itu dari sebuah operator secara in-place. Fungsi ini mengidentifikasinya dari struktur urutan normalnya — multiset mode kreasi yang cocok dengan multiset mode anihilasi — sehingga hanya valid pada operator yang sudah berurutan normal. Asumsi ini tidak diperiksa 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

Sekarang setelah kita mengelompokkan suku-suku dalam Hamiltonian, kita akan menentukan parameter berikut untuk membuat ansambel sirkuit:

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

Membuat sirkuit fermionik

Sekarang kita akan membuat sirkuit fermionik untuk setiap langkah waktu. Setiap sirkuit akan terdiri dari satu gate evolusi, dengan waktu evolusi yang telah kita deklarasikan sebelumnya. Operator evolusinya adalah Hamiltonian. Selanjutnya kita menjalankan pass transpiler pada sirkuit ini untuk membuat sirkuit qDRIFT.

Persiapan ansatz

Kita menyiapkan keadaan Hartree-Fock menggunakan kelas InitializeModes. Untuk nitrogen, prosesnya cukup dengan menerapkan gate X pada num_elec_a qubit pertama lalu pada num_elec_b qubit, yang keduanya sama dengan tujuh untuk nitrogen. Keadaan ini merepresentasikan tujuh elektron $\alpha$ dan tujuh elektron $\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)

Langkah 2: Optimalkan masalah untuk eksekusi perangkat keras kuantum

Sekarang setelah kita memiliki sirkuit kita, pertama kita akan menggunakan pass yang tersedia di qiskit-fermions untuk melakukan optimalisasi tingkat fermionik, diikuti dengan mentranspile sirkuit kita untuk backend pilihan. Karena ini adalah eksperimen simulator, kita akan melakukan ini terlebih dahulu untuk AerSimulator. Perhitungan bobot untuk setiap kelompok

Pada langkah ini, kita melakukan sampling qDRIFT dari suku-suku secara stokastik dengan probabilitas sebanding dengan koefisiennya dalam Hamiltonian. Pass transpiler qDRIFT melakukan ini untuk kita. Sekarang kita dapat membuat sirkuit yang lebih dangkal yang dapat dieksekusi di perangkat keras dengan lebih efisien meskipun konektivitas qubit terbatas, bahkan ketika Hamiltonian mengandung kopling jarak jauh dan suku yang lebih dari kuadratik. Setelah pengelompokan suku, pass ini menyampel operator berdasarkan bobotnya. Untuk setiap operator $h_i$, bobot $W_{h_i}$ didefinisikan sebagai berikut:

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

Optimalisasi fermionik dan native perangkat keras

Fungsi generate_preset_jw_pass_manager() mengembalikan MultiStagePassManager yang mengambil FermionicCircuit dan menghasilkan sirkuit akhir yang teroptimasi yang dapat kita transpile untuk dijalankan di perangkat keras kita. Kita mengganti tahap optimalisasi bawaannya dengan FermionicPassManager yang berisi pass QDriftTrotterization kita:

• Pass QDriftTrotterization menggunakan perhitungan bobot dan sampling secara internal untuk membuat sirkuit yang akan kita gunakan untuk sampling

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

Tahap-tahap lain dari MultiStagePassManager berjalan secara otomatis dan menangani pemetaan fermion-ke-qubit secara penuh:

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

• F2QSynth: Pass transpilasi untuk memetakan instruksi sirkuit 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 setelah kita selesai dengan optimalisasi tingkat fermionik, kita dapat mentranspile sirkuit untuk eksekusi di simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Langkah 3: Eksekusi menggunakan primitif Qiskit

Sekarang setelah kita memiliki sirkuit kita, kita dapat menjalankannya menggunakan primitif Qiskit di AerSimulator. Kita akan menggabungkan semua count dari berbagai sirkuit. Kita mengubahnya menjadi vektor boolean sebelum akhirnya melakukan post-processing 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: Post-process dan kembalikan hasil dalam format klasik yang diinginkan

Menggunakan bitstring untuk SQD

Sekarang kita dapat menjalankan skema diagonalisasi pada bitstring yang dipilih untuk menemukan eigenvalue terendah yang akan sesuai dengan energi keadaan dasar molekul. Kita membuat fungsi callback, mendeklarasikan okupansi awal, dan mengatur parameter sebelum akhirnya menjalankan skema diagonalisasi. Fungsi callback digunakan untuk mencetak iterasi saat ini dan estimasi eigenvalue saat ini pada setiap iterasi.

Akhirnya, untuk mendapatkan estimasi keadaan dasar, kita menambahkan nuclear_repulsion_energy ke energi hasilnya.

Catatan: Dimensi subruang tidak tetap di seluruh iterasi, bahkan pada simulator tanpa noise — setiap subsampel mengambil sekumpulan konfigurasi yang berbeda, dan langkah recovery membentuk ulang pool di antara iterasi, sehingga dimensi yang dilaporkan bervariasi dari satu subsampel ke subsampel berikutnya. Sampling tanpa noise tidak dengan sendirinya mematok dimensi subruang yang dipilih. Namun, eksekusi di perangkat keras cenderung memberikan subruang yang secara sistematis lebih besar, karena shot yang mengandung noise merusak simetri jumlah partikel dan recovery konfigurasi mengubahnya menjadi vektor basis tambahan. Karena itu, kita juga akan memperkenalkan langkah tambahan untuk pruning bitstring pada bagian perangkat keras.

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 perangkat keras

Contoh ini menggunakan 20 qubit (10 orbital spasial). Pilihan itu adalah kenyamanan untuk tutorial yang harus berjalan cepat, bukan batas atas yang mutlak untuk metode ini.

Biaya langkah klasik tidak ditentukan langsung oleh jumlah qubit. SQD mendiagonalisasi Hamiltonian yang diproyeksikan ke subruang yang dibentang oleh konfigurasi yang disampel, sehingga yang menentukan biaya klasik adalah dimensi subruang terpilih tersebut — di sini ditentukan oleh samples_per_batch, num_batches, dan berapa banyak konfigurasi berbeda yang benar-benar dihasilkan oleh sirkuit — bersama dengan aljabar linear jarang yang diperlukan untuk menerapkan Hamiltonian terproyeksi. Ruang CI penuh tumbuh secara kombinatorial dengan orbital dan elektron, tetapi subruang terpilih adalah irisan kecil dan dapat disesuaikan darinya, dan kita mengendalikan ukurannya secara langsung. Akibatnya, jumlah qubit dan kesulitan klasik dapat bervariasi secara agak independen: ruang orbital yang lebih luas yang disampel ke dalam subruang sedang bisa lebih murah daripada sistem yang lebih kecil yang didiagonalisasi pada subruang yang sangat besar.

Dalam praktiknya, ukuran sistem yang layak bergantung pada dimensi subruang yang Anda butuhkan untuk akurasi yang diinginkan dan pada memori serta core yang tersedia untuk eigensolver. Ruang orbital yang lebih besar biasanya memang memerlukan subruang yang lebih besar untuk mencapai akurasi kimia, dan itulah yang pada akhirnya mendorong penggunaan sumber daya terdistribusi — lihat qiskit-addon-sqd-hpc untuk menskalakan langkah ini. Alih-alih mengasumsikan batas tetap, pendekatan praktisnya adalah mengamati dimensi subruang yang dilaporkan dan konvergensi energi di seluruh iterasi lalu meningkatkan ukuran subruang hingga energi berhenti membaik atau Anda kehabisan memori yang tersedia.

Catatan: Karena kesalahan sampling akibat noise pada perangkat keras, subruang yang dibuat untuk diagonalisasi pada eksekusi perangkat keras akan lebih besar daripada yang kita dapatkan saat menggunakan simulator. Meskipun ini meningkatkan dimensi subruang yang ingin kita diagonalisasi, alur kerja tetap memberi kita jawaban yang akurat berkat ketahanan SQD terhadap noise.

Pruning string yang keliru

Di sini kita dapat memilih untuk melakukan langkah tambahan. Ketika kita memiliki semua bitstring dari eksekusi sirkuit, kita bisa menyaring bitstring yang tidak valid sebelum menjalankan SQD, atau melanjutkan tanpa pruning. Melewati pruning umumnya lebih disukai untuk eksekusi perangkat keras, karena hal ini membiarkan shot yang simetrinya rusak tetap tersedia untuk recovery konfigurasi, yang dapat memperbaikinya menjadi konfigurasi valid dan dengan demikian memperluas subruang alih-alih membuang shot tersebut begitu saja.

Karena nitrogen hanya dapat memiliki tujuh elektron $\alpha$ dan tujuh elektron $\beta$, bitstring apa pun yang memiliki lebih atau kurang dari tujuh angka 1 pada paruh pertama dan paruh kedua output dapat dibuang. Kita mendefinisikan fungsi yang memeriksa apakah bitstring valid, dan jika tidak, membuangnya. Setelah kita menyaring bitstring yang keliru, sisanya dikirim ke skema diagonalisasi. Gunakan flag PRUNE di bawah untuk beralih antara dua perilaku ini.

Perlu diingat bahwa pruning hanyalah salah satu dari beberapa pilihan yang membentuk subruang akhir, selain jumlah sirkuit, kumpulan waktu evolusi, dan penyaringan suku diagonal. Membandingkan eksekusi yang di-pruning dengan yang tidak di-pruning hanya informatif jika semua hal lain tetap sama; pendamping C++ membahas ini lebih rinci, karena ia melakukan postseleksi alih-alih recovery dan juga berbeda dalam parameter lainnya.

name = "assets/sqdrift/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 selanjutnya

Rekomendasi

Jika Anda merasa karya ini menarik, Anda mungkin tertarik dengan materi berikut:

• Diagonalisasi kuantum Krylov berbasis sampel dari model kisi fermionik - tutorial terkait yang menggunakan sirkuit evolusi waktu alih-alih ansatz variasional.
• Diagonalisasi kuantum berbasis sampel dari Hamiltonian kimia - tutorial tentang cara membangun sirkuit local unitary cluster Jastrow (LUCJ) untuk simulasi kimia kuantum.
• Paper SqDRIFT - literatur yang menjadi dasar tutorial ini. (Perhatikan bahwa beberapa optimalisasi yang dibahas dalam paper ini saat ini masih dalam pengembangan, dan tutorial ini dapat berubah di masa depan berdasarkan evolusi pustaka yang digunakan.)