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

Cauți versiunea C++?

Acest tutorial folosește Python. Pentru implementarea în C++, inclusiv codul sursă și instrucțiunile de compilare, consultă tutorialul C++ SqDRIFT.

Rezultate ale învățării

• Învață 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

• Învață 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 scopuri didactice.

Cerințe preliminare

• Citește prezentarea generală a Diagonalizării cuantice bazate pe eșantioane (SQD)

• Citește lecția Diagonalizare cuantică Krylov bazată pe eșantioane (SKQD)

Context

SqDRIFT este o variantă a SKQD care elimină nevoia de a alege un ansatz din care să eșantionezi șiruri de biți, folosind în schimb un ansamblu de circuite de evoluție temporală construite direct din Hamiltonianul țintă. Acest lucru se realizează prin subeșantionarea unor operatori de evoluție temporală mai mici din Hamiltonian, pe baza coeficienților săi, cunoscută sub numele de metoda de Trotterizare qDRIFT.

Acest tutorial folosește Qiskit Fermions pentru a crea circuite fermionice mai naturale pentru algoritmul qDRIFT, urmate de utilizarea pass-urilor de layout și sinteză fermionică înainte de a introduce circuitele în pipeline-ul tradițional Qiskit pentru execuția pe hardware.

Fie Hamiltonianul de forma:

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

unde, fără a pierde din generalitate, 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$, deci coeficienții $c_i$ sunt ponderi strict pozitive, în timp ce $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 anumit operator $V_k$, unde $k$ merge de la $1 \cdots K$ și semnifică al $k_{th}$ circuit SqDRIFT, definit ca:

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

Aici $n$ este numărul de operatori eșantionați per circuit, iar $K$ este numărul de circuite din ansamblu. Produsul se face peste cele $n$ trageri, nu peste toți cei $N$ termeni ai Hamiltonianului, iar deoarece termenii sunt trași cu revenire, 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 de termenul care a fost tras. Uniformitatea unghiului de pas este caracteristica esențială a qDRIFT: un coeficient influențează rezultatul prin cât de des este tras 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 trași din această distribuție. Deoarece $c_i$ sunt pozitivi și însumează $\lambda$, aceasta este o distribuție de probabilitate normalizată, iar valoarea așteptată a canalului rezultat, peste tragerile 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$.

(Articolul SqDRIFT scrie numărul de termeni ca $\mathcal{N}$ și lungimea secvenței ca $N$; noi folosim $N$ și $n$ aici pentru a-i menține clar distincți.)

Acest tutorial arată cum să generezi un ansamblu de astfel de circuite randomizate. După ce am creat aceste circuite, similar modului î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 (Reține 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 — 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,
)

Exemplu de 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 separare interatomică de 1,09 $\AA$, lungimea de echilibru a legăturii experimentale. Antetul său declară NORB=10, NELEC=14 și MS2=0: 10 orbitali spațiali (deci 20 de orbitali de spin și 20 de qubiți sub transformarea Jordan-Wigner), 14 electroni într-un singlet de spin, deci șapte electroni $\alpha$ și șapte $\beta$. Toți orbitalii primesc eticheta de simetrie 1, adică nu se exploatează nicio simetrie de grup punctual. Fiind un dump STO-3G pe întregul spațiu, niciun orbital nu este înghețat, iar spațiul de corelare este suficient de mic încât o energie de referință FCI exactă poate 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 poate diferi de cel livrat în privința fazei sau ordinii orbitalilor; energiile totale nu sunt afectate.

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

Mai întâi folosim cisolver oferit de pyscf pentru a obține energia de referință. Aceasta este energia reală a stării fundamentale a moleculei cu care lucrăm. Pentru asta, vom declara mai întâi norb și nelec, care reprezintă numărul de orbitali, respectiv numărul de electroni. Apoi declarăm h1e și h2e, care reprezintă integralele cu unul, respectiv două corpuri electronice. Toate acestea vor fi folosite mai târziu ș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 = "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

Î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. Aceștia vor fi folosiți mai târziu, înainte de pașii tradiționali de transpiler ai 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 astfel fixă. Acest lucru asigură reproductibilitatea operatorilor creați deoarece pasul QDriftTrotterization, pe care îl vom folosi mai târziu, eșantionează indici aleatorii pentru a crea operatorii qDRIFT.

În acest pas, exploatăm multiplele simetrii prezente în hamiltonianul structurii electronice, grupând termenii înrudiți cu coeficienți identici. Deși acest lucru schimbă distribuția coeficienților operatorului din care eșantionează protocolul qDRIFT, acest fapt nu afectează garanțiile sale de convergență. Esențial, gruparea termenilor înrudiți prin simetrie rezultă într-o anulare favorabilă a termenilor Pauli și într-o adâncime globală 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 realizează această grupare pentru noi.

Reține că group_terms_by_electronic_structure presupune termeni în ordine normală.

Filtrarea termenilor diagonali

Eliminăm termenii diagonali din hamiltonianul folosit pentru a genera circuitele, astfel încât cele $n$ sloturi de eșantionare qDRIFT să fie folosite pentru termenii care mută populația între configurații. Astfel de termeni sunt cel mai bine filtrați din hamiltonian în acest punct, înainte ca poarta Evolution să fie construită în pasul următor.

Termenii în discuție sunt cei diagonali în baza numărului de ocupare, adică produsele operatorilor 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 individuali de număr $n_i$, a căror evoluție în timp se reduce la rotații $Z$ pe un singur qubit;

• produsele de ordin superior, precum $n_i n_j$.

Luate separat, niciunul dintre acestea nu mută populația între configurațiile de număr de ocupare; ele acționează doar asupra fazelor configurațiilor deja prezente. Totuși, nu sunt inerte: aceste faze relative alimentează interferența generată de termenii de excitație mai târziu în circuit, astfel încât 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 circuitului, făcută pentru a concentra eșantionarea pe termenii de excitație, și nu un pas care lasă neatinsă distribuția eșantionată. Spre deosebire de gruparea prin simetrie de mai sus, care lasă intacte garanțiile de convergență ale qDRIFT, acest filtru schimbă operatorul care este evoluat. Prin urmare, circuitele nu mai aproximează evoluția sub hamiltonianul complet, iar limitele de eroare ale 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 a construi circuitele, în timp ce diagonalizarea clasică folosește ulterior 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 direct (in place). Ea îi identifică pe baza structurii lor în ordine normală — multisetul modurilor de creare care se potrivește cu multisetul modurilor de anihilare — astfel încât este validă doar pe un operator care este deja în ordine normală. Această presupunere nu este verificată la momentul execuției.

# 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 stabili următorii parametri pentru a genera ansamblul de circuite:

• Numărul de circuite de generat: num_circuits
• Lungimea fiecărui circuit în termeni de grupuri de excitație: num_exc
• Factorul pentru diferitele timpi 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 declarat anterior. 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 num_elec_b qubiți, ambele fiind egale cu șapte pentru azot. Această stare reprezintă cei șapte electroni $\alpha$ și cei ș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: Optimizarea problemei pentru execuția pe hardware cuantic

Acum că avem circuitele noastre, vom folosi mai întâi pașii disponibili în qiskit-fermions pentru a efectua optimizări la nivel fermionic, urmate de transpilarea circuitului nostru pentru backend-ul ales. Deoarece acesta este un experiment cu simulator, vom face mai întâi acest lucru pentru AerSimulator. Calculul ponderii pentru fiecare grup

În acest pas, efectuăm eșantionarea qDRIFT a termenilor stocastic, cu probabilități proporționale cu coeficienții lor din hamiltonian. Pasul de transpiler qDRIFT face acest lucru pentru noi. Putem crea acum circuite mai puțin adânci, care pot fi executate pe hardware mai eficient în ciuda conectivității limitate a qubiților, chiar și atunci când hamiltonianul conține cuplaje pe rază lungă și termeni de ordin mai mare decât pătratic. După gruparea termenilor, se eșantionează operatorii pe baza ponderilor lor. Pentru fiecare operator $h_i$, ponderea $W_{h_i}$ este definită după cum urmează:

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

Optimizări fermionice și native hardware-ului

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 pentru 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 maparea completă fermion-qubit:

• F2QLayout: Managerul de pași presetat aplică pasul TrivialF2QLayout, care mapează trivial $n$ biți fermionici la $n$ qubiți.

• F2QSynth: Un pas de transpilare pentru a mapa instrucțiunile de circuit bazate pe fermioni la cele 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 cu optimizările la nivel fermionic, putem transpila circuitele pentru execuție pe simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Pasul 3: Execuția folosind primitivele Qiskit

Acum că avem circuitele noastre, le putem rula folosind primitivele Qiskit pe AerSimulator. Vom combina toate numărătorile (counts) din diferite circuite. Le convertim în vectori booleeni înainte de a le post-procesa în 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: Post-procesarea și returnarea rezultatului în formatul clasic dorit

Folosirea șirurilor de biți pentru SQD

Putem rula acum 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 de callback, declarăm ocupările inițiale și setăm parametrii înainte de a rula, în final, schema de diagonalizare. Funcția de 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ă între iterații, chiar și pe simulatorul fără zgomot — fiecare subeșantion extrage un set diferit de configurații, iar pasul de recuperare remodelează grupul între iterații, astfel încât 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, însă, tinde să dea sistematic subspații mai mari, deoarece shot-urile zgomotoase 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 (pruning) ș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 trebuie 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, astfel încât 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 de fapt 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 electroni, dar subspațiul selectat este o felie mică, reglabilă din acesta, iar noi îi controlăm direct dimensiunea. Prin urmare, 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 unul foarte mare.

Prin urmare, în practică, dimensiunea fezabilă a sistemului depinde de dimensiunea subspațiului de care ai nevoie pentru precizia dorită și de memoria și nucleele disponibile pentru rezolvatorul de valori proprii. Spațiile orbitale mai mari necesită de obicei un subspațiu mai mare pentru a atinge precizia chimică, iar acest lucru este ceea ce motivează în cele din urmă resursele distribuite — vezi qiskit-addon-sqd-hpc pentru extinderea acestui pas. În loc să presupui un prag fix, abordarea practică este să urmărești dimensiunea raportată a subspațiului și convergența energiei de-a lungul iterațiilor și să crești dimensiunea subspațiului până când energia încetează să se îmbunătățească sau epuizezi memoria disponibilă.

Notă: Din cauza erorii de eșantionare provocate de zgomotul din hardware, subspațiul creat pentru diagonalizare în rularea pe hardware va fi mai mare decât cel obținut atunci când folosim simulatorul. Deși acest lucru 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 la 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 (pruning). Omiterea eliminării este în general preferabilă pentru rulările pe hardware, deoarece lasă shot-urile cu simetrie ruptă disponibile pentru recuperarea configurațiilor, care le poate repara în configurații valide și astfel poate lărgi subspațiul în loc să elimine acele shot-uri direct.

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

Ține minte 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ă restul rămâne fix; materialul însoțitor C++ discută acest lucru mai în detaliu, deoarece acesta postselectează în loc să recupereze și de asemenea diferă în privința celorlalți parametri.

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

Pașii următori

Recomandări

Dacă ai găsit această lucrare interesantă, poate te interesează următoarele materiale:

• 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.
• Articolul SqDRIFT - literatura pe care se bazează acest tutorial. (Reține că unele dintre optimizările discutate în această lucrare sunt în prezent în lucru, iar acest tutorial este supus schimbării în viitor pe baza evoluției bibliotecilor folosite.)