Pepenjuruan kuantum berasaskan sampel bagi Hamiltonian kimia

Anggaran penggunaan: kurang dari satu minit pada pemproses Heron r2 (NOTA: Ini adalah anggaran sahaja. Masa jalan anda mungkin berbeza.)

Hasil pembelajaran

• Cara menggunakan SQD Qiskit addon untuk menghampiri tenaga keadaan dasar sistem molekul menggunakan rentetan bit yang disampling daripada unit pemprosesan kuantum (QPU).

• Cara menggunakan ffsim untuk membina Circuit kluster unitari tempatan Jastrow (LUCJ) bagi simulasi kimia kuantum.

Prasyarat

• Kimia kuantum dan pengkuantuman kedua

• Menggunakan primitif Sampler untuk mengambil sampel daripada Circuit kuantum

Latar belakang

Dalam tutorial ini, kami menunjukkan cara memproses sampel kuantum yang berhingar untuk menghampiri keadaan dasar molekul nitrogen $\text{N}_2$ pada panjang ikatan keseimbangan, dengan menggunakan SQD Qiskit addon untuk melaksanakan algoritma pepenjuruan kuantum berasaskan sampel (SQD). Maklumat lanjut tentang perisian boleh didapati dalam dokumentasi yang berkaitan, termasuk contoh mudah untuk memulakan.

Kimia nitrogen menghubungkan contoh ini dengan kemampanan melalui pengeluaran ammonia untuk baja. Memahami cara pemangkin mengaktifkan nitrogen adalah sebahagian daripada penyelidikan ke arah kaedah pengeluaran yang kurang intensif tenaga. Mengira keadaan dasar bagi molekul nitrogen terasing memperkenalkan kaedah struktur elektronik yang relevan dengan cabaran yang lebih luas itu.

Tutorial ini disyorkan untuk pengguna yang biasa dengan kimia kuantum: khususnya, kebiasaan dengan mencari tenaga keadaan dasar molekul. Untuk panduan terperinci tentang aliran kerja, rujuk kursus algoritma pepenjuruan kuantum.

SQD ialah teknik untuk mencari nilai eigen dan vektor eigen bagi operator kuantum, seperti Hamiltonian sistem kuantum, dengan menggunakan pengkomputeran kuantum dan pengkomputeran klasik teragih bersama-sama. Pengkomputeran klasik teragih digunakan untuk memproses sampel yang diperoleh daripada pemproses kuantum, dan untuk mengunjur dan menyepelasi Hamiltonian sasaran dalam subruang yang dibentangnya. Aliran kerja berasaskan SQD mempunyai langkah-langkah berikut:

1. Pilih ansatz Circuit dan gunakannya pada komputer kuantum kepada keadaan rujukan (dalam kes ini, keadaan Hartree-Fock).

2. Ambil sampel rentetan bit daripada keadaan kuantum yang terhasil.

3. Jalankan prosedur pemulihan konfigurasi kendiri-konsisten pada rentetan bit untuk mendapatkan penghampiran kepada keadaan dasar.

SQD diketahui berfungsi dengan baik apabila eigenstate sasaran adalah jarang: fungsi gelombang disokong dalam set keadaan asas $\mathcal{S} = \{|x\rangle \}$ yang saiznya tidak meningkat secara eksponen dengan saiz masalah.

Kimia kuantum

Hamiltonian sistem molekul boleh 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}$ ialah nombor kompleks yang dipanggil integral molekul yang boleh dikira daripada spesifikasi molekul menggunakan program komputer. Dalam tutorial ini, kami mengira integral menggunakan pakej perisian PySCF.

Untuk butiran tentang bagaimana Hamiltonian molekul diterbitkan, rujuk buku teks kimia kuantum (contohnya, Modern Quantum Chemistry oleh Szabo dan Ostlund). Untuk penjelasan peringkat tinggi tentang cara masalah kimia kuantum dipetakan ke komputer kuantum, semak kuliah Mapping Problems to qubits dari Qiskit Global Summer School 2024.

Ansatz kluster unitari tempatan Jastrow (LUCJ)

SQD memerlukan ansatz Circuit kuantum untuk mengambil sampel. Dalam tutorial ini, kami akan menggunakan ansatz kluster unitari tempatan Jastrow (LUCJ) kerana gabungan motivasi fizikal dan keramahan perkakasannya. Kami akan menggunakan ffsim untuk membina Circuit ansatz.

Ansatz LUCJ menyesuaikan diri dengan QPU yang mempunyai kesambungan qubit terhad. Orbital spin dipetakan ke qubit supaya ansatz tidak memerlukan penghalaan dengan gate SWAP. Perkakasan IBM® mempunyai topologi qubit kekisi segi empat atau kekisi segi enam berat. Untuk kekisi segi enam berat, kami boleh menggunakan corak "zig-zag" yang digambarkan di bawah. Dalam corak ini, orbital berspin sama dipetakan ke qubit dengan topologi garisan (bulatan merah dan biru), dan sambungan antara orbital berspin berbeza wujud pada setiap orbital ruang ke-4, dengan sambungan difasilitasi oleh qubit ancilla (bulatan ungu).

Pemulihan konfigurasi kendiri-konsisten

Prosedur pemulihan konfigurasi kendiri-konsisten direka untuk mengekstrak sebanyak mungkin isyarat daripada sampel kuantum yang berhingar. Kerana Hamiltonian molekul memelihara nombor zarah dan spin Z, masuk akal untuk memilih ansatz Circuit yang juga memelihara simetri ini. Apabila digunakan pada keadaan Hartree-Fock, keadaan yang terhasil mempunyai nombor zarah dan spin Z yang tetap dalam tetapan tanpa hingar. Oleh itu, bahagian spin-$\alpha$ dan spin-$\beta$ bagi mana-mana rentetan bit yang disampling daripada keadaan ini harus mempunyai berat Hamming yang sama seperti dalam keadaan Hartree-Fock. Disebabkan kehadiran hingar dalam pemproses kuantum semasa, sesetengah rentetan bit yang diukur akan melanggar sifat ini. Bentuk pemilihan pasca yang mudah adalah membuang rentetan bit ini, tetapi ini membazir kerana rentetan bit tersebut mungkin masih mengandungi beberapa isyarat. Prosedur pemulihan kendiri-konsisten cuba memulihkan sebahagian isyarat itu dalam pasca-pemprosesan. Prosedur ini adalah berulang dan memerlukan sebagai input anggaran purata penghunian setiap orbital dalam keadaan dasar, yang pertama kali dikira daripada sampel mentah. Prosedur dijalankan dalam gelung, dan setiap iterasi mempunyai langkah-langkah berikut:

1. Untuk setiap rentetan bit yang melanggar simetri yang ditetapkan, balikkan bitnya dengan prosedur kebarangkalian yang direka untuk membawa rentetan bit lebih dekat kepada anggaran semasa purata penghunian orbital, untuk mendapatkan rentetan bit baharu.

2. Kumpulkan semua rentetan bit lama dan baharu yang memenuhi simetri, dan bahagikan subsampel daripadanya kepada kelompok tidak bertindih bersaiz tetap, yang dipilih terlebih dahulu.

3. Untuk setiap kelompok rentetan bit, unjurkan Hamiltonian ke dalam subruang yang dibentang oleh vektor asas yang sepadan (lihat bahagian sebelumnya untuk penerangan vektor asas ini), dan kira anggaran keadaan dasar bagi Hamiltonian yang diunjurkan pada komputer klasik. Setiap rentetan bit terdiri daripada separuh spin-$\alpha$ dan separuh spin-$\beta$, yang masing-masing kami panggil rentetan configuration interaction (CI). Subruang bagi sesuatu kelompok dibentang oleh vektor asas yang memasangkan setiap rentetan CI spin-$\alpha$ dalam kelompok dengan setiap rentetan CI spin-$\beta$ dalam kelompok.

4. Pangkas setiap kelompok supaya hanya rentetan CI dengan pemberat terbesar dalam anggaran keadaan dasar kelompok itu yang kekal. Gabungkan rentetan CI yang tinggal daripada semua kelompok, dan diagonalkan Hamiltonian sekali lagi dalam subruang yang digabungkan.

5. Kemas kini anggaran purata penghunian orbital dengan anggaran keadaan dasar daripada subruang yang digabungkan. Rentetan CI berpemberat tertinggi bagi subruang yang digabungkan juga dibawa ke dalam subruang lelaran seterusnya.

Operasi pangkas dan gabung dalam langkah 4 membentuk varian trim SQD bagi algoritma ini. Pendiagonalan satu kelompok menelan kos yang lebih rendah daripada pendiagonalan subruang yang memuatkan setiap rentetan CI calon. Oleh itu, dengan pemangkasan, SQD boleh memilih daripada kumpulan rentetan CI bersampel yang lebih besar daripada yang mampu didiagonalkan sekali gus.

Rajah aliran kerja SQD

Rajah berikut menunjukkan aliran kerja SQD. Rajah ini tidak menunjukkan operasi pangkas dan gabung daripada langkah 4.

Keperluan

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

• Qiskit SDK v1.0 atau lebih baru, dengan sokongan visualisasi

• Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)

• SQD Qiskit addon v0.14 atau lebih baru (pip install qiskit-addon-sqd)

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

Persediaan

# 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, kami akan mencari penghampiran kepada keadaan dasar molekul nitrogen berhampiran jarak ikatan keseimbangannya. Kami mula-mula menggunakan set basis STO-6G yang kecil supaya kami dapat mensimulasikan eksperimen dan memastikan ia berfungsi.

Langkah 1: Petakan input klasik ke masalah kuantum

Pertama, kami menentukan molekul dan sifat-sifatnya.

# 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 membina Circuit ansatz LUCJ, kami terlebih dahulu menjalankan pengiraan CCSD dalam sel kod berikut. Amplitud $t_1$ dan $t_2$ daripada pengiraan ini akan digunakan untuk memulakan 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

Kini, kami menggunakan ffsim untuk mencipta Circuit ansatz. Memandangkan molekul kami mempunyai keadaan Hartree-Fock cangkerang tertutup, kami menggunakan varian ansatz UCJ yang seimbang spin, UCJOpSpinBalanced. Kami menetapkan optimize=True dalam kaedah from_t_amplitudes untuk mendayakan pemfaktoran berganda "dimampatkan" bagi amplitud $t_2$ (lihat The local unitary cluster Jastrow (LUCJ) ansatz dalam dokumentasi ffsim untuk butiran).

Oleh kerana ansatz LUCJ menyesuaikan diri dengan kesambungan QPU yang tersedia, kami perlu memulakan backend QPU sebelum mencipta ansatz. Buat masa ini, kami akan mencipta backend generik dengan peta gandingan segi empat dan set gate yang menjadi hasil penguraian semula jadi ansatz LUCJ. Kemudian, kami akan menggunakan ffsim.qiskit.generate_lucj_pass_manager untuk mencipta pengurus laluan yang dikhususkan untuk mentranspilkan ansatz LUCJ ke backend yang diberikan mengikut susun atur "ladder" (lihat penjelasan di dokumentasi ffsim). Fungsi ini menggunakan heuristik pemarkahan untuk meminimumkan ralat yang berkaitan dengan susun atur yang dipilih, yang penting jika backend anda adalah QPU sebenar atau simulator dengan model hingar. Selain mengembalikan pengurus laluan, fungsi ini juga mengembalikan pasangan gandingan alpha-beta yang boleh dilaksanakan pada perkakasan. Jika tidak semua pasangan boleh dilaksanakan, ia akan mengeluarkan amaran.

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: Optimumkan untuk pelaksanaan perkakasan kuantum

Seterusnya, kami mengoptimumkan Circuit untuk perkakasan sasaran. Biasanya, langkah ini melibatkan pemulaan backend perkakasan dan pengurus laluan untuk backend tersebut. Walau bagaimanapun, oleh kerana ansatz LUCJ disesuaikan dengan kesambungan perkakasan, kami telah menjalankan tindakan ini dalam langkah sebelumnya. Yang tinggal hanyalah menjalankan pengurus laluan pada Circuit untuk mentranspilkannya kepada Circuit ISA yang boleh dilaksanakan terus pada 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: Jalankan menggunakan primitif Qiskit

Selepas mengoptimumkan Circuit untuk pelaksanaan perkakasan, kami sedia untuk menjalankannya pada perkakasan sasaran dan mengumpul sampel untuk anggaran tenaga keadaan dasar. Memandangkan kami hanya mempunyai satu Circuit, kami akan menggunakan mod pelaksanaan Job IBM Quantum Compute Service dan melaksanakan Circuit kami.

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 keputusan dalam format klasik yang dikehendaki

Metrik yang berguna untuk menilai kualiti output QPU ialah bilangan konfigurasi sah yang dikembalikan. Konfigurasi sah mempunyai nombor zarah dan spin Z yang betul, bermakna bahagian kanan rentetan bit mempunyai berat Hamming yang sama dengan bilangan elektron spin-atas, dan bahagian kiri mempunyai berat Hamming yang sama dengan bilangan elektron spin-bawah. Sel berikut mengira pecahan konfigurasi yang disampling yang sah.

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 rentetan bit adalah sah kerana kami mengambil sampel Circuit pada simulator tanpa hingar. Apabila dijalankan pada QPU berhingar, pecahannya akan kurang daripada satu, tetapi diharap ia lebih besar daripada pecahan yang dijangka jika rentetan bit disampling secara seragam rawak, yang dikira dalam 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

Kini, kami menganggarkan tenaga keadaan dasar Hamiltonian menggunakan fungsi diagonalize_fermionic_hamiltonian. Fungsi ini menjalankan prosedur pemulihan konfigurasi kendiri-konsisten untuk secara berulang memperhalusi sampel kuantum yang berhingar bagi meningkatkan anggaran tenaga. Kami menghantar fungsi panggil balik supaya kami dapat menyimpan keputusan pertengahan untuk analisis kemudian. Lihat dokumentasi API untuk penjelasan argumen bagi diagonalize_fermionic_hamiltonian.

Di sini, kami menggunakan argumen initial_occupancies kepada diagonalize_fermionic_hamiltonian untuk menentukan konfigurasi Hartree-Fock sebagai tekaan awal bagi penghunian orbital dalam keadaan dasar. Pendekatan ini munasabah untuk sistem di mana keadaan dasar mempunyai sokongan yang ketara pada konfigurasi Hartree-Fock, tetapi mungkin tidak sesuai dalam situasi lain, walaupun kaedah pengiraan lebih canggih mungkin menghasilkan tekaan awal yang lebih baik dalam kes tersebut. Menentukan initial_occupancies juga membolehkan pemulihan konfigurasi dijalankan walaupun tiada konfigurasi sah yang disampling, seperti yang mungkin berlaku apabila mengambil sampel Circuit yang besar pada QPU berhingar. Tanpa argumen ini, pemulihan konfigurasi akan gagal dan mengeluarkan ralat jika tiada konfigurasi sah yang diberikan.

Kami mengaktifkan trim SQD dengan menghantar TrimPolicy sebagai argumen policy. Setiap lelaran kemudian melakukan num_batches + 1 pendiagonalan: satu untuk setiap kelompok, dan satu untuk subruang yang digabungkan. Kami menggunakan dua kelompok, jadi setiap lelaran memerlukan tiga pendiagonalan. Callback hanya menerima hasil yang digabungkan, jadi setiap lelaran mencetak satu tenaga. Pilihan TrimPolicy ialah:

• trim_ratio: Pecahan rentetan CI setiap kelompok, bagi setiap sektor spin, yang bertahan selepas pemangkasan dan memasuki subruang yang digabungkan. Subruang yang digabungkan memuatkan sehingga kira-kira num_batches * trim_ratio kali ganda bilangan rentetan CI bagi satu kelompok, tolak rentetan yang sama dalam kelompok-kelompok tersebut. Argumen max_dim menghadkan setiap subruang, termasuk subruang yang digabungkan.

• carryover_threshold dan carryover_ratio: Pilihan ini menentukan rentetan CI subruang yang digabungkan yang dibawa ke lelaran seterusnya. Untuk dibawa ke hadapan, rentetan CI mesti muncul dalam konfigurasi yang amplitudnya melebihi carryover_threshold. Daripada rentetan tersebut, pecahan carryover_ratio dengan pemberat terbesar dikekalkan.

Sistem ini cukup kecil sehingga satu kelompok mengandungi kebanyakan rentetan CI yang mungkin, jadi kelompok-kelompok itu sangat bertindih. Oleh itu, kami menggunakan trim_ratio yang tinggi supaya subruang yang digabungkan tidak menjadi jauh lebih kecil daripada kelompok. Lihat panduan trim SQD untuk butiran lanjut tentang memilih pilihan 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

Visualisasi keputusan

Untuk setiap iterasi, plot pertama menunjukkan ralat tenaga yang dikira dalam subruang yang digabungkan. Dalam simulasi ini, kami sudah berada dalam jarak 1 mH daripada jawapan tepat selepas iterasi pertama (ketepatan kimia biasanya diterima sebagai 1 kcal/mol $\approx$ 1.6 mH). Ini ialah sistem yang kecil, dan kerana sampel tanpa hingar, pemulihan konfigurasi tidak diperlukan. Pada sistem yang lebih besar yang dijalankan pada QPU berhingar, beberapa iterasi pemulihan konfigurasi mungkin diperlukan, dan ketepatan akhir mungkin lebih teruk. Secara amnya, tenaga boleh ditingkatkan dengan membenarkan lebih banyak iterasi pemulihan konfigurasi atau dengan meningkatkan bilangan sampel setiap batch.

Plot kedua menunjukkan purata penghunian setiap orbital ruang selepas iterasi akhir. Kami dapat melihat bahawa kedua-dua elektron spin-atas dan spin-bawah menghuni lima orbital pertama dengan kebarangkalian yang tinggi dalam penyelesaian kami.

# 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 perkakasan skala besar

Sekarang, kami menjalankan contoh yang lebih besar pada perkakasan kuantum sebenar. Di sini, kami akan mendapatkan ruang aktif untuk molekul nitrogen daripada set basis cc-pVDZ. Kami juga akan membuat perubahan berikut:

• Rendahkan n_reps kepada 1 untuk mengurangkan hingar daripada ralat gate.
• Tetapkan kesambungan ansatz berdasarkan keluarga pemproses backend (segi enam berat untuk Heron, segi empat untuk Nighthawk).
• Naikkan max_dim kepada 1000 untuk membenarkan dimensi subruang yang lebih besar.
• Rendahkan trim_ratio kepada 0.5. Sistem ini mempunyai jauh lebih banyak rentetan CI yang mungkin daripada yang boleh ditampung oleh subruang. Dengan dua kelompok, nisbah 0.5 memastikan subruang yang digabungkan tidak lebih besar daripada satu kelompok.

Langkah 1-4

Di sini kami menyatukan semua langkah ke dalam satu aliran kerja pada skala yang lebih besar, yang kemudian dijalankan pada perkakasan kuantum sebenar.

# ------------------------------ 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 seterusnya

Cadangan

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

• Pepenjuruan kuantum Krylov berasaskan sampel bagi model kekisi fermionic - tutorial berkaitan menggunakan Circuit evolusi masa dan bukannya ansatz variasi
• Skala aliran kerja kimia SQD dengan penyelesai Dice - halaman yang menunjukkan cara menggunakan perisian Dice yang lebih cekap untuk pepenjuruan
• Dokumentasi API addon SQD - rujukan untuk fungsi diagonalize_fermionic_hamiltonian
• Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer - kertas kerja yang menjadi asas tutorial ini