Meningkatkan anggaran tenaga Hamiltonian kimia dengan SQD

Dalam tutorial ini, kita melaksanakan corak Qiskit yang menunjukkan cara memproses semula 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 penjajaran kuantum berasaskan sampel untuk memproses sampel yang diambil daripada ansatz Circuit kuantum 36-Qubit (dalam kes ini, Circuit LUCJ). Untuk mengambil kira kesan bunyi 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: Proses semula keputusan
- Gelung pemulihan konfigurasi yang konsisten sendiri
  - Proses semula keseluruhan set sampel bitstring, menggunakan pengetahuan awal tentang nombor zarah dan purata penghunian orbit yang dikira pada lelaran terkini.
  - Cipta kelompok subsampel daripada bitstring yang dipulihkan secara kebarangkalian.
  - Unjurkan dan jajarkan Hamiltonian molekul merentasi setiap subruang yang disampling.
  - Simpan tenaga keadaan dasar minimum yang dijumpai merentasi semua kelompok dan kemaskini purata penghunian orbit.

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 yang konsisten sendiri 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, penjajaran 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, Circuit kuantum menyediakan keadaan $|\Psi\rangle$ dalam peranti kuantum. Pengekodan Jordan-Wigner digunakan. Oleh itu, ahli asas pengkomputeran mewakili keadaan Fock (konfigurasi/penentu elektronik). Circuit disampling 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 yang konsisten sendiri digunakan. Dalam rangka kerja SQD, peranan peranti kuantum adalah untuk menghasilkan taburan kebarangkalian.

Langkah 1: Petakan masalah kepada Circuit 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) (Circuit kuantum) untuk menjana sampel daripada komputer kuantum bagi anggaran tenaga keadaan dasar.

Pertama, kita akan menentukan molekul dan sifat-sifatnya.

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

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 Circuit 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, orbit (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, orbit untuk spin berbeza mempunyai sambungan antara setiap orbit 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 Circuit kita untuk perkakasan sasaran. Kita perlu memilih peranti perkakasan yang hendak digunakan sebelum mengoptimumkan Circuit 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 Circuit 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 orbit dan kemudian menggabungkan putaran orbit, menghasilkan Gate yang lebih sedikit dalam Circuit akhir.
Jalankan pengurus laluan pada Circuit 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 Circuit untuk pelaksanaan perkakasan, kita bersedia untuk menjalankannya pada perkakasan sasaran dan mengumpul sampel untuk anggaran tenaga keadaan dasar. Memandangkan kita hanya mempunyai satu Circuit, kita akan menggunakan mod pelaksanaan Kerja Qiskit Runtime dan melaksanakan Circuit kita.

Nota: Kami telah mengkomen kod untuk menjalankan Circuit 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: Proses semula 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 adalah bunyi tulen. Isyarat di sini datang daripada pengetahuan a priori tentang struktur elektronik dan Hamiltonian molekul.

Plot kedua menunjukkan purata penghunian setiap orbit spatial selepas lelaran terakhir. Kita dapat lihat bahawa kedua-dua elektron spin-naik dan spin-turun menghuni lima orbit 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 Circuit kuantum​

Langkah 2: Optimumkan masalah​

Langkah 3: Jalankan eksperimen​

Langkah 4: Proses semula keputusan​

Langkah 1: Petakan masalah kepada Circuit kuantum

Langkah 2: Optimumkan masalah

Langkah 3: Jalankan eksperimen

Langkah 4: Proses semula keputusan