SQD สำหรับการประมาณพลังงานของ Chemistry Hamiltonian

ในบทเรียนนี้ เราจะนำ SQD มาใช้ประมาณ ground state energy ของโมเลกุล

โดยเฉพาะ เราจะอภิปรายหัวข้อต่อไปนี้โดยใช้แนวทาง Qiskit pattern แบบ $4$ ขั้นตอน:

1. ขั้นตอนที่ 1: แมปปัญหาไปยัง quantum circuit และตัวดำเนินการ
   • ตั้งค่า molecular Hamiltonian สำหรับ $N_2$
   • อธิบาย local unitary cluster Jastrow (LUCJ) [1] ที่ได้แรงบันดาลใจจากเคมีและเป็นมิตรกับฮาร์ดแวร์
2. ขั้นตอนที่ 2: เพิ่มประสิทธิภาพสำหรับฮาร์ดแวร์เป้าหมาย
   • เพิ่มประสิทธิภาพจำนวน gate และ layout ของ ansatz สำหรับการรันบนฮาร์ดแวร์
3. ขั้นตอนที่ 3: รันบนฮาร์ดแวร์เป้าหมาย
   • รัน circuit ที่เพิ่มประสิทธิภาพแล้วบน QPU จริงเพื่อสร้างตัวอย่างของ subspace
4. ขั้นตอนที่ 4: Post-process ผลลัพธ์
   • แนะนำ self-consistent configuration recovery loop [2]
     • Post-process ชุด bitstring sample ทั้งหมด โดยใช้ความรู้ล่วงหน้าเกี่ยวกับจำนวนอนุภาคและค่า average orbital occupancy ที่คำนวณในรอบล่าสุด
     • สร้าง batch ของ subsample จาก bitstring ที่กู้คืนแบบ probabilistic
     • ฉายและ diagonalize molecular Hamiltonian บน subspace ที่สุ่มตัวอย่างแต่ละตัว
     • บันทึก ground state energy ต่ำสุดที่พบใน batch ทั้งหมดและอัปเดต avg orbital occupancy

เราจะใช้ software package หลายตัวตลอดบทเรียน

• PySCF เพื่อนิยามโมเลกุลและตั้งค่า Hamiltonian
• ffsim package เพื่อสร้าง LUCJ ansatz
• Qiskit สำหรับ transpile ansatz เพื่อรันบนฮาร์ดแวร์
• Qiskit IBM Runtime เพื่อรัน circuit บน QPU และเก็บตัวอย่าง
• Qiskit addon SQD สำหรับ configuration recovery และการประมาณ ground state energy โดยใช้ subspace projection และ matrix diagonalization

1. แมปปัญหาไปยัง Quantum Circuit และตัวดำเนินการ

Molecular Hamiltonian

Molecular Hamiltonian มีรูปแบบทั่วไปดังนี้:

$$\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}$ คือ fermionic creation/annihilation operator ที่เกี่ยวข้องกับองค์ประกอบ basis set ที่ $p$ และ spin $\sigma$ ส่วน $h_{pr}$ และ $(pr|qs)$ คือ one- และ two-body electronic integral โดยใช้ pySCF เราจะนิยามโมเลกุลและคำนวณ one- และ two-body integral ของ Hamiltonian สำหรับ basis set 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

ในบทเรียนนี้ เราจะใช้การแปลง Jordan-Wigner (JW) เพื่อแมป fermionic wavefunction ไปยัง qubit wavefunction เพื่อให้สามารถเตรียมได้โดยใช้ quantum circuit การแปลง JW แมป Fock space ของ fermion ใน M spatial orbital ไปยัง Hilbert space ของ 2M qubit กล่าวคือ spatial orbital ถูกแบ่งออกเป็น spin orbital สองตัว ตัวหนึ่งเกี่ยวข้องกับอิเล็กตรอน spin up ($\alpha$) และอีกตัวเกี่ยวข้องกับ spin down ($\beta$) spin orbital สามารถถูกครอบครองหรือว่างได้ โดยทั่วไปเมื่อเราพูดถึงจำนวน orbital เราจะใช้จำนวน spatial orbital จำนวน spin orbital จะเป็นสองเท่า ใน quantum circuit เราจะแทน spin orbital แต่ละตัวด้วย qubit หนึ่งตัว ดังนั้น ชุด qubit จะแทน spin-up หรือ $\alpha$-orbital และชุดอื่นจะแทน spin-down หรือ $\beta$-orbital ตัวอย่างเช่น โมเลกุล $N_2$ สำหรับ basis set 6-31g มี $16$ spatial orbital (นั่นคือ $16$ $\alpha$ + $16$ $\beta$ = $32$ spin orbital) ดังนั้น เราต้องการ quantum circuit ขนาด $32$-qubit (เราอาจต้องการ ancilla qubit พิเศษตามที่จะอภิปรายต่อไป) qubit ถูกวัดใน computational basis เพื่อสร้าง bitstring ซึ่งแทน electronic configuration หรือ (Slater) determinant ตลอดบทเรียนนี้ เราจะใช้คำว่า bitstring, configuration และ determinant สลับกัน bitstring บอกเราเกี่ยวกับ electron occupancy ใน spin orbital: $1$ ในตำแหน่ง bit หมายความว่า spin orbital ที่สอดคล้องถูกครอบครอง ในขณะที่ $0$ หมายความว่า spin orbital นั้นว่าง เนื่องจากปัญหา electronic structure เป็น particle preserving มีเพียง spin orbital จำนวนคงที่เท่านั้นที่ต้องถูกครอบครอง โมเลกุล $N_2$ มีอิเล็กตรอน spin-up ($\alpha$) $5$ ตัวและ spin-down ($\beta$) $5$ ตัว ดังนั้น bitstring ใดๆ ที่แทน $\alpha$ และ $\beta$ orbital ต้องมี $1$ ห้าตัวแต่ละส่วนสำหรับโมเลกุล $N_2$

1.1 Quantum Circuit สำหรับการสร้างตัวอย่าง: LUCJ Ansatz

ในบทเรียนนี้ เราจะใช้ local unitary coupled cluster Jastrow (LUCJ) \[1\] ansatz สำหรับการเตรียมสถานะควอนตัมและการสุ่มตัวอย่างในเวลาต่อมา ก่อนอื่น เราจะอธิบาย building block ต่างๆ ของ UCJ ansatz เต็มรูปแบบและการประมาณที่ทำในเวอร์ชัน local จากนั้น โดยใช้ ffsim package เราจะสร้าง LUCJ ansatz และเพิ่มประสิทธิภาพโดยใช้ Qiskit transpiler สำหรับการรันบนฮาร์ดแวร์

UCJ ansatz มีรูปแบบดังต่อไปนี้ (สำหรับผลคูณของ $L$ layer หรือการทำซ้ำของ UCJ operator)

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

โดยที่ $\vert \Phi_{0} \rangle$ คือ reference state ซึ่งโดยทั่วไปถือว่าเป็น Hartree-Fock (HF) state เนื่องจาก Hartree-Fock state นิยามว่ามี orbital ที่มีหมายเลขต่ำสุดถูกครอบครอง การเตรียม HF state จะเกี่ยวข้องกับการใช้ X gate เพื่อตั้ง qubit ที่สอดคล้องกับ orbital ที่ถูกครอบครองเป็นหนึ่ง ตัวอย่างเช่น HF state preparation block สำหรับ spatial orbital $4$ ตัวและ spin up 2 ตัวและ spin down 2 ตัวอาจมีลักษณะดังนี้: การทำซ้ำครั้งเดียวของ UCJ operator ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ ประกอบด้วย diagonal Coulomb evolution ($e^{iJ^{(\mu)}}$) ที่ถูกประกบด้วย orbital rotation ($e^{K^{(\mu)}}$ และ $e^{-K^{(\mu)}}$) Orbital rotation block ทำงานบน spin species เดียว ($\alpha$ (up-spin)/$\beta$ (down-spin)) สำหรับ electron species แต่ละชนิด orbital rotation ประกอบด้วย layer ของ single-qubit $R_{z}$ gate ตามด้วยลำดับของ 2-qubit Given's rotation gate ($XX + YY$ gate)

2-qubit gate ทำงานบน spin-orbital ที่อยู่ติดกัน (qubit ที่เป็นเพื่อนบ้านใกล้เคียงที่สุด) ดังนั้นจึงสามารถนำไปใช้บน IBM® QPU ได้โดยไม่ต้องใช้ SWAP gate $e^{iJ^{(\mu)}}$ หรือที่รู้จักกันในชื่อ diagonal Coulomb operator ประกอบด้วยสาม block สองตัวทำงานบน spin sector เดียวกัน ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ และ $e^{iJ_{\beta \beta}^{(\mu)}}$) และหนึ่งตัวทำงานระหว่าง spin sector ทั้งสอง ($e^{iJ_{\alpha \beta}^{(\mu)}}$)

block ทั้งหมดใน $e^{iJ^{(\mu)}}$ ประกอบด้วย number-number gate $U_{nn}(\phi)$ [1] $U_{nn}(\phi)$ gate สามารถแยกย่อยเพิ่มเติมเป็น $R_{ZZ}(\frac{\phi}{2})$ gate ตามด้วย single-qubit $Rz(-\frac{\phi}{2})$ gate สองตัวที่ทำงานบน qubit แยกกันสองตัว

same-spin component ($J_{\alpha \alpha}$ และ $J_{\beta \beta}$) มี $U_{nn}$ gate ระหว่างคู่ qubit ทุกคู่ที่เป็นไปได้ อย่างไรก็ตาม เนื่องจาก superconducting QPU มีการเชื่อมต่อที่จำกัด qubit จึงต้องถูก swap เพื่อทำ gate ระหว่าง qubit ที่ไม่อยู่ติดกัน

ตัวอย่างเช่น พิจารณา $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (หรือ $e^{iJ_{\beta \beta}^{(\mu)}}$) block ต่อไปนี้สำหรับ spatial orbital $N = 4$ ตัว สำหรับการเชื่อมต่อ qubit แบบ linear gate สามตัวสุดท้ายไม่สามารถนำไปใช้โดยตรงได้เนื่องจากทำงานระหว่าง qubit ที่ไม่อยู่ติดกัน (ตัวอย่างเช่น Q0 และ Q2 ไม่ได้เชื่อมต่อโดยตรง) ดังนั้น เราต้องการ SWAP gate เพื่อทำให้อยู่ติดกัน (รูปต่อไปนี้แสดงตัวอย่างที่มี SWAP gate $3$ ตัว) ต่อไป $J_{\alpha \beta}$ ใช้ gate ระหว่าง orbital ที่มี index เดียวกันจาก spin sector ต่างกัน (ตัวอย่างเช่น ระหว่าง $0\alpha$ และ $0\beta$) เช่นเดียวกัน ถ้า qubit ไม่ได้อยู่ติดกันทางกายภาพบน QPU gate เหล่านี้ก็ต้องการ SWAP เช่นกัน จากการอภิปรายข้างต้น UCJ ansatz เผชิญกับอุปสรรคบางอย่างสำหรับการรันบนฮาร์ดแวร์เนื่องจากต้องใช้ SWAP gate เนื่องจาก qubit ที่ไม่อยู่ติดกัน เวอร์ชัน local ของ UCJ ansatz คือ LUCJ แก้ไขปัญหานี้โดยการลบ $U_{nn}$ บางตัวออกจาก diagonal Coulomb operator

ใน block ที่มี electron species เดียวกัน ($J_{\alpha \alpha}$ และ $J_{\beta \beta}$) เราเก็บเฉพาะ $U_{nn}$ gate ที่เข้ากันได้กับการเชื่อมต่อแบบ nearest-neighbor และลบ gate ระหว่าง qubit ที่ไม่อยู่ติดกันในเวอร์ชัน LUCJ รูปต่อไปนี้แสดง LUCJ block หลังจากลบ gate ที่ไม่อยู่ติดกัน ต่อไป เวอร์ชัน LUCJ ของ $J_{\alpha \beta}$ block ที่ทำงานระหว่าง electron species ต่างกันอาจมีรูปร่างต่างกันตาม device topology

ที่นี่เช่นกัน เวอร์ชัน LUCJ กำจัด gate ที่ไม่เข้ากัน รูปด้านล่างแสดงตัวแปรของ $J_{\alpha \beta}$ block สำหรับ qubit topology ต่างๆ รวมถึง grid, hexagonal, heavy-hex และ linear

• Square: เราสามารถมี $U_{nn}$ gate ระหว่าง $\alpha$ และ $\beta$ orbital ทั้งหมดโดยไม่ต้องใช้ SWAP ดังนั้นจึงไม่จำเป็นต้องลบ $U_{nn}$ gate ใดๆ
• Heavy-hex: $\alpha$-$\beta$ interaction ถูกเก็บไว้ระหว่าง orbital ที่มี index ทุก $4$-th (เช่น 0th, 4th และ 8th) และเป็นแบบ ancilla mediated กล่าวคือ เราต้องการ ancilla qubit ระหว่าง linear chain ที่แทน $\alpha$ และ $\beta$ orbital การจัดเรียงนี้ต้องการ SWAP จำนวนจำกัด
• Hexagonal: orbital ทุกอื่นๆ เช่น 0th, 2nd และ 4th indexed orbital จะกลายเป็น nearest neighbor เมื่อ $\alpha$ และ $\beta$ ถูกวาง layout เป็น linear chain ที่อยู่ติดกันสองเส้น
• Linear: มีเพียง $\alpha$ หนึ่งตัวและ $\beta$ หนึ่งตัวที่เชื่อมต่อกัน ซึ่งหมายความว่า $J_{\alpha \beta}$ block จะมีเพียง gate เดียว แม้ว่าการลบ gate จาก UCJ ansatz เพื่อสร้างเวอร์ชัน LUCJ จะทำให้เข้ากันได้กับฮาร์ดแวร์มากขึ้น แต่ ansatz ก็สูญเสีย expressivity บางส่วน ดังนั้น อาจต้องการการทำซ้ำมากขึ้น ($L$) ของ UCJ operator ที่แก้ไขแล้วเมื่อใช้ LUCJ ansatz

1.2 การกำหนดค่าเริ่มต้น LUCJ Ansatz

LUCJ เป็น ansatz แบบ parameterized และเราต้องกำหนดค่าเริ่มต้นพารามิเตอร์ก่อนการรันบนฮาร์ดแวร์ วิธีหนึ่งในการกำหนดค่าเริ่มต้น ansatz คือการใช้ t1 และ t2 amplitude จากวิธี classical coupled cluster singles and doubles (CCSD) โดยที่ t1 amplitude เป็น coefficient ของ single excitation operator และ t2 amplitude เป็นของ double excitation operator

หมายเหตุว่าแม้การกำหนดค่าเริ่มต้น LUCJ ansatz ด้วย t1 และ t2 amplitude จะให้ผลลัพธ์ที่ดีพอสมควร แต่พารามิเตอร์ 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]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 การสร้าง LUCJ Ansatz โดยใช้ ffsim

เราจะใช้ ffsim package เพื่อสร้างและกำหนดค่าเริ่มต้น ansatz ด้วย t1 และ t2 amplitude ที่คำนวณข้างต้น เนื่องจากโมเลกุลของเรามี closed-shell Hartree-Fock state เราจะใช้เวอร์ชัน spin-balanced ของ UCJ ansatz คือ UCJOpSpinBalanced

เนื่องจากฮาร์ดแวร์ IBM มี heavy-hex topology เราจะนำ zig-zag pattern ที่ใช้ใน [1] และอธิบายข้างต้นสำหรับ qubit interaction มาใช้ ในรูปแบบนี้ orbital (qubit) ที่มี spin เดียวกันจะเชื่อมต่อด้วย line topology (วงกลมสีแดงและน้ำเงิน) เนื่องจาก heavy-hex topology orbital สำหรับ spin ต่างกันมีการเชื่อมต่อระหว่าง orbital ทุก 4th ตัว นั่นคือ 0th, 4th, 8th และต่อๆ ไป (วงกลมสีม่วง)

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)

LUCJ ansatz ที่มี layer ซ้ำกันสามารถเพิ่มประสิทธิภาพได้โดยการรวม block ที่อยู่ติดกันบางตัวเข้าด้วยกัน พิจารณากรณีสำหรับ n_reps=2 orbital rotation block สองตัวตรงกลางสามารถรวมเป็น orbital rotation block เดียว ffsim package มี pass manager ชื่อ ffsim.qiskit.PRE_INIT เพื่อเพิ่มประสิทธิภาพ circuit โดยการรวม block ที่อยู่ติดกันดังกล่าว

2. เพิ่มประสิทธิภาพสำหรับฮาร์ดแวร์เป้าหมาย

ก่อนอื่น เราดึง backend ที่ต้องการ เราจะเพิ่มประสิทธิภาพ circuit สำหรับ backend และจากนั้นรัน circuit ที่เพิ่มประสิทธิภาพแล้วบน backend เดียวกันเพื่อสร้างตัวอย่างสำหรับ subspace

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")

ต่อไป เราแนะนำขั้นตอนต่อไปนี้เพื่อเพิ่มประสิทธิภาพ ansatz และทำให้เข้ากันได้กับฮาร์ดแวร์

• เลือก physical qubit (initial_layout) จากฮาร์ดแวร์เป้าหมายที่ยึดตาม zig-zag pattern (linear chain สองเส้นที่มี ancilla qubit คั่นกลาง) ที่อธิบายข้างต้น การวาง qubit ในรูปแบบนี้นำไปสู่ circuit ที่เข้ากันได้กับฮาร์ดแวร์อย่างมีประสิทธิภาพโดยมี gate น้อยลง
• สร้าง staged pass manager โดยใช้ฟังก์ชัน generate_preset_pass_manager จาก Qiskit พร้อม backend และ initial_layout ที่ต้องการ
• ตั้งค่า pre_init stage ของ staged pass manager เป็น ffsim.qiskit.PRE_INIT ffsim.qiskit.PRE_INIT รวม Qiskit transpiler pass ที่ decompose gate เป็น orbital rotation แล้วรวม orbital rotation ส่งผลให้มี gate น้อยลงใน circuit สุดท้าย
• รัน pass manager บน circuit ของคุณ

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. รันบนฮาร์ดแวร์เป้าหมาย

หลังจากเพิ่มประสิทธิภาพ circuit สำหรับการรันบนฮาร์ดแวร์แล้ว เราพร้อมที่จะรันบนฮาร์ดแวร์เป้าหมายและเก็บตัวอย่างสำหรับการประมาณ ground state energy เนื่องจากเรามีเพียง circuit เดียว เราจะใช้ Job execution mode ของ Qiskit Runtime และรัน circuit ของเรา

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. Post-process ผลลัพธ์

ส่วน post-processing ของ SQD workflow สามารถสรุปได้โดยใช้แผนภาพต่อไปนี้

การสุ่มตัวอย่าง LUCJ ansatz ใน computational basis สร้าง pool ของ noisy configuration $\tilde{\mathcal{\chi}}$ ซึ่งใช้ใน post-processing routine เกี่ยวข้องกับวิธีที่เรียกว่า configuration recovery (อภิปรายรายละเอียดในภายหลัง) เพื่อแก้ไข configuration ที่มีจำนวน electron ไม่ถูกต้องแบบ probabilistic Configuration ที่มีจำนวน electron ถูกต้องเท่านั้น $\tilde{\mathcal{\chi}}_{R}$ จะถูก subsample และกระจายเป็น batch หลายตัวตามความถี่ที่ปรากฏของแต่ละ configuration ที่ไม่ซ้ำกัน แต่ละ batch ของตัวอย่างนิยาม subspace ($\mathcal{S^{(k)}}$) ต่อไป molecular Hamiltonian, $H$, จะถูกฉายลงบน subspace:

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

Hamiltonian ที่ฉายแต่ละตัว $H_{\mathcal{S}^{(k)}}$ จะถูกส่งเข้า Eigensolver ซึ่งจะ diagonalize เพื่อคำนวณ eigenvalue และ eigenvector เพื่อสร้าง eigenstate ขึ้นใหม่ ในบทเรียนนี้ เราฉายและ diagonalize Hamiltonian โดยใช้ qiskit-addon-sqd package ซึ่งใช้วิธี Davidson จาก PySCF สำหรับ diagonalization

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

จากนั้นเราเก็บ eigenvalue (พลังงาน) ต่ำสุดจาก batch ต่างๆ และคำนวณ average orbital occupancy $\text{n}$ ด้วย ข้อมูล average occupancy ใช้ใน configuration recovery step เพื่อแก้ไข noisy configuration แบบ probabilistic

ต่อไป เราอธิบาย self-consistent configuration recovery loop อย่างละเอียดและแสดงตัวอย่างโค้ดที่เป็นรูปธรรมเพื่อนำขั้นตอนที่กล่าวถึงข้างต้นไปใช้สำหรับประมาณ ground state energy ของ $N_2$ Hamiltonian

4.1 Configuration Recovery: ภาพรวม

แต่ละ bit ใน bitstring (Slater determinant) แทน spin orbital ครึ่งขวาของ bitstring แทน spin-up orbital และครึ่งซ้ายแทน spin-down orbital 1 หมายความว่า orbital ถูกครอบครองโดยอิเล็กตรอน และ 0 หมายความว่า orbital นั้นว่าง เรารู้จำนวน particle ที่ถูกต้อง (ทั้ง up-spin electron และ down-spin electron) ล่วงหน้า สมมติว่าเรามี determinant $x$ ที่มีอิเล็กตรอน $N_x$ ตัว (นั่นคือ มี $1$ จำนวน $N_x$ ใน bitstring) จำนวน particle ที่ถูกต้องคือ $N$ ถ้า $N_x \neq N$ เราก็รู้ว่า bitstring ถูก noise ทำให้เสียหาย self-consistent configuration routine พยายามแก้ไข bitstring โดยการพลิก bit จำนวน $|N_x - N|$ แบบ probabilistic โดยใช้ข้อมูล average orbital occupancy average orbital occupancy ($n$) บอกเราว่า orbital มีแนวโน้มถูกครอบครองโดยอิเล็กตรอนมากแค่ไหน ถ้า $N_x < N$ เรามีอิเล็กตรอนน้อยกว่าและต้องพลิก $0$ บางตัวเป็น $1$ และในทางกลับกัน

ความน่าจะเป็นในการพลิกอาจเป็น $|x[i] - avg\_occupancy[i]|$ สำหรับ spin orbital ที่ i ใน [2] ผู้เขียนใช้ความน่าจะเป็นในการพลิกแบบ weighted โดยใช้ฟังก์ชัน modified ReLU

$$\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}$$

ที่นี่ $h$ นิยามตำแหน่ง "corner" ของฟังก์ชัน ReLU และพารามิเตอร์ $\delta$ นิยามค่าของฟังก์ชัน ReLU ที่ corner สำหรับ $\delta = 0$ $w$ จะกลายเป็นฟังก์ชัน ReLU แท้จริง และสำหรับ $\delta >0$ จะกลายเป็น modified ReLU ในเปเปอร์ ผู้เขียนใช้ $\delta = 0.01$ และ $h =$ จำนวน alpha (หรือ beta) particle/จำนวน alpha (หรือ beta) spin orbital $= N/M$ (filling factor)

average orbital occupancy ($n$) ไม่ทราบล่วงหน้า รอบแรกของการประมาณ ground state เริ่มต้นด้วย configuration ที่มีจำนวน particle ถูกต้องใน spin species ทั้งสองเท่านั้น หลังจากรอบแรก เรามีการประมาณ ground state และใช้การประมาณนั้น เราสามารถสร้างการเดาแรกของ $n$ ได้ การเดา $n$ นี้ใช้กู้คืน configuration รัน iteration ถัดไปของการประมาณ ground state และ refine การเดา $n$ แบบ self-consistent กระบวนการนี้ทำซ้ำจนกว่าจะตรงตาม stopping criterion

พิจารณาตัวอย่างต่อไปนี้สำหรับ $N = 2$ และ $x = |1000\rangle$ ($N_x = 1$) เราต้องพลิก $0$ หนึ่งตัวเป็น $1$ เพื่อแก้ไขจำนวน particle และตัวเลือกคือ 1100, 1010 และ 1001 ตาม probability ของการพลิก หนึ่งในตัวเลือกจะถูกเลือกเป็น recovered configuration (หรือ bitstring ที่มีจำนวน particle ถูกต้อง)

สมมติว่าในรอบแรก เรารัน batch สอง batch และ ground state ที่ประมาณได้จาก batch เหล่านั้นคือ:

$$\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}$$

โดยใช้ computational basis state และ amplitude เราสามารถคำนวณความน่าจะเป็นของ electron occupancy (สรุปสั้นๆ ว่า occupancy) ต่อ spin-orbital (qubit) (หมายเหตุว่า ความน่าจะเป็น = |amplitude|$^2$) ด้านล่างนี้ เราจัดตาราง occupancy ตาม qubit สำหรับแต่ละ bitstring ที่ปรากฏใน estimated ground state และคำนวณ orbital occupancy รวมสำหรับ batch หมายเหตุว่า ตาม Qiskit ordering convention bit ขวาสุดแทน qubit-0 (Q0) และ bit ซ้ายสุดแทน Q3

Occupancy (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

Occupancy (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

Occupancy (ค่าเฉลี่ยระหว่าง 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 (average)	0.49	0.51	0.35	0.65

โดยใช้ average orbital occupancy ที่คำนวณข้างต้น เราสามารถหาความน่าจะเป็นของการพลิกสำหรับ orbital ต่างๆ ใน configuration $x = \vert 1000 \rangle$ เนื่องจาก orbital ที่แทนด้วย Q3 ถูกครอบครองอยู่แล้วและไม่จำเป็นต้องพลิก เราตั้ง p(flip) เป็น $0$ สำหรับ orbital ที่เหลือซึ่งว่าง ความน่าจะเป็นของการพลิกคือ $\vert x[i] - \text{n}[i] \vert$ แต่ละตัว พร้อมกับ p(flip) เราคำนวณ probability weight ที่เกี่ยวข้องกับการพลิกโดยใช้ฟังก์ชัน modified ReLU ที่อธิบายข้างต้นด้วย

ความน่าจะเป็นของการพลิก ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

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

สุดท้าย โดยใช้ weighted probability ข้างต้น เราสามารถพลิกหนึ่งตัวจาก Q2, Q1 และ Q0 ที่ว่างอยู่ ตามค่าข้างต้น Q0 จะถูกพลิกมากที่สุด และ recovered configuration ที่เป็นไปได้คือ $\vert \text{1001} \rangle$ กระบวนการ self-consistent configuration recovery สมบูรณ์สามารถสรุปได้ดังนี้:

รอบแรก: สมมติว่า bitstring (configuration หรือ Slater determinant) ที่สร้างโดยควอนตัมคอมพิวเตอร์ประกอบเป็นชุด $\widetilde{\chi}$ ซึ่งรวมทั้ง configuration ที่มีจำนวน particle ถูกต้อง ($\widetilde{\chi}_{correct}$) และไม่ถูกต้อง ($\widetilde{\chi}_{incorrect}$) ใน spin sector แต่ละตัว

1. Configuration จาก ($\widetilde{\chi}_{correct}$) จะถูกสุ่มตัวอย่างเพื่อสร้าง batch $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ ของเวกเตอร์สำหรับ subspace projection จำนวน batch และตัวอย่างใน batch แต่ละตัวเป็นพารามิเตอร์ที่ผู้ใช้กำหนด ยิ่งมีตัวอย่างมากใน batch แต่ละตัว มิติ subspace ก็ยิ่งใหญ่และการ diagonalization ก็ยิ่งต้องการทรัพยากรการคำนวณมาก ในทางกลับกัน ตัวอย่างที่น้อยเกินไปอาจพลาด ground state support vector และนำไปสู่การประมาณที่ไม่ถูกต้อง
2. รัน eigenstate solver (นั่นคือ การฉายลงบน subspace และ diagonalization) บน batch และหา approximate eigenstate $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$
3. จาก approximate eigenstate สร้างการเดาแรกสำหรับ $n$

รอบถัดไป:

1. โดยใช้ $n$ แก้ไข configuration ที่มีจำนวน particle ผิดใน $\widetilde{\chi}_{incorrect}$ สมมติว่าเราตั้งชื่อว่า $\widetilde{\chi}_{correct\_new}$ จากนั้น $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ ประกอบเป็นชุด configuration ใหม่ที่มีจำนวน particle ถูกต้อง
2. $\widetilde{\chi}_{R}$ ถูกสุ่มตัวอย่างเพื่อสร้าง batch $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$
3. Eigenstate solver รันด้วย batch ใหม่และสร้างการประมาณ ground state ใหม่ $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$
4. จาก approximate eigenstate สร้างการเดาที่ refined ของ $n$
5. ถ้า stopping criterion ไม่ตรงตามเงื่อนไข ให้กลับไปขั้นตอน 2.1

4.2 การประมาณ Ground State

ก่อนอื่น เราจะแปลง counts เป็น bitstring matrix และ probability array สำหรับ post-processing

แต่ละแถวใน matrix แทน bitstring ที่ไม่ซ้ำกันหนึ่งตัว เนื่องจาก qubit ถูก index จากด้านขวาของ bitstring ใน Qiskit คอลัมน์ 0 แทน qubit N-1 และคอลัมน์ N-1 แทน qubit 0 โดยที่ N คือจำนวน qubit

alpha orbital ถูกแทนในช่วง column index (N, N/2] (ครึ่งขวา) และ beta orbital ถูกแทนในช่วง (N/2, 0] (ครึ่งซ้าย)

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)

มีตัวเลือกที่ผู้ใช้ควบคุมได้หลายตัวซึ่งสำคัญสำหรับเทคนิคนี้:

• iterations: จำนวนรอบของ self-consistent configuration recovery
• n_batches: จำนวน batch ของ configuration ที่ใช้โดยการเรียก eigenstate solver ต่างๆ
• samples_per_batch: จำนวน configuration ที่ไม่ซ้ำกันที่จะรวมใน batch แต่ละตัว
• max_davidson_cycles: จำนวนสูงสุดของ Davidson cycle ที่รันโดย eigensolver แต่ละตัว

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 postselect 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 การอภิปรายผลลัพธ์

กราฟแรกแสดงว่าหลังจากไม่กี่รอบ เราประมาณ ground state energy ได้ภายใน ~24 mH (ความแม่นยำทางเคมีโดยทั่วไปยอมรับที่ 1 kcal/mol $\approx$ 1.6 mH) กราฟที่สองแสดง average occupancy ของแต่ละ spatial orbital หลังจาก iteration สุดท้าย เราสามารถเห็นได้ว่าอิเล็กตรอน spin-up และ spin-down ทั้งสองครอบครอง orbital ห้าตัวแรกด้วยความน่าจะเป็นสูงในคำตอบของเรา

แม้ว่า ground state energy ที่ประมาณได้จะดีพอสมควร แต่ยังไม่อยู่ในขีดจำกัดความแม่นยำทางเคมี ($\pm \approx 1.6$ mH) ช่องว่างนี้อาจเกิดจาก subspace dimension ขนาดเล็กที่เราใช้ข้างต้นสำหรับการฉายและ diagonalization เนื่องจากเราใช้ samples_per_batch=500 subspace จึงถูกกางออกด้วย max $500$ เวกเตอร์ ซึ่งขาดเวกเตอร์จาก ground state support การเพิ่มพารามิเตอร์ samples_per_batch ควรปรับปรุงความแม่นยำแลกกับทรัพยากรการคำนวณแบบคลาสสิกและ runtime ที่มากขึ้น

# 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

แบบฝึกหัดสำหรับผู้อ่าน

ค่อยๆ เพิ่มพารามิเตอร์ samples_per_batch (ตัวอย่างเช่น จาก $1000$ ถึง $10000$ ด้วยขั้น $1000$; ขึ้นอยู่กับหน่วยความจำของคอมพิวเตอร์) และเปรียบเทียบ ground state energy ที่ประมาณได้

References

[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.