SqDRIFT algorithm para sa ground state estimation

Tinatayang paggamit: 180 seconds sa isang Heron r3 processor (TANDAAN: Isa lamang itong tantiya. Maaaring mag-iba ang runtime mo.)

Hinahanap ang bersyon ng C++?

Gumagamit ang tutorial na ito ng Python. Para sa C++ implementation, kasama ang source code at build instructions, tingnan ang C++ SqDRIFT tutorial.

Mga natutunan

• Matutunan kung paano gumawa ng mas mababang depth na circuits kumpara sa Trotterization

• Dumaan sa isang end-to-end workflow para sa ground state estimation gamit ang qDRIFT at SQD

• Matutunan kung paano gamitin ang qiskit-fermions kasabay ng ibang Qiskit addons para ipatupad ang ganitong workflow

Ipinapakita ang tutorial na ito bilang isang Python notebook para sa layuning pagtuturo.

Mga kinakailangan

• Basahin ang overview ng Sample-based quantum diagonalization (SQD)

• Basahin ang lesson na Sample-based Krylov Quantum Diagonalization (SKQD)

Kaligiran

Ang SqDRIFT ay isang variant ng SKQD na pumapalit sa pangangailangang pumili ng ansatz kung saan kukuha ng sample ng bitstrings gamit ang isang ensemble ng time-evolution circuits na direktang binuo mula sa target Hamiltonian. Nakakamit ito sa pamamagitan ng pag-subsample ng mas maliliit na time-evolution operators mula sa Hamiltonian batay sa mga coefficient nito, na kilala bilang ang qDRIFT Trotterization method.

Ginagamit ng tutorial na ito ang Qiskit Fermions para gumawa ng mas natural na fermionic circuits para sa algorithm na qDRIFT, na sinusundan ng paggamit ng fermionic layout at synthesis passes bago i-plug ang mga circuit sa tradisyunal na Qiskit pipeline para sa hardware execution.

Hayaang ang Hamiltonian ay sa anyong:

$$H = \sum_{i=1}^{N} c_i h_i$$

kung saan, nang hindi nawawalan ng generality, kinakailangan natin na $c_i > 0$ at ang pinakamalaking eigenvalue ng $h_i$ ay pantay, sa absolute value, sa $1$. Anumang may-sign o complex na prefactor ay naiuukol sa $h_i$, kaya ang mga coefficient na $c_i$ ay mahigpit na positive weights habang ang $h_i$ ay may dalang direksyon ng bawat term. Dito, ang $N$ ay ang bilang ng terms (o, pagkatapos ng grouping, ang bilang ng groups) sa Hamiltonian; ito ay property ng Hamiltonian at naiiba sa bilang ng operators na na-sample sa iisang circuit, na isinusulat bilang $n$ sa ibaba.

Para sa target time na $t$, ipinapatupad ng qDRIFT algorithm ang isang operator na $V_k$, kung saan ang $k$ ay mula $1 \cdots K$ at kumakatawan sa $k_{th}$ SqDRIFT circuit, na tinutukoy bilang:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Dito ang $n$ ay ang bilang ng sampled operators bawat circuit at ang $K$ ay ang bilang ng circuits sa ensemble. Ang produkto ay tumatakbo sa $n$ na pagkuha, hindi sa lahat ng $N$ Hamiltonian terms, at dahil kinukuha ang mga term nang may replacement, maaaring lumabas nang higit sa isang beses ang parehong $h_i$ sa iisang $V_k$.

Ang dami:

$$\lambda = \sum_{i=1}^{N} c_i$$

ay ang $L_1$ norm ng mga coefficient, kaya bawat isa sa $n$ na hakbang ay umuunlad sa parehong tagal na $\lambda t / n$ anuman ang term na nakuha. Ang pagkakapareho ng step angle ang katangiang tampok ng qDRIFT: ang isang coefficient ay nakaaapekto sa resulta sa pamamagitan ng gaano kadalas na nakukuha ang term nito, hindi sa pamamagitan ng kung gaano kalayo umiikot ang term na iyon. Ang mga index ay kinukuha mula sa distribution:

$$P[k_i] = \frac{c_i}{\lambda}$$

kaya ang series na $(k_1, \ldots, k_n)$ ay isang random sequence ng term indices na kinuha mula sa distribution na ito. Dahil positive ang mga $c_i$ at ang kabuuan nila ay $\lambda$, isa itong normalized probability distribution, at ang expectation ng resultang channel sa mga random draws ay tinatantya ang evolution sa ilalim ng $H$, na may error na bumababa habang lumalaki ang $n$. Tandaan na ang approximation error ay nakadepende sa $\lambda$ sa halip na sa bilang ng terms na $N$.

(Isinusulat ng SqDRIFT paper ang bilang ng terms bilang $\mathcal{N}$ at ang haba ng sequence bilang $N$; ginagamit natin ang $N$ at $n$ dito para panatilihing malinaw ang pagkakaiba ng dalawa.)

Ipinapakita ng tutorial na ito kung paano bumuo ng ensemble ng ganitong mga randomized circuits. Matapos nating gawin ang mga circuit na ito, katulad ng paggawa natin ng Krylov subspace para sa iba't ibang operators, kumukuha tayo ng sample ng bitstrings mula sa maraming ganitong operators na may iba't ibang time parameters. Sinisiguro nito ang mas mataas na overlap sa pagitan ng ground state vectors at ng sampled bitstrings.

Mga kailangan

Bago simulan ang tutorial na ito, siguraduhing na-install mo na ang

• Isang Python (>=3.10) virtual environment
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Tandaan na plural ang pangalan)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Puwede mong i-install ang lahat ng kinakailangang packages gamit ang:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Pag-setup

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Halimbawa ng simulator

Hakbang 1: I-map ang classical inputs papunta sa quantum problem

Pagbabasa at Paghahanda ng FCIDump

Para sa tutorial na ito, ilo-load natin ang electronic structure Hamiltonian para sa nitrogen (N2). May iba pang paraan din para gumawa ng fermionic operators. Sumangguni sa dokumentasyon sa qiskit_fermions.operators.library.

Tungkol sa FCIDump na ito. Inilalarawan ng file na N2_sto_3g ang isang molekula ng nitrogen ($N_2$) sa minimal STO-3G basis sa interatomic separation na 1.09 $\AA$, ang experimental equilibrium bond length. Idineklara ng header nito ang NORB=10, NELEC=14, at MS2=0: 10 spatial orbitals (kaya 20 spin orbitals, at 20 qubits sa ilalim ng Jordan-Wigner), 14 electrons sa isang spin singlet, kaya pito $\alpha$ at pitong $\beta$ electrons. Lahat ng orbitals ay binigyan ng symmetry label 1, ibig sabihin, walang point-group symmetry na ginagamit. Dahil isa itong full-space STO-3G dump, walang orbitals na na-freeze at maliit lang ang correlation space kaya maaaring kalkulahin nang classical ang eksaktong FCI reference energy para sa paghahambing, tulad ng ipinapakita sa susunod na cell.

Maaaring i-regenerate ang isang katumbas na file gamit ang PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Dahil nakadepende ang mga integral sa converged SCF orbitals, maaaring magkaiba ang isang na-regenerate na file sa orihinal na ibinigay sa orbital phase o ordering; hindi apektado ang total energies.

Pagkuha ng file. Hanapin ang FCIDump sa GitHub repository na ito. Maaari mong patakbuhin ang cell sa ibaba para kunin ito papunta sa lokasyong inaasahan ng natitirang bahagi ng tutorial.

Una, gagamitin natin ang cisolver na ibinibigay ng pyscf para makuha ang reference energy. Ito ang totoong ground state energy ng molekulang ginagamit natin. Para dito, ideklara muna natin ang norb at nelec, na siyang bilang ng mga orbitals at electrons, ayon sa pagkakasunod-sunod. Pagkatapos, ideklara natin ang h1e at h2e, na siyang one- at two-electron integrals ayon sa pagkakasunod-sunod. Gagamitin din ang lahat ng ito sa SQD sa ibang pagkakataon.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "assets/sqdrift/fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Pag-load ng Hamiltonian

Sa handa na ang kailangang data, babasahin natin ang Hamiltonian mula sa FCI file sa isang format na compatible sa qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Mga fermionic workflow gamit ang qiskit-fermions

Una nating i-map ang Hamiltonian papunta sa isang fermionic circuit model gamit ang qiskit-fermions, na nagbibigay ng transpiler passes at gates na partikular sa fermionic circuits. Gagamitin ang mga ito bago ang mga tradisyunal na transpiler passes ng Qiskit para sa workflow na ito.

Pagpapangkat ng termino (Term grouping)

Para matiyak ang reproducibility ng mga resulta, gagamitin muna natin ang canonical_order para i-sort ang mga termino batay lamang sa kanilang structure. Kaya naman nakapirmi ang order ng mga operator sa canon list. Tinitiyak nito ang reproducibility ng mga nagawang operator dahil ang QDriftTrotterization pass na gagamitin natin sa hinaharap ay nag-sample ng random indices para gumawa ng mga qDRIFT operator.

Sa hakbang na ito, sinasamantala natin ang maraming symmetry na naroroon sa electronic structure Hamiltonian sa pamamagitan ng pagpapangkat ng mga magkakaugnay na termino na may magkaparehong coefficient. Bagama't binabago nito ang distribution ng operator coefficient na pinagmumulan ng sample ng qDRIFT protocol, hindi ito nakakaapekto sa convergence guarantees nito. Mahalaga, ang pagpapangkat ng mga terminong magkaugnay dahil sa symmetry ay nagreresulta sa isang kanais-nais na pagkakansela ng mga Pauli term at sa mas maikling circuit depth kapag time-evolving ang isang state sa ilalim ng kanilang aksyon.

Ang qiskit-fermions ay nagbibigay ng group_terms_by_electronic_structure function na gumagawa ng pagpapangkat na ito para sa atin.

Tandaan na inaasahan ng group_terms_by_electronic_structure na ang mga termino ay normal ordering.

Pag-filter ng diagonal terms

Inaalis natin ang mga diagonal terms mula sa Hamiltonian na ginagamit para bumuo ng mga circuit, para magamit ang $n$ na qDRIFT sampling slots sa mga terminong naglilipat ng population sa pagitan ng mga configuration. Pinakamainam na i-filter out ang ganitong mga termino sa Hamiltonian sa puntong ito, bago buuin ang Evolution gate sa susunod na hakbang.

Ang mga terminong pinag-uusapan ay ang mga diagonal sa occupation-number basis, ibig sabihin, ang mga produkto ng number operators $a^\dagger_i a_i$. May tatlong uri ng termino na kabilang sa deskripsyong ito:

• ang constant energy offset, isang produkto ng zero number operators, na ang time evolution ay nag-aambag lamang ng global phase;

• ang individual number operators $n_i$, na ang time evolution ay nababawas sa single-qubit $Z$ rotations;

• ang higher-order products tulad ng $n_i n_j$.

Sa kanilang sarili, wala sa mga ito ang naglilipat ng population sa pagitan ng mga occupation-number configuration; umaaksyon lamang sila sa mga phase ng mga configuration na naroroon na. Hindi sila inert, gayunpaman: ang mga relative phase na iyon ay nagpapakain sa interference na nabubuo ng mga excitation term sa ibang pagkakataon sa circuit, kaya nababago ng pag-filter sa kanila ang evolution na aktwal na nabubuo at maaaring baguhin ang sampling distribution. Ito ay isang sinadyang approximation sa hakbang ng paggawa ng circuit, na ginawa upang ituon ang sampling sa mga excitation term, sa halip na isang hakbang na hindi nakaaapekto sa sampled distribution. Hindi katulad ng symmetry grouping sa itaas, na hindi humihipo sa convergence guarantees ng qDRIFT, binabago ng filter na ito ang operator na kino-convert. Kaya hindi na inaapproximate ng mga circuit ang evolution sa ilalim ng buong Hamiltonian, at ang qDRIFT error bounds ay nalalapat sa filtered operator sa halip na sa orihinal. Katanggap-tanggap ito dahil ang mga circuit ay heuristic lamang para sa sampling na ginagamit para magmungkahi ng mga configuration: walang terminong nawawala mula sa energy estimate mismo, dahil ang filter ay nalalapat lamang sa Hamiltonian na ginagamit para bumuo ng mga circuit, samantalang ginagamit ng classical diagonalization sa ibang pagkakataon ang buong Hamiltonian, kasama ang diagonal terms. Ang accuracy ng SQD ay nakadepende sa hakbang na iyon ng classical diagonalization, na nananatiling variational sa sampled subspace anuman ang paraan ng pagmumungkahi ng mga configuration.

Inaalis ng filter_diagonal_terms() function ang ganitong mga termino mula sa isang operator nang direkta (in place). Kinikilala nito ang mga ito mula sa kanilang normal-ordered structure — ang multiset ng creation modes na tumutugma sa multiset ng annihilation modes — kaya valid lamang ito sa isang operator na normal-ordered na. Hindi ito nasusuri sa runtime.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Ngayong napagpangkat na natin ang mga termino sa Hamiltonian, magdedesisyon tayo sa mga sumusunod na parameter para bumuo ng ensemble ng mga circuit:

• Ang bilang ng mga circuit na bubuuin: num_circuits
• Ang haba ng bawat circuit sa mga terminong excitation groups: num_exc
• Ang factor para sa iba't ibang evolution times: times

Paggawa ng mga fermionic circuit

Gagawa na tayo ngayon ng mga fermionic circuit para sa bawat time-step. Ang bawat circuit ay binubuo ng iisang evolution gate, kasama ang evolution time na idineklara natin kanina. Ang evolution operator ay ang Hamiltonian. Sa ibang pagkakataon, magpapatakbo tayo ng transpiler passes sa mga circuit na ito para bumuo ng mga qDRIFT circuit.

Paghahanda ng ansatz

Inihahanda natin ang Hartree-Fock state gamit ang InitializeModes class. Para sa nitrogen, ang proseso ay simpleng pag-apply ng X gates sa unang num_elec_a qubits at pagkatapos sa num_elec_b qubits, na pareho, pitong-pito, para sa nitrogen. Ang state na ito ay kumakatawan sa pito $\alpha$ at pitong $\beta$ electrons ng nitrogen.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Hakbang 2: I-optimize ang problema para sa pagsasagawa sa quantum hardware

Ngayong mayroon na tayong mga circuit, gagamitin muna natin ang mga pass na available sa qiskit-fermions para magsagawa ng fermionic level optimizations, susundan ng pag-transpile ng ating circuit para sa napiling backend. Dahil isa itong simulator experiment, gagawin muna natin ito para sa AerSimulator. Pagkalkula ng weight para sa bawat grupo

Sa hakbang na ito, isinasagawa natin ang qDRIFT sampling ng mga termino nang stochastic na may mga probabilidad na proportional sa kanilang mga coefficient sa Hamiltonian. Ginagawa ito para sa atin ng qDRIFT transpiler pass. Makakagawa na tayo ngayon ng mas mabababaw na circuit na mas episyenteng maipapatakbo sa hardware sa kabila ng limitadong qubit connectivity, kahit na naglalaman ang Hamiltonian ng long-range couplings at mga terminong mas mataas pa sa quadratic. Matapos ang term grouping, nagsa-sample ito ng mga operator batay sa kanilang mga weight. Para sa bawat operator $h_i$, ang weight $W_{h_i}$ ay nadedepine tulad ng sumusunod:

$$W_{h_i} = |c_i| / \lambda$$

Fermionic at hardware-native optimizations

Ang function na generate_preset_jw_pass_manager() ay nagbabalik ng isang MultiStagePassManager na kumukuha ng isang FermionicCircuit at gumagawa ng na-optimize na final circuit na maaari nating i-transpile para tumakbo sa ating hardware. Papalitan natin ang default nitong optimization stage ng isang FermionicPassManager na naglalaman ng ating QDriftTrotterization pass:

• Ginagamit ng QDriftTrotterization pass ang weight-calculation at sampling nang internal para bumuo ng mga circuit na gagamitin natin para sa sampling

• Ang RelabelModes pass ay isa pang optimization pass na maaaring gamitin para i-permute ang mga fermionic modes para i-optimize ang connectivity sa mga qubit at bawasan ang gate depth; magbasa pa sa API reference

Awtomatikong tatakbo ang natitirang mga stage ng MultiStagePassManager at ito ang humahawak sa buong fermion-to-qubit mapping:

• F2QLayout: Inilalapat ng preset pass manager ang TrivialF2QLayout pass, na simpleng nagma-map ng $n$ fermionic bits papunta sa $n$ qubits.

• F2QSynth: Isang transpilation pass para i-map ang mga fermion-based circuit instructions papunta sa mga qubit-based na instructions.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Ngayong tapos na tayo sa fermionic-level optimizations, maaari na nating i-transpile ang mga circuit para sa pagpapatakbo sa simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Hakbang 3: Isagawa gamit ang Qiskit primitives

Ngayong mayroon na tayong mga circuit, maaari na natin itong patakbuhin gamit ang Qiskit primitives sa AerSimulator. Pagsasamahin natin ang lahat ng counts mula sa iba't ibang circuit. Iko-convert natin ang mga ito papunta sa boolean vectors bago tuluyang mag-post-process gamit ang SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

Hakbang 4: I-post-process at ibalik ang resulta sa nais na classical format

Paggamit ng bitstrings para sa SQD

Maaari na nating patakbuhin ang diagonalization scheme sa mga napiling bitstrings para hanapin ang pinakamababang eigenvalue na tutumbas sa ground state energy ng molekula. Gagawa tayo ng callback function, magdedeklara ng initial occupancies, at itatakda ang mga parameter bago tuluyang patakbuhin ang diagonalization scheme. Ginagamit ang callback function para i-print ang kasalukuyang iteration at ang kasalukuyang eigenvalue estimate sa bawat iteration.

Sa wakas, para makuha ang ground state estimate, idadagdag natin ang nuclear_repulsion_energy sa resultang energy.

Tandaan: Hindi nakapirmi ang subspace dimension sa iba't ibang iteration, kahit sa noiseless simulator — bawat subsample ay kumukuha ng ibang set ng mga configuration, at binabago ng recovery step ang shape ng pool sa pagitan ng mga iteration, kaya nag-iiba ang naiulat na dimension mula sa isang subsample papunta sa susunod. Hindi mismo pinipirmi ng noiseless sampling ang dimension ng napiling subspace. Gayunpaman, ang hardware run ay may tendensyang magbigay ng sistematikong mas malalaking subspace, dahil sinisira ng noisy shots ang particle-number symmetry at ginagawang karagdagang basis vectors ng configuration recovery ang mga ito. Dahil dito, magpapakilala din tayo ng isa pang hakbang para sa pag-prune ng bitstrings sa hardware section.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Halimbawa sa hardware

Gumagamit ang halimbawang ito ng 20 qubits (10 spatial orbitals). Ang pagpiling iyon ay para sa kaginhawaan ng isang tutorial na dapat tumakbo nang mabilis, hindi isang mahigpit na hangganan sa method.

Hindi direktang itinatakda ng bilang ng qubit ang gastos ng classical step. Dini-diagonalize ng SQD ang Hamiltonian na na-project papunta sa subspace na sinasaklaw ng mga na-sample na configuration, kaya ang nagtutulak sa classical cost ay ang dimension ng napiling subspace na iyon — pinamamahalaan dito ng samples_per_batch, num_batches, at kung gaano karaming natatanging configuration ang aktwal na nabubuo ng mga circuit — kasama ang sparse linear algebra na kailangan para i-apply ang projected Hamiltonian. Lumalaki nang combinatorial ang buong CI space ayon sa mga orbital at electron, pero ang napiling subspace ay isang maliit, natu-tune na bahagi lamang nito, at direkta nating kinokontrol ang laki nito. Dahil dito, ang bilang ng mga qubit at ang classical difficulty ay maaaring magbago nang medyo independiyente sa isa't isa: ang mas malawak na orbital space na na-sample papunta sa isang katamtamang subspace ay maaaring mas mura kaysa sa isang mas maliit na sistema na dini-diagonalize sa napakalaking subspace.

Sa praktika, kung gayon, nakadepende ang posibleng laki ng sistema sa subspace dimension na kailangan mo para sa accuracy na gusto mo at sa memory at cores na available sa eigensolver. Karaniwang nangangailangan nga ang mas malalaking orbital space ng mas malaking subspace para makarating sa chemical accuracy, at iyan ang eventually nagbibigay-motibo para sa distributed resources — tingnan ang qiskit-addon-sqd-hpc para sa pag-scale out ng hakbang na ito. Sa halip na mag-assume ng nakapirming cutoff, ang praktikal na paraan ay panoorin ang naiulat na subspace dimension at ang energy convergence sa iba't ibang iteration at dagdagan ang laki ng subspace hanggang tumigil sa pag-improve ang energy o maubos ang available na memory.

Tandaan: Dahil sa sampling error mula sa noise sa hardware, ang subspace na nagawa para sa diagonalization sa hardware run ay magiging mas malaki kaysa sa nakukuha natin kapag gumagamit ng simulator. Habang pinapataas nito ang dimension ng subspace na gusto nating i-diagonalize, nagbibigay pa rin ang workflow ng tumpak na sagot dahil sa robustness ng SQD laban sa noise.

Pag-prune ng mga spurious strings

Dito, maaari nating piliing magsagawa ng karagdagang hakbang. Kapag nasa atin na ang lahat ng bitstrings mula sa mga pagpapatakbo ng circuit, maaari nating i-filter out ang mga invalid na bitstrings bago patakbuhin ang SQD, o magpatuloy nang walang pruning. Sa pangkalahatan, mas mainam na laktawan ang pruning para sa hardware runs, dahil iniiwan nitong available ang symmetry-broken shots para sa configuration recovery, na maaaring ayusin ang mga ito papunta sa valid na configurations at sa gayon palawakin ang subspace sa halip na direktang itapon ang mga shots na iyon.

Dahil pito lang $\alpha$ at pitong $\beta$ electrons ang maaaring meron ang nitrogen, anumang bitstrings na may higit o kulang sa pitong 1s sa una at pangalawang kalahati ng output ay maaaring itapon. Magdedepine tayo ng isang function na sumusuri kung valid ang mga bitstrings, at kung hindi, itinatapon ang mga ito. Kapag na-filter out na natin ang mga spurious na bitstrings, ang natitira ay ipinapadala papunta sa diagonalization scheme. Gamitin ang PRUNE flag sa ibaba para lumipat sa pagitan ng dalawang behavior.

Tandaan na ang pruning ay isa lamang sa ilang mga pagpipilian na humuhubog sa final subspace, kasama ang bilang ng mga circuit, ang set ng evolution times, at diagonal-term filtering. Ang paghahambing ng isang pruned run laban sa isang unpruned run ay makabuluhan lamang kung nananatiling fixed ang lahat ng iba pa; tinatalakay ito nang mas detalyado ng C++ companion, dahil nagpo-postselect ito sa halip na mag-recover at naiiba rin sa ibang mga parameter na iyon.

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Mga susunod na hakbang

Mga Rekomendasyon

Kung interesado ka sa gawaing ito, maaaring interesado ka rin sa mga sumusunod na materyal:

• Sample-based Krylov quantum diagonalization of a fermionic lattice model - isang kaugnay na tutorial na gumagamit ng time-evolution circuits sa halip na isang variational ansatz.
• Sample-based quantum diagonalization of a chemistry Hamiltonian - isang tutorial tungkol sa kung paano bumuo ng local unitary cluster Jastrow (LUCJ) circuit para sa quantum chemistry simulation.
• Ang SqDRIFT paper - ang literatura na pinagbasehan ng tutorial na ito. (Tandaan na ang ilan sa mga optimization na tinalakay sa paper na ito ay kasalukuyang isang work in progress, at maaaring magbago ang tutorial na ito sa hinaharap batay sa ebolusyon ng mga library na ginagamit.)