Meningkatkan anggaran tenaga Hamiltonian kimia dengan SQD

Dalam tutorial ini, kita melaksanakan corak Qiskit yang menunjukkan cara pasca-pemprosesan sampel kuantum berbising untuk mencari anggaran kepada keadaan dasar Hamiltonian kimia: molekul $N_2$ pada keseimbangan dalam set asas 6-31G. Kita akan mengikuti pendekatan pepenjuruan (diagonalisasi) kuantum berasaskan sampel untuk memproses sampel yang diambil daripada ansatz litar kuantum 36-qubit (dalam kes ini, litar LUCJ). Untuk mengambil kira kesan hingar kuantum, teknik pemulihan konfigurasi digunakan.

Corak ini boleh digambarkan dalam empat langkah:

Langkah 1: Petakan kepada masalah kuantum
- Jana ansatz untuk menganggar keadaan dasar
Langkah 2: Optimumkan masalah
- Transpile ansatz untuk backend
Langkah 3: Jalankan eksperimen
- Ambil sampel daripada ansatz menggunakan primitif Sampler
Langkah 4: Pasca-proses keputusan
- Gelung pemulihan konfigurasi swakonsisten
  - Pasca-proses keseluruhan set sampel bitstring, menggunakan pengetahuan awal tentang bilangan zarah dan purata penghunian orbital yang dikira pada lelaran terkini.
  - Cipta kelompok subsampel secara kebarangkalian daripada bitstring yang telah dipulihkan.
  - Unjurkan dan pepenjurukan Hamiltonian molekul bagi setiap subruang yang disampel.
  - Simpan tenaga keadaan dasar minimum yang dijumpai merentasi semua kelompok dan kemaskini purata penghunian orbital.

Untuk contoh ini, Hamiltonian elektron berinteraksi mengambil bentuk umum:

\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ adalah operator penciptaan/pemusnahan fermionic yang berkaitan dengan elemen set asas ke- $p$ dan spin $\sigma$ . $h_{pr}$ dan $(pr|qs)$ adalah kamiran elektronik satu badan dan dua badan.

Aliran kerja SQD dengan pemulihan konfigurasi swakonsisten digambarkan dalam rajah berikut.

Rajah SQD

SQD diketahui berfungsi baik apabila keadaan eigen sasaran adalah jarang: fungsi gelombang disokong dalam set keadaan asas $\mathcal{S} = \{|x\rangle \}$ yang saiznya tidak meningkat secara eksponen dengan saiz masalah. Dalam senario ini, pepenjuruan Hamiltonian yang diunjurkan ke dalam subruang yang ditentukan oleh $\mathcal{S}$ :

H_\mathcal{S} = P_\mathcal{S} H P_\mathcal{S} \textrm{ with } P_\mathcal{S} = \sum_{x \in \mathcal{S}} |x \rangle \langle x |;

menghasilkan anggaran yang baik kepada keadaan eigen sasaran. Peranan peranti kuantum adalah untuk menghasilkan sampel ahli $\mathcal{S}$ sahaja. Pertama, litar kuantum menyediakan keadaan $|\Psi\rangle$ dalam peranti kuantum. Pengekodan Jordan-Wigner digunakan. Oleh itu, ahli asas pengkomputeran mewakili keadaan Fock (konfigurasi/penentu elektronik). Litar disampel dalam asas pengkomputeran, menghasilkan set konfigurasi berbising $\tilde{\mathcal{X}}$ . Konfigurasi diwakili oleh bitstring. Set $\tilde{\mathcal{X}}$ kemudian dihantar ke blok pasca-pemprosesan klasik, di mana teknik pemulihan konfigurasi swakonsisten digunakan. Dalam rangka kerja SQD, peranan peranti kuantum adalah untuk menghasilkan taburan kebarangkalian.

Langkah 1: Petakan masalah kepada litar kuantum

Dalam tutorial ini, kita akan menganggar tenaga keadaan dasar molekul $N_2$ . Pertama, kita akan menentukan molekul dan sifat-sifatnya. Seterusnya, kita akan mencipta ansatz kluster Jastrow unitari setempat (LUCJ) (litar kuantum) untuk menjana sampel daripada komputer kuantum bagi anggaran tenaga keadaan dasar.

Pertama, kita akan menentukan molekul dan sifat-sifatnya.

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

warnings.filterwarnings("ignore")

import pyscf
import pyscf.cc
import pyscf.mcscf

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    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()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
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), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570775
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

Seterusnya, kita akan mencipta ansatz. Ansatz LUCJ adalah litar kuantum berparameter, dan kita akan memulakannya dengan amplitud t2 dan t1 yang diperoleh daripada pengiraan CCSD.

# 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) = -109.0398256929734  E_corr = -0.2045891221988311

Kita akan menggunakan pakej ffsim untuk mencipta dan memulakan ansatz dengan amplitud t2 dan t1 yang dikira di atas. Memandangkan molekul kita mempunyai keadaan Hartree-Fock cangkerang tertutup, kita akan menggunakan varian seimbang spin ansatz UCJ, UCJOpSpinBalanced.

Memandangkan perkakasan IBM sasaran kita mempunyai topologi heksagon berat, kita akan menggunakan corak zig-zag untuk interaksi qubit. Dalam corak ini, orbital (diwakili oleh qubit) dengan spin yang sama disambungkan dengan topologi garis (bulatan merah dan biru) di mana setiap garis mengambil bentuk zig-zag disebabkan oleh sambungan heksagon berat perkakasan sasaran. Sekali lagi, disebabkan oleh topologi heksagon berat, orbital untuk spin berbeza mempunyai sambungan antara setiap orbital ke-4 (0, 4, 8, dll.) (bulatan ungu).

lucj_ansatz

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 1
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, 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(num_orbitals, nelec), qubits)

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

Langkah 2: Optimumkan masalah

Seterusnya, kita akan mengoptimumkan litar kita untuk perkakasan sasaran. Kita perlu memilih peranti perkakasan yang hendak digunakan sebelum mengoptimumkan litar kita. Kita akan menggunakan backend palsu 127-qubit daripada qiskit_ibm_runtime untuk mensimulasikan peranti sebenar. Untuk menjalankan pada QPU sebenar, gantikan sahaja backend palsu dengan backend sebenar. Semak dokumen Qiskit IBM Runtime untuk maklumat lanjut.

from qiskit_ibm_runtime.fake_provider import FakeSherbrooke

backend = FakeSherbrooke()

Seterusnya, kami mengesyorkan langkah-langkah berikut untuk mengoptimumkan ansatz dan menjadikannya serasi dengan perkakasan.

Pilih qubit fizikal (initial_layout) daripada perkakasan sasaran yang mematuhi corak zig-zag yang diterangkan di atas. Menyusun qubit dalam corak ini menghasilkan litar yang cekap dan serasi dengan perkakasan dengan gate yang lebih sedikit.
Jana pengurus laluan berperingkat menggunakan fungsi generate_preset_pass_manager daripada Qiskit dengan pilihan backend dan initial_layout anda.
Tetapkan peringkat pre_init pengurus laluan berperingkat anda kepada ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT termasuk laluan transpiler Qiskit yang menguraikan gate kepada putaran orbital dan kemudian menggabungkan putaran orbital, menghasilkan gate yang lebih sedikit dalam litar akhir.
Jalankan pengurus laluan pada litar anda.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]
initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Gate counts (w/o pre-init passes): OrderedDict({'rz': 4420, 'sx': 3432, 'ecr': 1366, 'x': 239, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 2460, 'sx': 2156, 'ecr': 730, 'x': 71, 'measure': 32, 'barrier': 1})

Langkah 3: Jalankan eksperimen

Selepas mengoptimumkan litar untuk pelaksanaan perkakasan, kita bersedia untuk menjalankannya pada perkakasan sasaran dan mengumpul sampel untuk anggaran tenaga keadaan dasar. Memandangkan kita hanya mempunyai satu litar, kita akan menggunakan mod pelaksanaan job Qiskit Runtime dan melaksanakan litar kita.

Nota: Kami telah mengkomen kod untuk menjalankan litar pada QPU dan meninggalkannya sebagai rujukan pengguna. Daripada menjalankan pada perkakasan sebenar dalam panduan ini, kita hanya akan menjana sampel rawak yang diambil daripada taburan seragam.

import numpy as np
from qiskit_addon_sqd.counts import generate_bit_array_uniform

# from qiskit_ibm_runtime import SamplerV2 as Sampler

# sampler = Sampler(mode=backend)
# job = sampler.run([isa_circuit], shots=10_000)
# primitive_result = job.result()
# pub_result = primitive_result[0]
# bit_array = pub_result.data.meas

rng = np.random.default_rng(24)
bit_array = generate_bit_array_uniform(10_000, num_orbitals * 2, rand_seed=rng)

Langkah 4: Pasca-proses keputusan

Sekarang, kita jalankan algoritma SQD menggunakan fungsi diagonalize_fermionic_hamiltonian. Lihat dokumentasi API untuk penjelasan argumen fungsi ini.

Penyelesai yang disertakan dalam tambahan SQD menggunakan pelaksanaan CI terpilih PySCF, khususnya pyscf.fci.selected_ci.kernel_fixed_space. Contoh di bawah juga menunjukkan cara menghantar argumen kata kunci kepada fungsi tersebut melalui penyelesai yang disertakan. Di sini kita menghantar argumen max_cycle.

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 = 1
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# 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=num_orbitals,
    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,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

Iteration 1
	Subsample 0
		Energy: -105.45358671756313
		Subspace dimension: 5476
Iteration 2
	Subsample 0
		Energy: -107.95172900082163
		Subspace dimension: 249001
Iteration 3
	Subsample 0
		Energy: -108.97460330369815
		Subspace dimension: 339889
Iteration 4
	Subsample 0
		Energy: -109.02739376648793
		Subspace dimension: 440896
Iteration 5
	Subsample 0
		Energy: -109.030972328451
		Subspace dimension: 597529

Sekarang, kita plot keputusan.

Plot pertama menunjukkan bahawa selepas beberapa lelaran kita menganggar tenaga keadaan dasar dalam julat ~16 mH (ketepatan kimia biasanya diterima sebagai 1 kcal/mol $\approx$ 1.6 mH). Ingat, sampel kuantum dalam demo ini ialah hingar tulen. Isyarat di sini datang daripada pengetahuan a priori tentang struktur elektronik dan Hamiltonian molekul.

Plot kedua menunjukkan purata penghunian setiap orbital spatial selepas lelaran terakhir. Kita dapat lihat bahawa kedua-dua elektron spin-naik dan spin-turun menghuni lima orbital pertama dengan kebarangkalian tinggi dalam penyelesaian kita.

import matplotlib.pyplot as plt

# 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 - exact_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})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Exact energy: -109.04667 Ha
SQD energy: -109.03097 Ha
Absolute error: 0.01570 Ha

Plot output

Langkah 1: Petakan masalah kepada litar kuantum​

Langkah 2: Optimumkan masalah​

Langkah 3: Jalankan eksperimen​

Langkah 4: Pasca-proses keputusan​

Langkah 1: Petakan masalah kepada litar kuantum

Langkah 2: Optimumkan masalah

Langkah 3: Jalankan eksperimen

Langkah 4: Pasca-proses keputusan