Algoritmul SqDRIFT pentru estimarea stării fundamentale

Estimare de utilizare: 180 de secunde pe un procesor Heron r3 (NOTĂ: aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Rezultatele învățării

• Află cum să creezi circuite cu adâncime mai mică în comparație cu trotterizarea

• Parcurge un flux de lucru complet pentru estimarea stării fundamentale folosind qDRIFT și SQD

• Află cum să folosești qiskit-fermions împreună cu alte addon-uri Qiskit pentru a implementa un astfel de flux de lucru

Acest tutorial este prezentat ca un notebook Python, în scop didactic.

Cerințe preliminare

• Citește prezentarea generală Sample-based quantum diagonalization (SQD)

• Citește lecția Sample-based Krylov Quantum Diagonalization (SKQD)

Context

SqDRIFT este o variantă a SKQD care înlocuiește necesitatea de a alege un ansatz din care să se eșantioneze șiruri de biți cu un ansamblu de circuite de evoluție în timp construite direct din Hamiltonianul țintă. Aceasta se realizează prin subeșantionarea unor operatori de evoluție în timp mai mici din Hamiltonian pe baza coeficienților săi, metodă cunoscută ca trotterizarea qDRIFT.

Acest tutorial folosește Qiskit Fermions pentru a crea circuitele fermionice mai naturale pentru algoritmul qDRIFT, urmate de folosirea pașilor de aranjare (layout) și sinteză fermionice, înainte de a integra circuitele în fluxul Qiskit tradițional pentru execuția pe hardware.

Fie Hamiltonianul de forma:

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

unde, fără pierderea generalității, cerem ca $c_i > 0$ și ca cea mai mare valoare proprie a lui $h_i$ să fie egală, în valoare absolută, cu $1$. Orice prefactor cu semn sau complex este absorbit în $h_i$, astfel încât coeficienții $c_i$ sunt ponderi strict pozitive, iar $h_i$ poartă direcția fiecărui termen. Aici $N$ este numărul de termeni (sau, după grupare, numărul de grupuri) din Hamiltonian; este o proprietate a Hamiltonianului și este distinct de numărul de operatori eșantionați într-un singur circuit, notat $n$ mai jos.

Algoritmul qDRIFT realizează apoi, pentru timpul țintă $t$, un operator $V_k$, unde $k$ merge de la $1 \cdots K$ și semnifică al $k_{th}$ circuit SqDRIFT, definit astfel:

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

Aici $n$ este numărul de operatori eșantionați pe circuit, iar $K$ este numărul de circuite din ansamblu. Produsul se întinde peste cele $n$ extrageri, nu peste toți cei $N$ termeni ai Hamiltonianului, și, deoarece termenii sunt extrași cu înlocuire, același $h_i$ poate apărea de mai multe ori într-un singur $V_k$.

Cantitatea:

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

este norma $L_1$ a coeficienților, deci fiecare dintre cei $n$ pași evoluează pe aceeași durată $\lambda t / n$, indiferent ce termen a fost extras. Uniformitatea unghiului pasului este trăsătura caracteristică a qDRIFT: un coeficient influențează rezultatul prin cât de des este extras termenul său, nu prin cât de mult este rotit acel termen. Indicii sunt eșantionați din distribuția:

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

deci seria $(k_1, \ldots, k_n)$ este o secvență aleatoare de indici de termeni extrași din această distribuție. Deoarece $c_i$ sunt pozitivi și au suma $\lambda$, aceasta este o distribuție de probabilitate normalizată, iar valoarea așteptată a canalului rezultat peste extragerile aleatoare aproximează evoluția sub $H$, cu o eroare care scade pe măsură ce $n$ crește. Observă că eroarea de aproximare depinde de $\lambda$, nu de numărul de termeni $N$.

(Lucrarea SqDRIFT notează numărul de termeni cu $\mathcal{N}$ și lungimea secvenței cu $N$; aici folosim $N$ și $n$ pentru a păstra cele două clar distincte.)

Acest tutorial arată cum se generează un ansamblu de astfel de circuite randomizate. După ce am creat aceste circuite, similar cu modul în care creăm un subspațiu Krylov pentru operatori diferiți, eșantionăm șiruri de biți din mai mulți astfel de operatori cu parametri de timp diferiți. Aceasta asigură o suprapunere mai mare între vectorii stării fundamentale și șirurile de biți eșantionate.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalat

• Un mediu virtual Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Observă că numele este la plural)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Poți instala toate pachetele necesare cu:

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

Configurare

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

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "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")

# 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,
)

Exemplu cu simulator

Pasul 1: Mapează intrările clasice la o problemă cuantică

Citirea și pregătirea FCIDump

Pentru acest tutorial, vom încărca Hamiltonianul de structură electronică pentru azot (N2). Există și alte moduri de a crea operatori fermionici. Consultă documentația la qiskit_fermions.operators.library.

Despre acest FCIDump. Fișierul N2_sto_3g descrie o moleculă de azot ($N_2$) în baza minimală STO-3G, la o distanță interatomică de 1.09 $\AA$, lungimea experimentală de echilibru a legăturii. Antetul său declară NORB=10, NELEC=14 și MS2=0: 10 orbitali spațiali (deci 20 de spin-orbitali și 20 de qubiți sub Jordan-Wigner), 14 electroni într-un singlet de spin, deci șapte electroni $\alpha$ și șapte $\beta$. Tuturor orbitalilor li se dă eticheta de simetrie 1, adică nu se exploatează nicio simetrie de grup punctual. Fiind un dump STO-3G în spațiul complet, niciun orbital nu este înghețat, iar spațiul de corelație este suficient de mic pentru ca o energie FCI exactă de referință să poată fi calculată clasic pentru comparație, așa cum se arată în celula următoare.

Un fișier echivalent poate fi regenerat cu 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")

Deoarece integralele depind de orbitalii SCF convergenți, un fișier regenerat ar putea diferi de cel livrat prin faza sau ordinea orbitalilor; energiile totale nu sunt afectate.

Obținerea fișierului. Găsești FCIDump în acest depozit GitHub. Poți rula celula de mai jos pentru a-l descărca în locația pe care o așteaptă restul tutorialului.

Mai întâi folosim cisolver furnizat de pyscf pentru a obține energia de referință. Aceasta este energia reală a stării fundamentale a moleculei cu care lucrăm. Pentru aceasta vom declara mai întâi norb și nelec, adică numărul de orbitali și, respectiv, numărul de electroni. Apoi declarăm h1e și h2e, adică integralele cu un electron și, respectiv, cu doi electroni. Toate acestea vor fi folosite ulterior și pentru SQD.

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

Încărcarea Hamiltonianului

Cu datele necesare pregătite, citim Hamiltonianul din fișierul FCI într-un format compatibil cu qiskit-fermions

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

Fluxuri de lucru fermionice cu qiskit-fermions

Vom mapa mai întâi Hamiltonianul într-un model de circuit fermionic folosind qiskit-fermions, care oferă pași de transpiler și porți specifice circuitelor fermionice. Acestea vor fi folosite ulterior înaintea pașilor tradiționali ai transpiler-ului Qiskit pentru acest flux de lucru.

Gruparea termenilor

Pentru a asigura reproductibilitatea rezultatelor, folosim mai întâi canonical_order pentru a sorta termenii doar pe baza structurii lor. Ordinea operatorilor din lista canon este deci fixă. Aceasta asigură reproductibilitatea operatorilor creați, deoarece pasul QDriftTrotterization pe care îl vom folosi ulterior eșantionează indici aleatori pentru a crea operatorii qDRIFT.

În acest pas exploatăm numeroasele simetrii prezente în Hamiltonianul de structură electronică, grupând termeni înrudiți cu coeficienți identici. Deși aceasta schimbă distribuția coeficienților operatorului din care eșantionează protocolul qDRIFT, nu afectează garanțiile sale de convergență. În mod esențial, gruparea termenilor legați prin simetrie duce la o anulare favorabilă a termenilor Pauli și la o adâncime totală mai mică a circuitului atunci când o stare este evoluată în timp sub acțiunea lor.

qiskit-fermions oferă funcția group_terms_by_electronic_structure care face această grupare pentru noi.

Observă că group_terms_by_electronic_structure presupune termeni cu ordonare normală.

Filtrarea termenilor diagonali

Eliminăm termenii diagonali din Hamiltonianul folosit pentru generarea circuitelor, astfel încât cele $n$ sloturi de eșantionare qDRIFT să fie consumate de termeni care mută populația între configurații. Cel mai bine este ca astfel de termeni să fie filtrați din Hamiltonian în acest punct, înainte ca poarta Evolution să fie construită în pasul următor.

Termenii în cauză sunt cei diagonali în baza numerelor de ocupare, adică produsele de operatori de număr $a^\dagger_i a_i$. Trei tipuri de termeni se încadrează în această descriere:

• decalajul constant de energie, un produs de zero operatori de număr, a cărui evoluție în timp contribuie doar cu o fază globală;

• operatorii de număr individuali $n_i$, a căror evoluție în timp se reduce la rotații $Z$ pe un singur qubit;

• produsele de ordin superior, cum ar fi $n_i n_j$.

Singuri, niciunul dintre aceștia nu mută populația între configurațiile de numere de ocupare; acționează doar asupra fazelor configurațiilor deja prezente. Totuși, nu sunt inerți: aceste faze relative alimentează interferența generată de termenii de excitație mai târziu în circuit, deci filtrarea lor schimbă evoluția generată efectiv și poate schimba distribuția de eșantionare. Aceasta este o aproximare deliberată în pasul de generare a circuitelor, făcută pentru a concentra eșantionarea pe termenii de excitație, și nu un pas care lasă distribuția eșantionată neatinsă. Spre deosebire de gruparea pe simetrii de mai sus, care lasă intacte garanțiile de convergență qDRIFT, acest filtru schimbă operatorul evoluat. Circuitele nu mai aproximează deci evoluția sub Hamiltonianul complet, iar limitele de eroare qDRIFT se aplică operatorului filtrat, nu celui original. Acest lucru este acceptabil aici deoarece circuitele sunt doar o euristică de eșantionare folosită pentru a propune configurații: niciun termen nu se pierde din estimarea energiei în sine, deoarece filtrul se aplică doar Hamiltonianului folosit pentru construirea circuitelor, în timp ce diagonalizarea clasică de mai târziu folosește Hamiltonianul complet, inclusiv termenii diagonali. Precizia SQD depinde de acel pas clasic, care rămâne variațional în subspațiul eșantionat, indiferent de modul în care au fost propuse configurațiile.

Funcția filter_diagonal_terms() elimină astfel de termeni dintr-un operator pe loc. Îi identifică după structura lor în ordine normală — multimulțimea modurilor de creare care coincide cu multimulțimea modurilor de anihilare — deci este validă doar pe un operator care este deja în ordine normală. Această presupunere nu este verificată la rulare.

# 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

Acum că am grupat termenii din Hamiltonian, vom decide asupra următorilor parametri pentru a genera ansamblul de circuite:

• Numărul de circuite de generat: num_circuits
• Lungimea fiecărui circuit în grupuri de excitații: num_exc
• Factorul pentru diferitele timpuri de evoluție: times

Crearea circuitelor fermionice

Vom crea acum circuite fermionice pentru fiecare dintre pașii de timp. Fiecare circuit va consta dintr-o singură poartă de evoluție, cu timpul de evoluție pe care l-am declarat mai devreme. Operatorul de evoluție este Hamiltonianul. Mai târziu rulăm pași de transpiler pe aceste circuite pentru a crea circuite qDRIFT.

Pregătirea ansatzului

Pregătim starea Hartree-Fock folosind clasa InitializeModes. Pentru azot, procesul constă pur și simplu în aplicarea porților X pe primii num_elec_a qubiți și apoi pe cei num_elec_b qubiți, ambele fiind egale cu șapte pentru azot. Această stare reprezintă cei șapte electroni $\alpha$ și șapte $\beta$ ai azotului.

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

Pasul 2: Optimizează problema pentru execuția pe hardware cuantic

Acum că avem circuitele, vom folosi mai întâi pașii disponibili în qiskit-fermions pentru a face optimizări la nivel fermionic, apoi vom transpila circuitul pentru backend-ul ales. Deoarece este un experiment cu simulator, vom face mai întâi aceasta pentru AerSimulator. Calculul ponderilor pentru fiecare grup

În acest pas, efectuăm eșantionarea qDRIFT a termenilor în mod stochastic, cu probabilități proporționale cu coeficienții lor din Hamiltonian. Pasul de transpilare qDRIFT face acest lucru pentru noi. Acum putem crea circuite mai puțin adânci, care pot fi executate mai eficient pe hardware în ciuda conectivității limitate a qubiților, chiar și atunci când Hamiltonianul conține cuplaje la distanță mare și termeni de grad mai mare decât pătratic. După gruparea termenilor, operatorii sunt eșantionați pe baza ponderilor lor. Pentru fiecare operator $h_i$, ponderea $W_{h_i}$ este definită astfel:

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

Deoarece termenii au fost grupați în Pasul 1, fiecare $h_i$ de aici este un grup întreg: $c_i$ este coeficientul absolut mediu al termenilor din grupul $i$, iar fiecare termen din grup este evoluat cu coeficientul său redus la semn.

Optimizări fermionice și native hardware

Funcția generate_preset_jw_pass_manager() returnează un MultiStagePassManager care primește un FermionicCircuit și produce un circuit final optimizat pe care îl putem transpila pentru a rula pe hardware-ul nostru. Înlocuim etapa sa implicită de optimizare cu un FermionicPassManager care conține pasul nostru QDriftTrotterization:

• Pasul QDriftTrotterization folosește intern calculul ponderilor și eșantionarea pentru a genera circuitele pe care le vom folosi la eșantionare

• Pasul RelabelModes este un alt pas de optimizare care poate fi folosit pentru a permuta modurile fermionice, pentru a optimiza conectivitatea între qubiți și a reduce adâncimea porților; citește mai multe în referința API

Etapele rămase ale MultiStagePassManager rulează automat și se ocupă de întreaga mapare fermion-qubit:

• F2QLayout: Pass manager-ul predefinit aplică pasul TrivialF2QLayout, care mapează trivial $n$ biți fermionici pe $n$ qubiți.

• F2QSynth: Un pas de transpilare care mapează instrucțiunile de circuit bazate pe fermioni în instrucțiuni bazate pe qubiți.

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

Acum că am terminat optimizările la nivel fermionic, putem transpila circuitele pentru execuția pe simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Pasul 3: Execută folosind primitivele Qiskit

Acum că avem circuitele, le putem rula folosind primitivele Qiskit pe AerSimulator. Vom combina toate numărătorile din circuitele diferite. Le convertim în vectori booleeni înainte de postprocesarea finală cu 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

Pasul 4: Postprocesează și returnează rezultatul în formatul clasic dorit

Folosirea șirurilor de biți pentru SQD

Acum putem rula schema de diagonalizare pe șirurile de biți selectate, pentru a găsi cea mai mică valoare proprie, care va corespunde energiei stării fundamentale a moleculei. Creăm o funcție callback, declarăm ocupările inițiale și setăm parametrii înainte de a rula în final schema de diagonalizare. Funcția callback este folosită pentru a afișa iterația curentă și estimarea curentă a valorii proprii la fiecare iterație.

În final, pentru a obține estimarea stării fundamentale, adăugăm nuclear_repulsion_energy la energia rezultată.

Notă: Dimensiunea subspațiului nu este fixă de-a lungul iterațiilor, nici măcar pe simulatorul fără zgomot — fiecare subeșantion extrage un set diferit de configurații, iar pasul de recuperare remodelează setul între iterații, astfel că dimensiunea raportată variază de la un subeșantion la altul. Eșantionarea fără zgomot nu fixează, prin ea însăși, dimensiunea subspațiului selectat. Rularea pe hardware tinde însă să producă subspații sistematic mai mari, deoarece măsurătorile afectate de zgomot rup simetria numărului de particule, iar recuperarea configurațiilor le transformă în vectori de bază suplimentari. Din acest motiv, vom introduce și un alt pas pentru eliminarea șirurilor de biți în secțiunea despre hardware.

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

Exemplu pe hardware

Acest exemplu folosește 20 de qubiți (10 orbitali spațiali). Această alegere este o comoditate pentru un tutorial care ar trebui să ruleze rapid, nu o limită strictă a metodei.

Costul pasului clasic nu este stabilit direct de numărul de qubiți. SQD diagonalizează Hamiltonianul proiectat pe subspațiul generat de configurațiile eșantionate, deci ceea ce determină costul clasic este dimensiunea acelui subspațiu selectat — guvernată aici de samples_per_batch, num_batches și de câte configurații distincte produc efectiv circuitele — împreună cu algebra liniară rară necesară pentru a aplica Hamiltonianul proiectat. Spațiul CI complet crește combinatorial cu numărul de orbitali și de electroni, dar subspațiul selectat este o felie mică și reglabilă din el, iar dimensiunea lui o controlăm direct. În consecință, numărul de qubiți și dificultatea clasică pot varia oarecum independent: un spațiu orbital mai larg, eșantionat într-un subspațiu modest, poate fi mai ieftin decât un sistem mai mic diagonalizat pe un subspațiu foarte mare.

În practică, așadar, dimensiunea fezabilă a sistemului depinde de dimensiunea subspațiului de care ai nevoie pentru acuratețea dorită și de memoria și nucleele disponibile pentru solverul de valori proprii. Spațiile orbitale mai mari necesită de obicei un subspațiu mai mare pentru a atinge acuratețea chimică, iar acest lucru motivează în cele din urmă resursele distribuite — vezi qiskit-addon-sqd-hpc pentru scalarea acestui pas. În loc să presupui o limită fixă, abordarea practică este să urmărești dimensiunea raportată a subspațiului și convergența energiei de-a lungul iterațiilor și să mărești dimensiunea subspațiului până când energia nu se mai îmbunătățește sau până când epuizezi memoria disponibilă.

Notă: Din cauza erorii de eșantionare provocate de zgomotul hardware-ului, subspațiul creat pentru diagonalizare în rularea pe hardware va fi mai mare decât cel obținut pe simulator. Deși crește dimensiunea subspațiului pe care vrem să-l diagonalizăm, fluxul de lucru ne oferă totuși un răspuns precis datorită robusteții SQD față de zgomot.

Eliminarea șirurilor false

Aici putem alege să efectuăm un pas suplimentar. Când avem toate șirurile de biți din execuțiile circuitelor, putem fie să filtrăm șirurile de biți invalide înainte de a rula SQD, fie să mergem mai departe fără eliminare. Renunțarea la eliminare este în general preferabilă pentru rulările pe hardware, deoarece lasă măsurătorile cu simetria ruptă disponibile pentru recuperarea configurațiilor, care le poate repara în configurații valide și astfel lărgește subspațiul, în loc să renunțe pur și simplu la acele măsurători.

Deoarece azotul poate avea doar șapte electroni $\alpha$ și șapte electroni $\beta$, orice șiruri de biți care au mai mult sau mai puțin de șapte de 1 în prima și în a doua jumătate a ieșirii pot fi eliminate. Definim o funcție care verifică dacă șirurile de biți sunt valide și, dacă nu, le elimină. După ce filtrăm șirurile false, restul sunt trimise în schema de diagonalizare. Folosește flag-ul PRUNE de mai jos pentru a comuta între cele două comportamente.

Ține cont că eliminarea este doar una dintre mai multele alegeri care modelează subspațiul final, alături de numărul de circuite, setul de timpi de evoluție și filtrarea termenilor diagonali. Compararea unei rulări cu eliminare cu una fără eliminare este informativă doar dacă tot restul este menținut fix; versiunea C++ a acestui tutorial discută acest lucru mai în detaliu, deoarece postselectează în loc să recupereze și diferă și în ceilalți parametri.

name = "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

Pașii următori

Recomandări

Dacă ți-a plăcut această lucrare, te-ar putea interesa următorul material:

• Diagonalizarea cuantică Krylov bazată pe eșantioane a unui model de rețea fermionică - un tutorial înrudit care folosește circuite de evoluție în timp în loc de un ansatz variațional.
• Diagonalizarea cuantică bazată pe eșantioane a unui Hamiltonian chimic - un tutorial despre cum se construiește un circuit local unitary cluster Jastrow (LUCJ) pentru simularea chimiei cuantice.
• Lucrarea SqDRIFT - literatura pe care se bazează acest tutorial. (Reține că unele dintre optimizările discutate în această lucrare sunt în prezent în curs de dezvoltare, iar acest tutorial poate suferi modificări în viitor, în funcție de evoluția bibliotecilor folosite.)