SQD untuk anggaran tenaga Hamiltonian kimia

Dalam pelajaran ini, kita akan menerapkan SQD untuk menganggar tenaga keadaan tanah sebuah molekul.

Khususnya, kita akan membincangkan topik-topik berikut menggunakan pendekatan corak Qiskit 4 langkah:

1. Langkah 1: Petakan masalah kepada litar kuantum dan operator
   • Sediakan Hamiltonian molekul untuk $N_2$.
   • Terangkan ansatz cluster Jastrow unitari setempat (LUCJ) yang terinspirasi kimia dan mesra perkakasan [1]
2. Langkah 2: Optimumkan untuk perkakasan sasaran
   • Optimumkan bilangan gate dan susun atur ansatz untuk pelaksanaan perkakasan
3. Langkah 3: Laksanakan pada perkakasan sasaran
   • Jalankan litar yang dioptimumkan pada QPU sebenar untuk menjana sampel subruang.
4. Langkah 4: Pasca-proses keputusan
   • Perkenalkan gelung pemulihan konfigurasi konsisten diri [2]
     • Pasca-proses set penuh sampel bitstring, menggunakan pengetahuan terdahulu tentang bilangan zarah dan purata penghunian orbital yang dikira pada iterasi terbaru.
     • Buat batch subsampel secara kebarangkalian daripada bitstring yang dipulihkan.
     • Unjurkan dan diagonalisasikan Hamiltonian molekul ke atas setiap subruang yang disampling.
     • Simpan tenaga keadaan tanah minimum yang dijumpai merentasi semua batch dan kemas kini purata penghunian orbital.

Kita akan menggunakan beberapa pakej perisian sepanjang pelajaran.

• PySCF untuk mendefinisikan molekul dan menyediakan Hamiltonian.
• pakej ffsim untuk membina ansatz LUCJ.
• Qiskit untuk mentranspil ansatz untuk pelaksanaan perkakasan.
• Qiskit IBM Runtime untuk melaksanakan litar pada QPU dan mengumpul sampel.
• Qiskit addon SQD pemulihan konfigurasi dan anggaran tenaga keadaan tanah menggunakan unjuran subruang dan diagonalisasi matriks.

1. Petakan masalah kepada litar kuantum dan operator

Hamiltonian Molekul

Hamiltonian molekul mengambil bentuk generik:

$$\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}$ ialah operator penciptaan/pemusnahan fermionis yang berkaitan dengan elemen set asas ke-$p$ dan spin $\sigma$. $h_{pr}$ dan $(pr|qs)$ ialah kamiran elektronik satu dan dua jasad. Menggunakan pySCF, kita akan mendefinisikan molekul dan mengira kamiran satu dan dua jasad Hamiltonian untuk set asas 6-31g.

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

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# 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)]],  # Two N atoms 1 angstrom apart
    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)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

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

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

Dalam pelajaran ini, kita akan menggunakan transformasi Jordan-Wigner (JW) untuk memetakan fungsi gelombang fermionis kepada fungsi gelombang Qubit supaya ia boleh disediakan menggunakan litar kuantum. Transformasi JW memetakan ruang Fock fermion dalam M orbital ruang ke ruang Hilbert 2M Qubit, iaitu, orbital ruang dibahagikan kepada dua orbital spin, satu yang berkaitan dengan elektron spin atas ($\alpha$) dan satu lagi spin bawah ($\beta$). Orbital spin boleh dihuni atau tidak dihuni. Biasanya, apabila kita merujuk kepada bilangan orbital, kita menggunakan bilangan orbital ruang. Bilangan orbital spin akan menjadi dua kali ganda. Dalam litar kuantum, kita akan mewakili setiap orbital spin dengan satu Qubit. Oleh itu, satu set Qubit akan mewakili orbital spin atas atau $\alpha$, dan set lain akan mewakili orbital spin bawah atau $\beta$. Sebagai contoh, molekul $N_2$ untuk set asas 6-31g mempunyai $16$ orbital ruang (iaitu, $16$ $\alpha$ + $16$ $\beta$ = $32$ orbital spin). Oleh itu, kita memerlukan litar kuantum $32$-Qubit (kita mungkin memerlukan Qubit ancilla tambahan seperti yang dibincangkan kemudian). Qubit diukur dalam asas pengkomputeran untuk menjana bitstring, yang mewakili konfigurasi elektronik atau determinan (Slater). Sepanjang pelajaran ini, kita akan menggunakan istilah bitstring, konfigurasi, dan determinan secara bergantian. Bitstring memberitahu kita penghunian elektron dalam orbital spin: $1$ dalam kedudukan bit bermakna orbital spin yang sepadan dihuni, manakala $0$ bermakna orbital spin kosong. Kerana masalah struktur elektronik memelihara zarah, hanya bilangan tetap orbital spin yang mesti dihuni. Molekul $N_2$ mempunyai $5$ elektron spin atas ($\alpha$) dan $5$ elektron spin bawah ($\beta$). Oleh itu, sebarang bitstring yang mewakili orbital $\alpha$ dan $\beta$ mesti mempunyai lima $1$ masing-masing untuk molekul $N_2$.

1.1 Litar kuantum untuk penjanaan sampel: Ansatz LUCJ

Dalam pelajaran ini, kita akan menggunakan ansatz cluster Jastrow unitari setempat (LUCJ) \[1\] untuk penyediaan keadaan kuantum dan pensampelan berikutnya. Pertama, kita akan menjelaskan blok-blok bangunan yang berbeza bagi ansatz UCJ penuh dan penghampiran yang dibuat dalam versi setempatnya. Seterusnya, dengan menggunakan pakej ffsim, kita akan membina ansatz LUCJ dan mengoptimumkannya menggunakan Transpiler Qiskit untuk pelaksanaan perkakasan.

Ansatz UCJ mempunyai bentuk berikut (untuk produk $L$ lapisan atau ulangan operator UCJ.)

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

di mana, $\vert \Phi_{0} \rangle$ ialah keadaan rujukan, biasanya diambil sebagai keadaan Hartree-Fock (HF). Kerana keadaan Hartree-Fock ditakrifkan sebagai mempunyai orbital bernombor terendah yang dihuni, penyediaan keadaan HF akan melibatkan penerapan gate X untuk menetapkan Qubit yang sepadan dengan orbital yang dihuni kepada satu. Sebagai contoh, blok penyediaan keadaan HF untuk 4 orbital ruang dan 2 spin atas dan 2 spin bawah mungkin kelihatan seperti berikut: Satu ulangan operator UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ terdiri daripada evolusi Coulomb pepenjuru ($e^{iJ^{(\mu)}}$) yang dikelilingi oleh putaran orbital ($e^{K^{(\mu)}}$ dan $e^{-K^{(\mu)}}$). Blok putaran orbital berfungsi pada satu spesies spin ($\alpha$ (spin atas)/$\beta$ (spin bawah)). Untuk setiap spesies elektron, putaran orbital terdiri daripada lapisan gate $R_{z}$ satu-qubit diikuti oleh urutan gate putaran Given's 2-qubit ($XX + YY$ gate).

Gate 2-qubit bertindak pada orbital spin bersebelahan (Qubit jiran terdekat), dan oleh itu, boleh dilaksanakan pada QPU IBM® tanpa memerlukan gate SWAP. $e^{iJ^{(\mu)}}$, juga dikenali sebagai operator Coulomb pepenjuru, terdiri daripada tiga blok. Dua daripadanya berfungsi pada sektor spin yang sama ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ dan $e^{iJ_{\beta \beta}^{(\mu)}}$), dan satu berfungsi antara dua sektor spin ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Semua blok dalam $e^{iJ^{(\mu)}}$ terdiri daripada gate nombor-nombor $U_{nn}(\phi)$ [1]. Gate $U_{nn}(\phi)$ boleh diurai lebih lanjut kepada gate $R_{ZZ}(\frac{\phi}{2})$ diikuti oleh dua gate $Rz(-\frac{\phi}{2})$ satu-qubit yang bertindak pada dua Qubit berasingan.

Komponen spin yang sama ($J_{\alpha \alpha}$ dan $J_{\beta \beta}$) mempunyai gate $U_{nn}$ antara semua pasangan Qubit yang mungkin. Namun, kerana QPU superkonduktor mempunyai sambungan yang terhad, Qubit mesti ditukar untuk merealisasikan gate antara Qubit yang tidak bersebelahan.

Sebagai contoh, pertimbangkan blok $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (atau $e^{iJ_{\beta \beta}^{(\mu)}}$) berikut untuk $N = 4$ orbital ruang. Untuk sambungan Qubit linear, tiga gate terakhir tidak boleh dilaksanakan secara langsung kerana mereka berfungsi antara Qubit yang tidak bersebelahan (sebagai contoh, Q0 dan Q2 tidak bersambung secara langsung). Oleh itu, kita memerlukan gate SWAP untuk menjadikannya bersebelahan (rajah berikut menunjukkan contoh dengan $3$ gate SWAP). Seterusnya, $J_{\alpha \beta}$ melaksanakan gate antara orbital berindeks sama dari sektor spin yang berbeza (sebagai contoh, antara $0\alpha$ dan $0\beta$). Begitu juga, jika Qubit tidak bersampingan secara fizikal pada QPU, gate ini juga akan memerlukan SWAP. Daripada perbincangan di atas, ansatz UCJ menghadapi beberapa cabaran untuk pelaksanaan perkakasan kerana ia memerlukan gate SWAP disebabkan interaksi Qubit yang tidak bersebelahan. Varian setempat ansatz UCJ, LUCJ, menangani cabaran ini dengan mengeluarkan beberapa $U_{nn}$ daripada operator Coulomb pepenjuru.

Dalam blok spesies elektron yang sama, $J_{\alpha \alpha}$ dan $J_{\beta \beta}$), kita hanya menyimpan gate $U_{nn}$ yang serasi dengan sambungan jiran terdekat dan mengeluarkan gate antara Qubit yang tidak bersebelahan dalam versi LUCJ. Rajah berikut menunjukkan blok LUCJ selepas penyingkiran gate yang tidak bersebelahan. Seterusnya, versi LUCJ blok $J_{\alpha \beta}$ yang berfungsi antara spesies elektron yang berbeza boleh mengambil bentuk berbeza berdasarkan topologi peranti.

Di sini juga, versi LUCJ menyingkirkan gate yang tidak serasi. Rajah di bawah menunjukkan varian blok $J_{\alpha \beta}$ untuk topologi Qubit yang berbeza termasuk grid, heksagonal, heavy-hex, dan linear.

• Petak: kita boleh mempunyai gate $U_{nn}$ antara semua orbital $\alpha$ dan $\beta$ tanpa sebarang SWAP, dan oleh itu, tidak perlu mengeluarkan sebarang gate $U_{nn}$.
• Heavy-hex: Interaksi $\alpha$-$\beta$ dikekalkan antara setiap orbital berindeks ke-$4$ (seperti ke-0, ke-4, dan ke-8) dan dimediasi ancilla, iaitu, kita memerlukan Qubit ancilla antara rantaian linear yang mewakili orbital $\alpha$ dan $\beta$. Susunan ini memerlukan bilangan SWAP yang terhad.
• Heksagonal: Setiap orbital yang lain, seperti orbital berindeks ke-0, ke-2, dan ke-4, menjadi jiran terdekat apabila $\alpha$ dan $\beta$ disusun dalam dua rantaian linear bersebelahan.
• Linear: Hanya satu orbital $\alpha$ dan satu $\beta$ bersambung, bermakna blok $J_{\alpha \beta}$ hanya akan mempunyai satu gate. Walaupun mengeluarkan gate daripada ansatz UCJ untuk membina versi LUCJ menjadikannya lebih serasi dengan perkakasan, ansatz kehilangan beberapa kebolehunyahan. Oleh itu, lebih banyak ulangan ($L$) operator UCJ yang dimodifikasi mungkin diperlukan apabila menggunakan ansatz LUCJ.

1.2 Permulaan ansatz LUCJ

LUCJ ialah ansatz berparameter, dan kita perlu memulakan parameter sebelum pelaksanaan perkakasan. Satu cara untuk memulakan ansatz ialah dengan menggunakan amplitud t1 dan t2 daripada kaedah cluster tunggal dan ganda berpasangan (CCSD) klasik, di mana amplitud t1 ialah pekali operator rangsangan tunggal dan amplitud t2 adalah untuk operator rangsangan ganda.

Perhatikan bahawa walaupun memulakan ansatz LUCJ dengan amplitud t1 dan t2 menghasilkan keputusan yang baik, parameter ansatz mungkin memerlukan pengoptimuman lebih lanjut.

# 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]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Membina ansatz LUCJ menggunakan ffsim

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

Kerana perkakasan IBM mempunyai topologi heavy-hex, kita akan menggunakan corak zig-zag yang digunakan dalam [1] dan dijelaskan di atas untuk interaksi Qubit. Dalam corak ini, orbital (Qubit) dengan spin yang sama disambungkan dengan topologi garis (bulatan merah dan biru). Disebabkan topologi heavy-hex, orbital untuk spin yang berbeza mempunyai sambungan antara setiap orbital ke-4, iaitu, ke-0, ke-4, ke-8, dan seterusnya (bulatan ungu).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
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()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

Ansatz LUCJ dengan lapisan berulang boleh dioptimumkan dengan menggabungkan beberapa blok bersebelahan. Pertimbangkan kes untuk n_reps=2. Dua blok putaran orbital di tengah boleh digabungkan menjadi satu blok putaran orbital tunggal. Pakej ffsim mempunyai pengurus laluan bernama ffsim.qiskit.PRE_INIT untuk mengoptimumkan litar dengan menggabungkan blok bersebelahan sedemikian.

2. Optimumkan untuk perkakasan sasaran

Pertama, kita ambil Backend pilihan kita. Kita akan mengoptimumkan litar kita untuk Backend tersebut, dan kemudian melaksanakan litar yang dioptimumkan pada Backend yang sama untuk menjana sampel bagi subruang.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Seterusnya, kita cadangkan langkah-langkah berikut untuk mengoptimumkan ansatz dan menjadikannya serasi dengan perkakasan.

• Pilih Qubit fizikal (initial_layout) dari perkakasan sasaran yang mengikut corak zig-zag (dua rantaian linear dengan Qubit ancilla di antara mereka) yang diterangkan di atas. Menyusun Qubit dalam corak ini menghasilkan litar yang cekap dan serasi perkakasan dengan lebih sedikit gate.
• Jana pengurus laluan berperingkat menggunakan fungsi generate_preset_pass_manager dari Qiskit dengan backend dan initial_layout pilihan anda.
• Tetapkan peringkat pre_init pengurus laluan berperingkat anda kepada ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT merangkumi laluan Transpiler Qiskit yang mengurai gate kepada putaran orbital dan kemudian menggabungkan putaran orbital, menghasilkan lebih sedikit gate 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': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Laksanakan pada perkakasan sasaran

Setelah mengoptimumkan litar untuk pelaksanaan perkakasan, kita bersedia untuk menjalankannya pada perkakasan sasaran dan mengumpul sampel untuk anggaran tenaga keadaan tanah. Kerana kita hanya mempunyai satu litar, kita akan menggunakan mod pelaksanaan Job Qiskit Runtime dan melaksanakan litar kita.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Pasca-proses keputusan

Bahagian pasca-pemprosesan aliran kerja SQD boleh diringkaskan menggunakan gambar rajah berikut.

Pensampelan ansatz LUCJ dalam asas pengkomputeran menjana kolam konfigurasi berhingar $\tilde{\mathcal{\chi}}$, yang digunakan dalam rutin pasca-pemprosesan. Ia melibatkan kaedah yang dipanggil (perincian dibincangkan kemudian) pemulihan konfigurasi untuk membetulkan konfigurasi dengan bilangan elektron yang salah secara kebarangkalian. Konfigurasi hanya dengan bilangan elektron yang betul $\tilde{\mathcal{\chi}}_{R}$ kemudian disubsampling dan diagihkan ke dalam pelbagai batch berdasarkan kekerapan kemunculan setiap konfigurasi unik. Setiap batch sampel mendefinisikan subruang ($\mathcal{S^{(k)}}$). Seterusnya, Hamiltonian molekul, $H$, diunjurkan ke atas subruang:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ dengan } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

Setiap Hamiltonian yang diunjurkan $H_{\mathcal{S}^{(k)}}$ kemudian dimasukkan ke dalam Penyelesai Eigen, di mana ia didiagonalisasikan untuk mengira nilai eigen dan vektor eigen bagi membina semula keadaan eigen. Dalam pelajaran ini, kita mengunjurkan dan mendiagonalisasikan Hamiltonian menggunakan pakej qiskit-addon-sqd yang menggunakan kaedah Davidson dari PySCF untuk diagonalisasi.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

Kita kemudian mengumpul nilai eigen terendah (tenaga) dari batch, dan juga mengira purata penghunian orbital, $\text{n}$. Maklumat penghunian purata digunakan dalam langkah pemulihan konfigurasi untuk membetulkan konfigurasi berhingar secara kebarangkalian.

Seterusnya, kita menjelaskan gelung pemulihan konfigurasi konsisten diri secara terperinci dan menunjukkan contoh kod konkrit untuk melaksanakan langkah-langkah yang disebutkan di atas bagi menganggar tenaga keadaan tanah Hamiltonian $N_2$.

4.1 Pemulihan konfigurasi: gambaran keseluruhan

Setiap bit dalam bitstring (determinan Slater) mewakili orbital spin. Separuh kanan bitstring mewakili orbital spin atas, dan separuh kiri mewakili orbital spin bawah. 1 bermakna orbital dihuni oleh elektron, dan 0 bermakna orbital kosong. Kita tahu bilangan zarah yang betul (kedua-dua elektron spin atas dan elektron spin bawah) terlebih dahulu. Andaikan kita mempunyai determinan $x$ dengan $N_x$ elektron (iaitu, terdapat $N_x$ nombor 1 dalam bitstring) di dalamnya. Bilangan zarah yang betul ialah $N$. Jika $N_x \neq N$, maka kita tahu bahawa bitstring dirosakkan oleh hingar. Rutin konfigurasi konsisten diri cuba membetulkan bitstring dengan membalikkan $|N_x - N|$ bit secara kebarangkalian dengan memanfaatkan maklumat penghunian orbital purata. Penghunian orbital purata ($n$) memberitahu kita seberapa besar kemungkinan orbital dihuni oleh elektron. Jika $N_x < N$, kita mempunyai lebih sedikit elektron dan perlu membalikkan beberapa 0 kepada 1 dan sebaliknya.

Kebarangkalian membalikkan boleh menjadi $|x[i] - avg\_occupancy[i]|$ untuk orbital spin ke-i. Dalam [2], penulis menggunakan kebarangkalian membalikkan berbobot menggunakan fungsi ReLU yang dimodifikasi.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

Di sini $h$ mentakrifkan lokasi "sudut" fungsi ReLU, dan parameter $\delta$ mentakrifkan nilai fungsi ReLU pada sudut. Untuk $\delta = 0$, $w$ menjadi fungsi ReLU sebenar, dan untuk $\delta >0$, ia menjadi ReLU yang dimodifikasi. Dalam makalah tersebut, penulis menggunakan $\delta = 0.01$ dan $h =$ bilangan zarah alfa (atau beta)/bilangan orbital spin alfa (atau beta) $= N/M$ (faktor pengisian).

Penghunian orbital purata ($n$) tidak diketahui terlebih dahulu. Iterasi pertama anggaran keadaan tanah bermula dengan konfigurasi yang hanya mempunyai bilangan zarah yang betul dalam kedua-dua spesies spin. Selepas iterasi pertama, kita mempunyai anggaran keadaan tanah, dan menggunakan anggaran tersebut, kita boleh membina tekaan pertama $n$. Tekaan $n$ ini digunakan untuk memulihkan konfigurasi, menjalankan iterasi anggaran keadaan tanah seterusnya, dan menghalusi tekaan $n$ secara konsisten diri. Proses ini berulang sehingga kriteria berhenti dipenuhi.

Pertimbangkan contoh berikut untuk $N = 2$ dan $x = |1000\rangle$ ($N_x = 1$). Kita perlu membalikkan salah satu daripada 0 kepada 1 untuk membetulkannya untuk bilangan zarah, dan pilihannya ialah 1100, 1010, dan 1001. Berdasarkan kebarangkalian membalikkan, salah satu daripada pilihan tersebut akan dipilih sebagai konfigurasi yang dipulihkan (atau bitstring dengan bilangan zarah yang betul).

Andaikan dalam iterasi pertama kita menjalankan dua batch, dan keadaan tanah yang dianggar daripada mereka ialah:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

Menggunakan keadaan asas pengkomputeran dan amplitudnya, kita boleh mengira kebarangkalian penghunian elektron (ringkasnya penghunian) per orbital spin (Qubit) (perhatikan bahawa kebarangkalian = |amplitud|$^2$). Di bawah kita menyusun penghunian mengikut Qubit untuk setiap bitstring yang muncul dalam keadaan tanah yang dianggar dan mengira penghunian orbital keseluruhan untuk setiap batch. Perhatikan bahawa, menurut konvensyen susunan Qiskit, bit paling kanan mewakili qubit-0 (Q0), dan bit paling kiri mewakili Q3.

Penghunian (Batch0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

Penghunian (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

Penghunian (purata merentasi batch)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (purata)	0.49	0.51	0.35	0.65

Menggunakan penghunian orbital purata yang dikira di atas, kita boleh mencari kebarangkalian membalikkan untuk orbital berbeza dalam konfigurasi $x = \vert 1000 \rangle$. Kerana orbital yang diwakili oleh Q3 sudah dihuni dan tidak perlu dibalikkan, kita tetapkan p(balik) kepada $0$. Untuk orbital yang tidak dihuni yang selebihnya, kebarangkalian membalikkan ialah $\vert x[i] - \text{n}[i] \vert$ masing-masing. Bersama-sama dengan p(balik), kita juga mengira pemberat kebarangkalian yang berkaitan dengan membalikkan menggunakan fungsi ReLU yang dimodifikasi yang diterangkan di atas.

Kebarangkalian membalikkan ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

	Q3	Q2	Q1	Q0
p(balik) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(balik))	0	0.03	0.007	0.31

Akhirnya, menggunakan kebarangkalian berbobot di atas, kita boleh membalikkan salah satu daripada orbital Q2, Q1, dan Q0 yang tidak dihuni. Berdasarkan nilai di atas, Q0 paling berkemungkinan dibalikkan, dan konfigurasi yang dipulihkan yang mungkin ialah $\vert \text{1001} \rangle$. Proses pemulihan konfigurasi konsisten diri yang lengkap boleh diringkaskan seperti berikut:

Iterasi pertama: Andaikan bitstring (konfigurasi atau determinan Slater) yang dijana oleh komputer kuantum membentuk set $\widetilde{\chi}$, yang merangkumi kedua-dua konfigurasi dengan bilangan zarah yang betul ($\widetilde{\chi}_{correct}$) dan tidak betul ($\widetilde{\chi}_{incorrect}$) dalam setiap sektor spin.

1. Konfigurasi dari ($\widetilde{\chi}_{correct}$) disampling secara rawak untuk mencipta batch $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ vektor untuk unjuran subruang. Bilangan batch dan sampel dalam setiap batch adalah parameter yang ditentukan pengguna. Semakin besar bilangan sampel dalam setiap batch, semakin besar dimensi subruang dan semakin menuntut pengiraan diagonalisasinya. Di sisi lain, bilangan sampel yang terlalu kecil mungkin terlepas vektor sokongan keadaan tanah dan mengakibatkan anggaran yang salah.
2. Jalankan penyelesai keadaan eigen (iaitu, unjuran ke atas subruang dan diagonalisasi) pada batch dan dapatkan keadaan eigen anggaran. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. Daripada keadaan eigen anggaran bina tekaan pertama untuk $n$.

Iterasi berikutnya:

1. Menggunakan $n$ betulkan konfigurasi dengan bilangan zarah yang salah dalam $\widetilde{\chi}_{incorrect}$. Andaikan kita menamakannya $\widetilde{\chi}_{correct\_new}$. Kemudian, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ membentuk set konfigurasi baharu dengan bilangan zarah yang betul.
2. $\widetilde{\chi}_{R}$ disampling untuk mencipta batch $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Penyelesai keadaan eigen dijalankan dengan batch baharu dan menjana anggaran keadaan tanah baharu $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Daripada keadaan eigen anggaran bina tekaan yang dihalusi untuk $n$.
5. Jika kriteria berhenti tidak dipenuhi, kembali ke langkah 2.1.

4.2 Anggaran keadaan tanah

Pertama, kita akan mengubah kiraan kepada matriks bitstring dan array kebarangkalian untuk pasca-pemprosesan.

Setiap baris dalam matriks mewakili satu bitstring unik. Kerana Qubit diindeks dari kanan bitstring dalam Qiskit, lajur 0 mewakili qubit N-1, dan lajur N-1 mewakili qubit 0, di mana N ialah bilangan Qubit.

Orbital alfa diwakili dalam julat indeks lajur (N, N/2] (separuh kanan), dan orbital beta diwakili dalam julat lajur (N/2, 0] (separuh kiri).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Terdapat beberapa pilihan yang dikawal pengguna yang penting untuk teknik ini:

• iterations: Bilangan iterasi pemulihan konfigurasi konsisten diri
• n_batches: Bilangan batch konfigurasi yang digunakan oleh panggilan berbeza kepada penyelesai keadaan eigen
• samples_per_batch: Bilangan konfigurasi unik untuk disertakan dalam setiap batch
• max_davidson_cycles: Bilangan kitaran Davidson maksimum yang dijalankan oleh setiap penyelesai eigen

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We post-select here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_davidson=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Perbincangan keputusan

Plot pertama menunjukkan bahawa selepas beberapa iterasi kita menganggar tenaga keadaan tanah dalam ~24 mH (ketepatan kimia biasanya diterima sebagai 1 kcal/mol $\approx$ 1.6 mH). Plot kedua menunjukkan purata penghunian setiap orbital ruang selepas iterasi terakhir. Kita dapat melihat bahawa kedua-dua elektron spin atas dan spin bawah menghuni lima orbital pertama dengan kebarangkalian tinggi dalam penyelesaian kita.

Walaupun tenaga keadaan tanah yang dianggar adalah baik, ia tidak dalam had ketepatan kimia ($\pm \approx 1.6$ mH). Jurang ini boleh dikaitkan dengan dimensi subruang kecil yang kita gunakan di atas untuk unjuran dan diagonalisasi. Kerana kita menggunakan samples_per_batch=500, subruang dicakupi oleh maksimum $500$ vektor, yang ketinggalan vektor dari sokongan keadaan tanah. Meningkatkan parameter samples_per_batch harus meningkatkan ketepatan dengan mengorbankan lebih banyak sumber pengkomputeran klasik dan masa jalan.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

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

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

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-6)
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.02234 Ha
Absolute error: 0.02434 Ha

Latihan untuk pembaca

Tingkatkan parameter samples_per_batch secara beransur-ansur (sebagai contoh, dari $1000$ hingga $10000$ dengan langkah $1000$; bergantung pada memori komputer anda) dan bandingkan tenaga keadaan tanah yang dianggar.

Rujukan

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.