Penyepelasan 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

Selepas mengikuti tutorial ini, pengguna harus memahami:

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

Kami mengesyorkan pengguna agar membiasakan diri dengan topik-topik berikut sebelum mengikuti tutorial ini:

Kimia kuantum dan kuantisasi kedua
Menggunakan primitif Sampler untuk mengambil sampel daripada Circuit kuantum

Latar belakang

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

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 penyepelasan 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:

Pilih ansatz Circuit dan gunakannya pada komputer kuantum kepada keadaan rujukan (dalam kes ini, keadaan Hartree-Fock).
Ambil sampel rentetan bit daripada keadaan kuantum yang terhasil.
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 enam berat, di mana kami boleh menggunakan corak "zig-zag", seperti 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).

Rajah pemetaan qubit untuk ansatz LUCJ pada kekisi segi enam berat

Pemulihan konfigurasi kendiri-konsisten

Prosedur pemulihan konfigurasi kendiri-konsisten direka untuk mengekstrak sebanyak mungkin isyarat daripada sampel kuantum yang bermasalah. 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 bunyi. 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 bunyi 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:

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.
Kumpulkan semua rentetan bit lama dan baharu yang memenuhi simetri, dan subsampel subset bersaiz tetap, dipilih terlebih dahulu.
Untuk setiap subset 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.
Kemas kini anggaran purata penghunian orbital dengan anggaran keadaan dasar yang mempunyai tenaga paling rendah.

Rajah aliran kerja SQD

Aliran kerja SQD digambarkan dalam rajah berikut:

Rajah aliran kerja algoritma SQD

Keperluan

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

Qiskit SDK v1.0 atau lebih baru, dengan sokongan visualisasi
Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)
SQD Qiskit addon v0.11 atau lebih baru (pip install qiskit-addon-sqd)
ffsim v0.0.75 atau lebih baru (pip install ffsim)

Persediaan

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

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 asas 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)], ["N", (1.0, 0, 0)]],
    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.464957764796
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.5933309085008  E_corr = -0.1283731437052354

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 enam berat dan set gate yang secara semula jadi terurai oleh 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 "zig-zag" yang diterangkan dalam bahagian latar belakang tentang ansatz LUCJ. 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 bunyi. 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 = 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
coupling_map = CouplingMap.from_heavy_hex(3)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    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,
    # Limit the number of optimization iterations to prevent the code cell
    # from running too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# 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': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, '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 Qiskit Runtime dan melaksanakan Circuit kami.

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

Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout

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: {valid_fraction}")

Fraction of sampled configurations that are valid: 1.0

Semua rentetan bit adalah sah kerana kami mengambil sampel Circuit pada simulator tanpa bunyi. Apabila dijalankan pada QPU bermasalah, 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 bermasalah 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 bermasalah. Tanpa argumen ini, pemulihan konfigurasi akan gagal dan mengeluarkan ralat jika tiada konfigurasi sah yang diberikan.

from functools import partial

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

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# 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=max_cycle)

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\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,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

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

Iteration 1
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Iteration 2
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Final energy: -108.59275573641656
Final energy error: 0.0032316145694579745

Visualisasi keputusan

Plot pertama menunjukkan bahawa 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 adalah sistem yang kecil, dan kerana sampel adalah tanpa bunyi, pemulihan konfigurasi tidak diperlukan. Pada sistem yang lebih besar yang dijalankan pada QPU bermasalah, beberapa iterasi pemulihan konfigurasi mungkin diperlukan, dan ketepatan akhir mungkin lebih rendah. Secara umum, 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]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# 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_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
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()

Output of the previous code cell

Contoh perkakasan skala besar

Kini, kami menjalankan contoh yang lebih besar pada perkakasan kuantum sebenar. Di sini, kami akan memperoleh ruang aktif untuk molekul nitrogen daripada set asas cc-pVDZ.

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)], ["N", (1.0, 0, 0)]],
    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, min_num_qubits=133
)
print(f"Using backend {backend.name}")

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    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,
    # Limit the number of optimization iterations to prevent the code cell
    # from running too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# 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: {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

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# 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=max_cycle)

# 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,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {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]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# 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_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
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_boston

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': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.02124
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.13889134249762
		Subspace dimension: 120409
	Subsample 1
		Energy: -109.11785470455858
		Subspace dimension: 110889
	Subsample 2
		Energy: -109.13234360554011
		Subspace dimension: 130321
Iteration 2
	Subsample 0
		Energy: -109.16392179579177
		Subspace dimension: 223729
	Subsample 1
		Energy: -109.16281938332986
		Subspace dimension: 223729
	Subsample 2
		Energy: -109.16955816711932
		Subspace dimension: 233289
Iteration 3
	Subsample 0
		Energy: -109.17905772999075
		Subspace dimension: 324900
	Subsample 1
		Energy: -109.17532445048462
		Subspace dimension: 357604
	Subsample 2
		Energy: -109.1733168689756
		Subspace dimension: 348100
Iteration 4
	Subsample 0
		Energy: -109.18437778820451
		Subspace dimension: 474721
	Subsample 1
		Energy: -109.18450164209159
		Subspace dimension: 476100
	Subsample 2
		Energy: -109.18493571190754
		Subspace dimension: 487204
Iteration 5
	Subsample 0
		Energy: -109.18616522497996
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.18652868888333
		Subspace dimension: 644809
	Subsample 2
		Energy: -109.18753326484406
		Subspace dimension: 585225
Final energy: -109.18753326484406
Final energy error: 0.040495951813099396

Output of the previous code cell

Langkah seterusnya

Cadangan

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

Penyepelasan 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 penyepelasan
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

Hasil pembelajaran​

Prasyarat​

Latar belakang​

Kimia kuantum​

Ansatz kluster unitari tempatan Jastrow (LUCJ)​

Pemulihan konfigurasi kendiri-konsisten​

Rajah aliran kerja SQD​

Keperluan​

Persediaan​

Contoh simulator skala kecil​

Langkah 1: Petakan input klasik ke masalah kuantum​

Langkah 2: Optimumkan untuk pelaksanaan perkakasan kuantum​

Langkah 3: Jalankan menggunakan primitif Qiskit​

Langkah 4: Pasca-proses dan kembalikan keputusan dalam format klasik yang dikehendaki​

Visualisasi keputusan​

Contoh perkakasan skala besar​

Langkah 1-4​

Langkah seterusnya​