Sample-based quantum diagonalization ng isang chemistry Hamiltonian

Tantiya ng paggamit: mas mababa sa isang minuto sa isang Heron r2 processor (PAALALA: Ito ay tantiya lamang. Ang iyong runtime ay maaaring mag-iba.)

Mga resulta ng pag-aaral

Paano gamitin ang SQD Qiskit addon upang tantiyahin ang ground state energy ng isang molecular system gamit ang mga bitstrings na na-sample mula sa isang quantum processing unit (QPU).
Paano gamitin ang ffsim upang bumuo ng local unitary cluster Jastrow (LUCJ) circuit para sa quantum chemistry simulation.

Mga kinakailangan

Quantum chemistry at second quantization
Paggamit ng Sampler primitive upang mag-sample mula sa quantum circuits

Background

Sa tutorial na ito, ipapakita natin kung paano mag-post-process ng maingay na quantum samples upang tantiyahin ang ground state ng nitrogen molecule $\text{N}_2$ sa equilibrium bond length, sa pamamagitan ng paggamit ng SQD Qiskit addon upang ipatupad ang sample-based quantum diagonalization (SQD) algorithm. Mas maraming detalye sa software ay matatagpuan sa kaukulang dokumentasyon, kabilang ang isang simple example upang makapagsimula.

Ikinokonekta ng nitrogen chemistry ang halimbawang ito sa sustainability sa pamamagitan ng produksyon ng ammonia para sa mga fertilizer. Ang pag-unawa kung paano ina-activate ng mga catalyst ang nitrogen ay bahagi ng research tungo sa hindi gaanong energy-intensive na mga paraan ng produksyon. Ang pagkalkula ng ground state ng isang nakahiwalay na nitrogen molecule ay nagpapakilala ng mga electronic-structure method na may kaugnayan sa mas malawak na hamong iyon.

Ang tutorial na ito ay inirerekomenda para sa mga gumagamit na pamilyar sa quantum chemistry: partikular, ang pagkilala sa paghahanap ng ground state energies ng isang molekula. Para sa detalyadong walkthrough ng workflow, sumangguni sa quantum diagonalization algorithm course.

Ang SQD ay isang teknik para sa paghahanap ng eigenvalues at eigenvectors ng quantum operators, tulad ng isang quantum system Hamiltonian, sa pamamagitan ng paggamit ng quantum at distributed classical computing nang magkasama. Ang classical distributed computing ay ginagamit upang mag-proseso ng mga samples na nakuha mula sa isang quantum processor, at upang mag-project at mag-diagonalize ng isang target Hamiltonian sa isang subspace na binubuo (spanned) ng mga samples na iyon. Ang isang SQD-based workflow ay may sumusunod na hakbang:

Pumili ng circuit ansatz at ilapat ito sa isang quantum computer sa isang reference state (sa kasong ito, ang Hartree-Fock state).
Mag-sample ng mga bitstrings mula sa resultang quantum state.
Patakbuhin ang self-consistent configuration recovery procedure sa mga bitstrings upang makuha ang approximation sa ground state.

Ang SQD ay kilala na gumagana nang maayos kapag ang target eigenstate ay sparse: ang suporta (support) ng wave function ay nasa isang set ng basis states $\mathcal{S} = \{|x\rangle \}$ na ang laki ay hindi tumataas nang exponentially sa laki ng problema.

Quantum chemistry

Ang Hamiltonian ng isang molecular system ay maaaring isulat bilang

\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},

kung saan ang $h_{pr}$ at $h_{prqs}$ ay mga complex numbers na tinatawag na molecular integrals na maaaring kalkulahin mula sa specification ng molekula gamit ang isang computer program. Sa tutorial na ito, kinakalkula natin ang mga integrals gamit ang PySCF software package.

Para sa mga detalye tungkol sa kung paano na-derive ang molecular Hamiltonian, kumonsulta sa isang textbook sa quantum chemistry (halimbawa, Modern Quantum Chemistry ni Szabo at Ostlund). Para sa isang high-level explanation kung paano ang quantum chemistry problems ay nima-map sa quantum computers, tingnan ang lecture na Mapping Problems to Qubits mula sa Qiskit Global Summer School 2024.

Local unitary cluster Jastrow (LUCJ) ansatz

Ang SQD ay nangangailangan ng isang quantum circuit ansatz upang kumuha ng mga samples. Sa tutorial na ito, gagamitin natin ang local unitary cluster Jastrow (LUCJ) ansatz dahil sa kombinasyon nito ng physical motivation at hardware-friendliness. Gagamitin natin ang ffsim upang itayo ang ansatz circuit.

Ang LUCJ ansatz ay nag-aadapt sa mga QPU na may limitadong qubit connectivity. Ang mga spin orbital ay nima-map sa mga qubit upang ang ansatz ay hindi nangangailangan ng routing gamit ang mga SWAP gate. Ang IBM® hardware ay may square lattice o heavy-hex lattice qubit topology. Para sa heavy-hex lattice, maaari tayong gumamit ng "zig-zag" pattern, na inilalarawan sa ibaba. Sa pattern na ito, ang mga orbitals na may parehong spin ay nima-map sa mga qubits na may line topology (pula at asul na bilog), at ang koneksyon sa pagitan ng mga orbitals na may ibang spin ay naroroon sa bawat ika-4 na spatial orbital, na ang koneksyon ay pinadadali ng isang ancilla qubit (lila na bilog).

Qubit mapping diagram for the LUCJ ansatz on a heavy-hex lattice

Self-consistent na pagbawi ng configuration

Ang self-consistent configuration recovery procedure ay dinisenyo upang kunin ang pinakamaraming signal hangga't maaari mula sa maingay na quantum samples. Dahil ang molecular Hamiltonian ay nag-conserve ng particle number at spin Z, may katuturan na pumili ng circuit ansatz na nag-conserve din ng mga simetriya na ito. Kapag inilapat sa Hartree-Fock state, ang resultang state ay may fixed particle number at spin Z sa noiseless setting. Samakatuwid, ang spin- $\alpha$ at spin- $\beta$ halves ng anumang bitstring na na-sample mula sa state na ito ay dapat magkaroon ng parehong Hamming weight tulad ng sa Hartree-Fock state. Dahil sa presensya ng noise sa kasalukuyang quantum processors, ang ilang nasukat na bitstrings ay lalabag sa property na ito. Ang isang simpleng anyo ng postselection ay magtatapon ng mga bitstrings na ito, ngunit sayang ito dahil ang mga bitstrings na iyon ay maaaring maglaman pa ng ilang signal. Ang self-consistent recovery procedure ay sinusubukang mabawi ang ilan sa signal na iyon sa post-processing. Ang procedure ay iterative at nangangailangan bilang input ng isang tantiya ng average occupancies ng bawat orbital sa ground state, na unang kinakalkula mula sa raw samples. Ang procedure ay pinapatakbo sa isang loop, at bawat iteration ay may sumusunod na mga hakbang:

Para sa bawat bitstring na lumalabag sa tinukoy na mga simetriya, i-flip ang mga bits nito gamit ang isang probabilistic procedure na dinisenyo upang dalhin ang bitstring nang mas malapit sa kasalukuyang tantiya ng average orbital occupancies, upang makakuha ng bagong bitstring.
Kolektahin ang lahat ng lumang at bagong bitstrings na nakakatugon sa mga simetriya, at mag-subsample ng mga subsets ng fixed size, na pinili nang maaga.
Para sa bawat subset ng bitstrings, i-project ang Hamiltonian sa subspace na spanned ng kaukulang basis vectors (tingnan ang nakaraang seksyon para sa paglalarawan ng mga basis vectors na ito), at kalkulahin ang ground state estimate ng projected Hamiltonian sa isang classical computer.
I-update ang tantiya ng average orbital occupancies gamit ang ground state estimate na may pinakamababang energy.

Workflow diagram ng SQD

Ang SQD workflow ay inilalarawan sa sumusunod na diagram:

Workflow diagram of the SQD algorithm

Requirements

Bago simulan ang tutorial na ito, tiyaking mayroon kang sumusunod na naka-install:

Qiskit SDK v1.0 o mas bago, na may visualization support
Qiskit Runtime v0.22 o mas bago (pip install qiskit-ibm-runtime)
SQD Qiskit addon v0.12 o mas bago (pip install qiskit-addon-sqd)
ffsim v0.0.75 o mas bago (pip install ffsim)

Setup

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"ffsim": "ffsim", "matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Maliit na halimbawa ng simulator

Sa tutorial na ito, hahanapin natin ang approximation sa ground state ng nitrogen molecule malapit sa equilibrium bond distance nito. Gagamitin muna natin ang isang maliit na STO-6G basis set upang masimulate natin ang eksperimento at masiguro na gumagana ito.

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

Una, tutukuyin natin ang molekula at ang mga katangian nito.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="sto-6g",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

converged SCF energy = -108.464957764796
CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000
norb = 8
nelec = (5, 5)

Bago itayo ang LUCJ ansatz circuit, magsasagawa muna tayo ng CCSD calculation sa sumusunod na code cell. Ang $t_1$ at $t_2$ amplitudes mula sa calculation na ito ay gagamitin upang i-initialize ang mga parameters ng 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]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052352

Ngayon, gagamitin natin ang ffsim upang likhain ang ansatz circuit. Dahil ang ating molekula ay may closed-shell Hartree-Fock state, gagamitin natin ang spin-balanced variant ng UCJ ansatz, UCJOpSpinBalanced. Itatakda natin ang optimize=True sa from_t_amplitudes method upang paganahin ang "compressed" double-factorization ng $t_2$ amplitudes (tingnan ang The local unitary cluster Jastrow (LUCJ) ansatz sa documentation ng ffsim para sa mga detalye).

Dahil ang LUCJ ansatz ay nag-aadapt sa available connectivity ng QPU, kailangan nating i-initialize ang QPU backend bago likhain ang ansatz. Sa ngayon, lilikha tayo ng generic backend na may square lattice coupling map at gate set kung saan natural na nade-decompose ang LUCJ ansatz. Pagkatapos, gagamitin natin ang ffsim.qiskit.generate_lucj_pass_manager upang lumikha ng pass manager na iniangkop para sa pag-transpile ng LUCJ ansatz sa ibinigay na backend ayon sa "ladder" layout (tingnan ang paliwanag sa dokumentasyon ng ffsim). Ang function na ito ay gumagamit ng scoring heuristic upang mabawasan ang mga error na nauugnay sa napiling layout, na mahalaga kung ang iyong backend ay isang tunay na QPU, o isang simulator na may noise model. Bukod sa pagbabalik ng pass manager, nagbabalik din ang function na ito ng alpha-beta coupling pairs na maaaring ipatupad sa hardware. Kung hindi lahat ng pairs ay maaaring ipatupad, maglalabas ito ng babala.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 4
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
coupling_map = CouplingMap.from_grid(norb, 2)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

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

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

Hakbang 2: I-optimize para sa quantum hardware execution

Susunod, io-optimize natin ang circuit para sa isang target hardware. Karaniwang kasama sa hakbang na ito ang pag-initialize ng hardware backend at isang pass manager para sa backend na iyon. Gayunpaman, dahil ang LUCJ ansatz ay nag-aadapt sa hardware connectivity, ginawa na natin ang mga aksyon na ito sa nakaraang hakbang. Ang natitira na lang ay patakbuhin ang pass manager sa circuit upang ma-transpile ito sa isang ISA circuit na maaaring direktang i-execute sa QPU.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})

Hakbang 3: I-execute gamit ang Qiskit primitives

Pagkatapos i-optimize ang circuit para sa hardware execution, handa na tayong patakbuhin ito sa target hardware at mangolekta ng samples para sa ground state energy estimation. Dahil mayroon lamang tayong isang circuit, gagamitin natin ang IBM Quantum Compute Service Job execution mode at i-execute ang ating circuit.

rng = np.random.default_rng(12345)
sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

primitive_result = job.result()
pub_result = primitive_result[0]

Hakbang 4: Mag-post-process at ibalik ang resulta sa gustong classical format

Ang isang kapaki-pakinabang na sukatan upang husgahan ang kalidad ng output ng QPU ay ang bilang ng mga valid na configuration na ibinalik. Ang isang valid na configuration ay may tamang particle number at spin Z, na nangangahulugang ang kanang kalahati ng bitstring ay may Hamming weight na katumbas ng bilang ng spin-up electrons, at ang kaliwang kalahati ay may Hamming weight na katumbas ng bilang ng spin-down electrons. Ang sumusunod na cell ay kinakalkula ang fraction ng mga na-sample na configuration na valid.

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
    )

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")

Fraction of sampled configurations that are valid: 1.0

Lahat ng bitstring ay valid dahil sinesample natin ang circuit sa isang noiseless simulator. Kapag nagpapatakbo sa isang maingay na QPU, ang fraction ay magiging mas mababa sa isa, pero sana mas malaki ito kaysa sa fraction na inaasahan mo kung ang mga bitstring ay sinesample nang pantay-pantay at random, na kinukwenta sa sumusunod na cell.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Ngayon, tatantiyahin natin ang ground state energy ng Hamiltonian gamit ang diagonalize_fermionic_hamiltonian function. Ginagawa ng function na ito ang self-consistent configuration recovery procedure upang iteratively pag-refine ang maingay na quantum samples upang mapabuti ang energy estimate. Ipapasa natin ang isang callback function upang mai-save natin ang intermediate results para sa paglaon na pagsusuri. Tingnan ang API documentation para sa mga paliwanag ng mga arguments sa diagonalize_fermionic_hamiltonian.

Dito, ginagamit natin ang initial_occupancies argument ng diagonalize_fermionic_hamiltonian upang itakda ang Hartree-Fock configuration bilang paunang hula para sa orbital occupancies sa ground state. Makatuwiran ang paraang ito para sa mga sistemang may malaking suporta ang ground state sa Hartree-Fock configuration, ngunit maaaring hindi ito angkop sa ibang sitwasyon, bagaman ang mas advanced na computational methods ay maaaring magbigay ng mas mahusay na paunang hula sa mga kasong iyon. Ang pagtukoy ng initial_occupancies ay nagpapahintulot din sa configuration recovery na tumakbo kahit walang valid na configurations na na-sample, gaya ng maaaring mangyari kapag nagsa-sample ng malaking circuit sa isang maingay na QPU. Kung wala ang argument na ito, mabibigo ang configuration recovery at magtataas ng error kapag walang naibigay na valid na configuration.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 50
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

Iteration 1
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Iteration 2
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Final energy: -108.59589540253083
Final energy error: 9.19484552071026e-05

I-visualize ang mga resulta

Ang unang plot ay nagpapakita na sa simulation na ito ay nasa loob na tayo ng 1 mH ng eksaktong sagot pagkatapos ng unang iteration (ang chemical accuracy ay karaniwang tinatanggap na 1 kcal/mol $\approx$ 1.6 mH). Gayunman, maliit na sistema lang ito, at dahil walang ingay ang mga sample, hindi na kailangan ang configuration recovery. Sa isang mas malaking sistema na pinatakbo sa isang maingay na QPU, maaaring kailanganin ang maraming configuration recovery iterations, at ang panghuling katumpakan ay maaaring mas masama. Sa pangkalahatan, ang energy ay maaaring mapabuti sa pamamagitan ng pagpapahintulot ng mas maraming configuration recovery iterations o sa pamamagitan ng pagtaas ng bilang ng sample sa bawat batch.

Ang pangalawang plot ay nagpapakita ng average occupancy ng bawat spatial orbital pagkatapos ng huling iteration. Makikita natin na parehong ang spin-up at spin-down electrons ay nag-occupy ng unang limang orbitals na may mataas na probability sa ating mga solusyon.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-5, 1e-2)
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})

plt.tight_layout()
plt.show()

Output of the previous code cell

Halimbawa ng malaking sukat sa tunay na hardware

Ngayon, magpapatakbo tayo ng mas malaking halimbawa sa tunay na quantum hardware. Dito, kukuha tayo ng active space para sa nitrogen molecule mula sa cc-pVDZ basis set. Gagawin din natin ang mga sumusunod na pagbabago:

Babaan ang n_reps sa 1 para mabawasan ang noise mula sa gate errors
Itatakda ang ansatz connectivity batay sa backend processor family (heavy-hex para sa Heron, square para sa Nighthawk)
Dagdagan ang max_dim sa 1000 para payagan ang mas malalaking subspace dimension.

Mga Hakbang 1-4

Dito pinagsasama natin ang lahat ng hakbang sa isang solong workflow sa mas malaking sukat, na pagkatapos ay pinapatakbo sa tunay na quantum hardware.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="cc-pvdz",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(f"Using backend {backend.name}")

# Create pass manager
connectivity_map = {"Heron": "heavy-hex", "Nighthawk": "square"}
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity_map[
        backend.configuration().processor_type["family"]
    ],
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

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

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

# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 1000
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the
# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-3, 1e-1)
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})

plt.tight_layout()
plt.show()

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450041
Using backend ibm_boston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.01677
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.17922180754668
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.17713541195594
		Subspace dimension: 962361
	Subsample 2
		Energy: -109.17407767345108
		Subspace dimension: 968256
Iteration 2
	Subsample 0
		Energy: -109.1873108379026
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.1876970140745
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.18255107152576
		Subspace dimension: 1000000
Iteration 3
	Subsample 0
		Energy: -109.19052148187927
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19316375313544
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19087348337145
		Subspace dimension: 1000000
Iteration 4
	Subsample 0
		Energy: -109.19374847699348
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19517039328333
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19422416032815
		Subspace dimension: 1000000
Iteration 5
	Subsample 0
		Energy: -109.19580768221108
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19608553581007
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19548131833729
		Subspace dimension: 1000000
Final energy: -109.19608553581007
Final energy error: 0.03194368084709254

Output of the previous code cell

Mga susunod na hakbang

Recommendations

Kung nakita mong kawili-wili ang gawaing ito, maaaring interesado ka sa 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 variational ansatz
Scale SQD chemistry workflows with Dice solver - isang pahina na nagpapakita kung paano gamitin ang mas episyenteng Dice software para sa diagonalization
SQD addon API documentation - sanggunian para sa diagonalize_fermionic_hamiltonian function
Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer - ang papel na batayan ng tutorial na ito

Mga resulta ng pag-aaral​

Mga kinakailangan​

Background​

Quantum chemistry​

Local unitary cluster Jastrow (LUCJ) ansatz​

Self-consistent na pagbawi ng configuration​

Workflow diagram ng SQD​

Requirements​

Setup​

Maliit na halimbawa ng simulator​

Hakbang 1: I-map ang mga classical inputs sa isang quantum problem​

Hakbang 2: I-optimize para sa quantum hardware execution​

Hakbang 3: I-execute gamit ang Qiskit primitives​

Hakbang 4: Mag-post-process at ibalik ang resulta sa gustong classical format​

I-visualize ang mga resulta​

Halimbawa ng malaking sukat sa tunay na hardware​

Mga Hakbang 1-4​

Mga susunod na hakbang​