אלגוריתם SqDRIFT להערכת מצב היסוד

הערכת שימוש: 180 שניות על מעבד Heron r3 (הערה: זו הערכה בלבד. זמן הריצה שלכם עשוי להיות שונה.)

מטרות הלמידה

למדו איך ליצור מעגלים בעומק קטן יותר בהשוואה ל-Trotterization
עברו על זרימת עבודה מקצה לקצה להערכת מצב היסוד באמצעות qDRIFT ו-SQD
למדו איך להשתמש ב-qiskit-fermions יחד עם תוספי Qiskit אחרים כדי לממש זרימת עבודה כזו

המדריך הזה מוצג כמחברת Python למטרות הוראה.

דרישות קדם

קראו את הסקירה של דיאגונליזציה קוונטית מבוססת דגימות (SQD)
קראו את השיעור דיאגונליזציה קוונטית של Krylov מבוססת דגימות (SKQD)

רקע

SqDRIFT היא גרסה של SKQD שמחליפה את הצורך לבחור ansatz שממנו דוגמים מחרוזות ביטים באנסמבל של מעגלי התפתחות בזמן שנבנים ישירות מההמילטוניאן של היעד. זה מושג על ידי דגימת משנה של אופרטורי התפתחות בזמן קטנים יותר מההמילטוניאן לפי המקדמים שלו, והשיטה הזו ידועה כשיטת הטרוטריזציה qDRIFT.

המדריך הזה משתמש ב-Qiskit Fermions כדי ליצור את המעגלים הפרמיוניים הטבעיים יותר עבור האלגוריתם qDRIFT, ואחר כך משתמש במעברי פריסה וסינתזה פרמיוניים לפני חיבור המעגלים לצנרת Qiskit המסורתית להרצה בחומרה.

יהי ההמילטוניאן מהצורה:

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

כאשר, בלי הגבלת הכלליות, אנו דורשים $c_i > 0$ ושהערך העצמי הגדול ביותר של $h_i$ יהיה שווה, בערכו המוחלט, ל- $1$ . כל גורם כפלי מסומן או מרוכב נבלע ב- $h_i$ , כך שהמקדמים $c_i$ הם משקלים חיוביים בהחלט בעוד ש- $h_i$ נושאים את הכיוון של כל איבר. כאן $N$ הוא מספר האיברים (או, לאחר קיבוץ, מספר הקבוצות) בהמילטוניאן; זו תכונה של ההמילטוניאן והיא שונה ממספר האופרטורים הנדגמים למעגל בודד, שנכתב להלן $n$ .

אלגוריתם qDRIFT ממש אז, עבור זמן היעד $t$ , אופרטור כלשהו $V_k$ , כאשר $k$ רץ מ- $1 \cdots K$ ומציין את מעגל SqDRIFT ה- $k_{th}$ , המוגדר כך:

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

כאן $n$ הוא מספר האופרטורים הנדגמים לכל מעגל ו- $K$ הוא מספר המעגלים באנסמבל. המכפלה רצה על $n$ ההגרלות ולא על כל $N$ איברי ההמילטוניאן, ומכיוון שהאיברים נדגמים עם החזרה, אותו $h_i$ יכול להופיע יותר מפעם אחת ב- $V_k$ בודד.

הגודל:

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

הוא הנורמה $L_1$ של המקדמים, ולכן כל אחד מ- $n$ הצעדים מתפתח לאותה משך $\lambda t / n$ בלי קשר לאיזה איבר הוגרל. אחידות זווית הצעד היא המאפיין של qDRIFT: מקדם משפיע על התוצאה דרך כמה פעמים האיבר שלו מוגרל, ולא דרך כמה האיבר הזה מסתובב. האינדקסים נדגמים מההתפלגות:

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

ולכן הסדרה $(k_1, \ldots, k_n)$ היא רצף אקראי של אינדקסי איברים שהוגרל מההתפלגות הזו. מכיוון ש- $c_i$ חיוביים ומסתכמים ל- $\lambda$ , זו התפלגות הסתברות מנורמלת, והתוחלת של הערוץ המתקבל על פני ההגרלות האקראיות מקרבת התפתחות תחת $H$ , עם שגיאה שקטנה ככל ש- $n$ גדל. שימו לב ששגיאת הקירוב תלויה ב- $\lambda$ ולא במספר האיברים $N$ .

(מאמר SqDRIFT כותב את מספר האיברים כ- $\mathcal{N}$ ואת אורך הרצף כ- $N$ ; כאן אנו משתמשים ב- $N$ וב- $n$ כדי לשמור על הבחנה ברורה בין השניים.)

המדריך הזה מראה איך ליצור אנסמבל של מעגלים אקראיים כאלה. אחרי שיצרנו את המעגלים האלה, בדומה לאופן שבו אנו יוצרים תת-מרחב Krylov לאופרטורים שונים, אנו דוגמים מחרוזות ביטים ממספר אופרטורים כאלה עם פרמטרי זמן שונים. זה מבטיח חפיפה גבוהה יותר בין וקטורי מצב היסוד לבין מחרוזות הביטים שנדגמו.

דרישות

לפני שמתחילים את המדריך הזה, ודאו שהתקנתם

סביבה וירטואלית של Python (>=3.10)
pip>=25.1
qiskit ~= 2.5
qiskit-fermions==0.1.0 (שימו לב שהשם הוא ברבים)
numpy
pyscf
qiskit-aer
qiskit-ibm-runtime
qiskit-addon-sqd

אפשר להתקין את כל החבילות הנדרשות באמצעות:

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

הגדרה

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

דוגמה עם סימולטור

שלב 1: מיפוי קלטים קלאסיים לבעיה קוונטית

קריאה והכנה של ה-FCIDump

במדריך הזה נטען את ההמילטוניאן של המבנה האלקטרוני של חנקן (N2). יש גם דרכים אחרות ליצור אופרטורים פרמיוניים. עיינו בתיעוד ב-qiskit_fermions.operators.library.

על ה-FCIDump הזה. הקובץ N2_sto_3g מתאר מולקולת חנקן ( $N_2$ ) בבסיס STO-3G מינימלי במרחק בין-אטומי של 1.09 $\AA$ , אורך הקשר בשיווי משקל הניסויי. הכותרת שלו מצהירה NORB=10, NELEC=14 ו-MS2=0: 10 אורביטלים מרחביים (ולכן 20 אורביטלי ספין, ו-20 קיוביטים תחת Jordan-Wigner), 14 אלקטרונים בסינגלט ספין, כלומר שבעה אלקטרוני $\alpha$ ושבעה אלקטרוני $\beta$ . לכל האורביטלים ניתנה תווית סימטריה 1, כלומר לא מנוצלת סימטריית חבורת נקודה. מכיוון שזה dump מלא של STO-3G, אין אורביטלים קפואים ומרחב המתאם קטן מספיק כדי לחשב קלאסית אנרגיית ייחוס FCI מדויקת להשוואה, כפי שמוצג בתא הבא.

אפשר ליצור מחדש קובץ מקביל באמצעות 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")

מכיוון שהאינטגרלים תלויים באורביטלי SCF שהתכנסו, קובץ שנוצר מחדש עשוי להיות שונה מזה שסופק בפאזה או בסדר האורביטלים; האנרגיות הכוללות לא מושפעות.

השגת הקובץ. מצאו את ה-FCIDump ב-מאגר GitHub הזה. אפשר להריץ את התא שלהלן כדי להוריד אותו למיקום שבו שאר המדריך מצפה לו.

קודם נשתמש ב-cisolver שמספק pyscf כדי לקבל את אנרגיית הייחוס. זו אנרגיית מצב היסוד האמיתית של המולקולה שאיתה אנו עובדים. לשם כך נצהיר קודם על norb ו-nelec, שהם מספר האורביטלים ומספר האלקטרונים בהתאמה. אחר כך נצהיר על h1e ו-h2e, שהם האינטגרלים החד-אלקטרוניים והדו-אלקטרוניים בהתאמה. את כולם נשתמש מאוחר יותר גם ב-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

טעינת ההמילטוניאן

כשהנתונים הדרושים מוכנים, אנו קוראים את ההמילטוניאן מקובץ ה-FCI בפורמט התואם ל-qiskit-fermions

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

זרימות עבודה פרמיוניות עם qiskit-fermions

נמפה קודם את ההמילטוניאן למודל מעגל פרמיוני באמצעות qiskit-fermions, שמספק מעברי Transpiler ושערים ייעודיים למעגלים פרמיוניים. אלה ישמשו מאוחר יותר לפני מעברי ה-Transpiler המסורתיים של Qiskit בזרימת העבודה הזו.

קיבוץ איברים

כדי להבטיח שחזוריות של התוצאות, אנו משתמשים קודם ב-canonical_order כדי למיין את האיברים לפי המבנה שלהם בלבד. לכן סדר האופרטורים ברשימה canon קבוע. זה מבטיח שחזוריות של האופרטורים שנוצרו, כי המעבר QDriftTrotterization שבו נשתמש בהמשך דוגם אינדקסים אקראיים כדי ליצור את אופרטורי qDRIFT.

בשלב הזה אנו מנצלים את הסימטריות הרבות שקיימות בהמילטוניאן של המבנה האלקטרוני על ידי קיבוץ איברים קשורים עם מקדמים זהים. קיבוץ כזה משנה את התפלגות המקדמים של האופרטור שממנה qDRIFT דוגם, אך הוא לא פוגע בערבויות ההתכנסות שלו. חשוב לציין שקיבוץ איברים הקשורים בסימטריה גורם לביטול רצוי של איברי פאולי ולעומק מעגל כולל קצר יותר בעת התפתחות בזמן של מצב תחת פעולתם.

qiskit-fermions מספקת את הפונקציה group_terms_by_electronic_structure שעושה בשבילנו את הקיבוץ הזה.

שימו לב ש-group_terms_by_electronic_structure מניחה איברים ב-סדר נורמלי.

סינון איברים אלכסוניים

אנו מסירים את האיברים האלכסוניים מההמילטוניאן שמשמש ליצירת המעגלים, כדי שמשבצות הדגימה $n$ של qDRIFT יוקדשו לאיברים שמזיזים אוכלוסייה בין תצורות. עדיף לסנן איברים כאלה מההמילטוניאן בנקודה הזו, לפני שבונים את השער Evolution בשלב הבא.

האיברים הנדונים הם אלה שאלכסוניים בבסיס מספרי התפוסה, כלומר מכפלות של אופרטורי מספר $a^\dagger_i a_i$ . שלושה סוגי איברים נכללים בתיאור הזה:

היסט האנרגיה הקבוע, מכפלה של אפס אופרטורי מספר, שההתפתחות בזמן שלו תורמת רק פאזה גלובלית;
אופרטורי המספר הבודדים $n_i$ , שההתפתחות בזמן שלהם מצטמצמת לסיבובי $Z$ של קיוביט יחיד;
מכפלות מסדר גבוה יותר כמו $n_i n_j$ .

בפני עצמם, אף אחד מהם לא מזיז אוכלוסייה בין תצורות של מספרי תפוסה; הם פועלים רק על הפאזות של התצורות שכבר קיימות. הם אינם אינרטיים, עם זאת: הפאזות היחסיות האלה מזינות את האינטרפרנציה שנוצרת על ידי איברי העירור בהמשך המעגל, ולכן סינונם משנה את ההתפתחות שנוצרת בפועל ועשוי לשנות את התפלגות הדגימה. זו קירוב מכוון בשלב יצירת המעגלים, שנועד למקד את הדגימה באיברי עירור, ולא שלב שמשאיר את ההתפלגות הנדגמת ללא שינוי. בניגוד לקיבוץ הסימטריה שלמעלה, שמשאיר את ערבויות ההתכנסות של qDRIFT שלמות, המסנן הזה משנה את האופרטור שמתפתח. לכן המעגלים כבר לא מקרבים התפתחות תחת ההמילטוניאן המלא, וחסמי השגיאה של qDRIFT חלים על האופרטור המסונן ולא על המקורי. זה מקובל כאן כי המעגלים הם רק היוריסטיקת דגימה שמשמשת להצעת תצורות: שום איבר לא אובד מהערכת האנרגיה עצמה, כי המסנן חל רק על ההמילטוניאן שמשמש לבניית המעגלים, בעוד שהאלכסון הקלאסי בהמשך משתמש בהמילטוניאן המלא, כולל האיברים האלכסוניים. הדיוק של SQD תלוי בשלב הקלאסי הזה, שנשאר וריאציוני בתת-המרחב שנדגם בלי קשר לאופן שבו הוצעו התצורות.

הפונקציה filter_diagonal_terms() מסירה איברים כאלה מאופרטור במקום. היא מזהה אותם לפי המבנה הנורמלי-מסודר שלהם — הרב-קבוצה של מודי היצירה זהה לרב-קבוצה של מודי ההשמדה — ולכן היא תקפה רק על אופרטור שכבר בסדר נורמלי. ההנחה הזו לא נבדקת בזמן ריצה.

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

עכשיו, אחרי שקיבצנו את האיברים בהמילטוניאן, נחליט על הפרמטרים הבאים כדי ליצור את אנסמבל המעגלים:

מספר המעגלים ליצירה: num_circuits
אורך כל מעגל במונחי קבוצות עירור: num_exc
הגורם לזמני ההתפתחות השונים: times

יצירת מעגלים פרמיוניים

עכשיו ניצור מעגלים פרמיוניים לכל אחד מצעדי הזמן. כל מעגל יכלול שער התפתחות יחיד, עם זמן ההתפתחות שהצהרנו עליו קודם. אופרטור ההתפתחות הוא ההמילטוניאן. בהמשך נריץ על המעגלים האלה מעברי Transpiler כדי ליצור מעגלי qDRIFT.

הכנת ansatz

אנו מכינים את מצב Hartree-Fock באמצעות המחלקה InitializeModes. עבור חנקן, התהליך הוא פשוט הפעלת שערי X על num_elec_a הקיוביטים הראשונים ואז על num_elec_b הקיוביטים, ששניהם שווים לשבעה עבור חנקן. המצב הזה מייצג את שבעת אלקטרוני $\alpha$ ושבעת אלקטרוני $\beta$ של חנקן.

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

שלב 2: אופטימיזציה של הבעיה להרצה בחומרה קוונטית

עכשיו, כשיש לנו את המעגלים, נשתמש קודם במעברים הזמינים ב-qiskit-fermions כדי לבצע אופטימיזציות ברמה הפרמיונית, ואחר כך נבצע טרנספילציה של המעגל ל-Backend שבחרנו. מכיוון שזה ניסוי עם סימולטור, נעשה זאת קודם עבור AerSimulator. חישוב משקל לכל קבוצה

בשלב הזה אנחנו מבצעים דגימת qDRIFT של איברים באופן סטוכסטי, בהסתברויות שפרופורציונליות למקדמים שלהם בהמילטוניאן. מעבר ה-Transpiler של qDRIFT עושה את זה בשבילנו. עכשיו אפשר ליצור מעגלים רדודים יותר, שאפשר להריץ בצורה יעילה יותר על החומרה למרות קישוריות הקיוביטים המוגבלת, גם כשההמילטוניאן מכיל צימודים ארוכי טווח ואיברים מסדר גבוה מריבועי. לאחר קיבוץ האיברים, הדגימה של האופרטורים נעשית לפי המשקלים שלהם. לכל אופרטור $h_i$ , המשקל $W_{h_i}$ מוגדר כך:

W_{h_i} = |c_i| / \lambda

מכיוון שהאיברים קובצו בשלב 1, כל $h_i$ כאן הוא קבוצה שלמה: $c_i$ הוא הממוצע של הערכים המוחלטים של המקדמים של האיברים בקבוצה $i$ , וכל איבר בקבוצה מתפתח כשהמקדם שלו מצטמצם לסימן שלו בלבד.

אופטימיזציות פרמיוניות ואופטימיזציות מקוריות לחומרה

הפונקציה generate_preset_jw_pass_manager() מחזירה MultiStagePassManager שמקבל FermionicCircuit ומפיק מעגל סופי מותאם, שאפשר לבצע לו transpile כדי להריץ אותו על החומרה שלנו. אנחנו מחליפים את שלב האופטימיזציה ברירת המחדל שלו ב-FermionicPassManager שמכיל את המעבר QDriftTrotterization שלנו:

המעבר QDriftTrotterization משתמש בחישוב המשקלים ובדגימה באופן פנימי כדי ליצור את המעגלים שבהם נשתמש לדגימה
המעבר RelabelModes הוא מעבר אופטימיזציה נוסף שאפשר להשתמש בו כדי לבצע פרמוטציה של המודים הפרמיוניים, כדי לשפר את הקישוריות בין הקיוביטים ולהקטין את עומק השערים; למידע נוסף ראו את התיעוד של ה-API

שאר השלבים של MultiStagePassManager רצים אוטומטית ומטפלים במיפוי המלא מפרמיונים לקיוביטים:

F2QLayout: מנהל המעברים המוגדר מראש מפעיל את המעבר TrivialF2QLayout, שממפה באופן טריוויאלי $n$ ביטים פרמיוניים ל- $n$ קיוביטים.
F2QSynth: מעבר transpilation שממפה הוראות של מעגלים מבוססי פרמיונים להוראות מבוססות קיוביטים.

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

עכשיו, אחרי שסיימנו את האופטימיזציות ברמה הפרמיונית, אפשר לבצע transpile למעגלים כדי להריץ אותם על הסימולטור.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

שלב 3: הרצה באמצעות פרימיטיבים של Qiskit

עכשיו, כשיש לנו את המעגלים, אפשר להריץ אותם באמצעות פרימיטיבים של Qiskit על AerSimulator. נשלב את כל הספירות ממעגלים שונים. נמיר אותן לוקטורים בוליאניים, ורק אז נבצע עיבוד-המשך באמצעות 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

שלב 4: עיבוד-המשך והחזרת התוצאה בפורמט הקלאסי הרצוי

שימוש ב-bitstrings עבור SQD

עכשיו אפשר להריץ את תכנית האלכסון על ה-bitstrings שנבחרו כדי למצוא את הערך העצמי הנמוך ביותר, שמתאים לאנרגיית מצב היסוד של המולקולה. אנחנו יוצרים פונקציית callback, מגדירים תפוסות התחלתיות וקובעים את הפרמטרים, ורק אז מריצים את תכנית האלכסון. פונקציית ה-callback משמשת להדפסת האיטרציה הנוכחית ואומדן הערך העצמי הנוכחי בכל איטרציה.

לבסוף, כדי לקבל את אומדן מצב היסוד, אנחנו מוסיפים את nuclear_repulsion_energy לאנרגיה שהתקבלה.

הערה: ממד תת-המרחב אינו קבוע בין האיטרציות, גם בסימולטור נטול רעש — כל תת-דגימה שולפת קבוצה אחרת של קונפיגורציות, ושלב השחזור משנה את המאגר בין האיטרציות, ולכן הממד המדווח משתנה מתת-דגימה אחת לאחרת. דגימה נטולת רעש לבדה אינה קובעת את ממד תת-המרחב הנבחר. לעומת זאת, בהרצה על חומרה נוטים לקבל תתי-מרחבים גדולים יותר באופן שיטתי, כי shots רועשים שוברים את סימטריית מספר החלקיקים, ושחזור הקונפיגורציות הופך אותם לווקטורי בסיס נוספים. בגלל זה נציג בחלק של החומרה גם שלב נוסף של גיזום bitstrings.

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

דוגמה על חומרה

הדוגמה הזו משתמשת ב-20 קיוביטים (10 אורביטלים מרחביים). הבחירה הזו נוחה לשיעור שאמור לרוץ מהר, והיא לא תקרה קשיחה לשיטה.

העלות של השלב הקלאסי לא נקבעת ישירות לפי מספר הקיוביטים. SQD מבצע אלכסון של ההמילטוניאן המוטל על תת-המרחב שנפרש על ידי הקונפיגורציות השנדגמו, ולכן מה שקובע את העלות הקלאסית הוא הממד של תת-המרחב שנבחר — שנקבע כאן לפי samples_per_batch, num_batches וכמה קונפיגורציות שונות המעגלים באמת מייצרים — יחד עם האלגברה הלינארית הדלילה שנדרשת כדי להפעיל את ההמילטוניאן המוטל. מרחב ה-CI המלא גדל קומבינטורית עם מספר האורביטלים והאלקטרונים, אבל תת-המרחב שנבחר הוא פרוסה קטנה וניתנת לכוונון ממנו, ואנחנו שולטים בגודל שלו ישירות. לכן אפשר לשנות את מספר הקיוביטים ואת הקושי הקלאסי באופן בלתי תלוי במידה מסוימת: מרחב אורביטלים רחב יותר שנדגם לתוך תת-מרחב צנוע יכול להיות זול יותר ממערכת קטנה יותר שמבצעים לה אלכסון על פני תת-מרחב גדול מאוד.

בפועל, גודל המערכת האפשרי תלוי בממד תת-המרחב שנדרש לדיוק הרצוי, ובזיכרון ובליבות שזמינים לפותר הערכים העצמיים. מרחבי אורביטלים גדולים יותר דורשים בדרך כלל תת-מרחב גדול יותר כדי להגיע לדיוק כימי, וזה מה שבסופו של דבר מניע לשימוש במשאבים מבוזרים — ראו qiskit-addon-sqd-hpc להרחבת השלב הזה. במקום להניח סף קבוע, הגישה המעשית היא לעקוב אחרי ממד תת-המרחב המדווח ואחרי התכנסות האנרגיה לאורך האיטרציות, ולהגדיל את תת-המרחב עד שהאנרגיה מפסיקה להשתפר או עד שהזיכרון הזמין אוזל.

הערה: בגלל שגיאת דגימה שנובעת מהרעש בחומרה, תת-המרחב שנוצר לאלכסון בהרצה על החומרה יהיה גדול יותר ממה שמתקבל בסימולטור. זה מגדיל את הממד של תת-המרחב שרוצים לבצע לו אלכסון, אבל תהליך העבודה עדיין נותן תשובה מדויקת בזכות העמידות של SQD לרעש.

גיזום של מחרוזות מזויפות

כאן אפשר לבחור לבצע שלב נוסף. אחרי שיש לנו את כל ה-bitstrings מהרצות המעגלים, אפשר לסנן את ה-bitstrings הלא תקינים לפני הרצת SQD, או להמשיך בלי גיזום. בהרצות על חומרה עדיף בדרך כלל לדלג על הגיזום, כי כך ה-shots ששברו סימטריה נשארים זמינים לשחזור הקונפיגורציות, שיכול לתקן אותם לקונפיגורציות תקינות ובכך להרחיב את תת-המרחב במקום לזרוק את ה-shots האלה לגמרי.

מכיוון שלחנקן יכולים להיות רק שבעה אלקטרונים מסוג $\alpha$ ושבעה מסוג $\beta$ , אפשר להשליך כל bitstring שיש בו יותר או פחות משבע ספרות 1 בחצי הראשון או בחצי השני של הפלט. אנחנו מגדירים פונקציה שבודקת אם ה-bitstrings תקינים, ואם לא, משליכה אותם. אחרי שמסננים את ה-bitstrings המזויפים, השאר נשלחים לתכנית האלכסון. השתמשו בדגל PRUNE שלמטה כדי לעבור בין שתי ההתנהגויות.

חשוב לזכור שגיזום הוא רק אחת מכמה בחירות שמעצבות את תת-המרחב הסופי, לצד מספר המעגלים, קבוצת זמני ההתפתחות וסינון איברים אלכסוניים. השוואה בין הרצה עם גיזום להרצה בלי גיזום מועילה רק אם כל השאר נשאר קבוע; הגרסה ב-C++ של השיעור הזה דנה בזה בפירוט רב יותר, כי היא מבצעת postselection במקום שחזור, וגם שונה בפרמטרים האחרים האלה.

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

השלבים הבאים

המלצות

אם העבודה הזו עניינה אתכם, ייתכן שגם החומרים הבאים יעניינו אתכם:

Sample-based Krylov quantum diagonalization of a fermionic lattice model - שיעור קשור שמשתמש במעגלי התפתחות בזמן במקום ansatz וריאציוני.
Sample-based quantum diagonalization of a chemistry Hamiltonian - שיעור על בניית מעגל local unitary cluster Jastrow (LUCJ) לסימולציה כימית קוונטית.
המאמר SqDRIFT - הספרות שהשיעור הזה מבוסס עליה. (שימו לב שחלק מהאופטימיזציות שנדונות במאמר הזה הן כרגע עבודה בתהליך, והשיעור הזה עשוי להשתנות בעתיד בהתאם להתפתחות הספריות שבהן נעשה שימוש.)

מטרות הלמידה​

דרישות קדם​

רקע​

דרישות​

הגדרה​

דוגמה עם סימולטור​

שלב 1: מיפוי קלטים קלאסיים לבעיה קוונטית​

שלב 2: אופטימיזציה של הבעיה להרצה בחומרה קוונטית​

שלב 3: הרצה באמצעות פרימיטיבים של Qiskit​

שלב 4: עיבוד-המשך והחזרת התוצאה בפורמט הקלאסי הרצוי​

דוגמה על חומרה​

השלבים הבאים​