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

Python

הורידו את מחברת Python

הורידו את מחברת Python כדי להריץ את המדריך בעצמכם. המחברת מורידה את קלט ה-FCIDump אם הוא עדיין לא קיים.

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

תוצאות למידה

• למדו כיצד ליצור מעגלים בעומק קטן יותר בהשוואה ל-Trotterization

• עברו על זרימת עבודה מקצה לקצה להערכת מצב יסוד באמצעות qDRIFT ו-SQD

• למדו כיצד להשתמש ב-qiskit-fermions יחד עם תוספי Qiskit אחרים כדי ליישם זרימת עבודה כזו

דרישות קדם

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

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

רקע

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

שני המימושים משתמשים ב-Qiskit Fermions כדי לייצג את ההמילטוניאן המולקולרי ולבנות מעגלים עבור אלגוריתם ה-qDRIFT. מימוש ה-Python משתמש במעברי (passes) transpiler פרמיוניים, ואילו מימוש ה-C++ דוגם קבוצות במפורש וממפה אותן לאופרטורי קיוביטים.

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

$$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$-י, המוגדר כ:

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

הגדרה

# 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: מיפוי קלטים קלאסיים לבעיה קוונטית

שלבי ה-Python שלהלן משתמשים בסימולטור Aer. דוגמת החומרה באה אחריהם.

קריאה והכנה של ה-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}")

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

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, המספקת מעברי טרנספיילר ושערים ייעודיים למעגלים פרמיוניים. אלה ישמשו מאוחר יותר לפני מעברי הטרנספיילר המסורתיים של 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))

5060

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

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

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

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

הכנת האנזץ

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

• מעבר ה-QDriftTrotterization משתמש בחישוב המשקלים ובדגימה באופן פנימי כדי ליצור את המעגלים שבהם נשתמש לדגימה

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

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

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

• F2QSynth: מעבר טרנספילציה למיפוי הוראות מעגל מבוססות-פרמיונים להוראות מבוססות-קיוביטים.

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

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

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: עיבוד סופי והחזרת התוצאה בפורמט הקלאסי הרצוי

שימוש במחרוזות ביטים עבור SQD

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

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

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

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 כלפי רעש.

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

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

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

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

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

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

C++

קבלו את קוד ה-C++

קוד המקור ב-C++, פרויקט ה-CMake והוראות הבנייה נמצאים בתיקייה docs/tutorials/assets/sqdrift/cpp של מאגר התיעוד של Qiskit. קלט ה-FCIDump נמצא בתיקייה השכנה fcidump_files. שכפלו את המאגר, או הורידו את שתי התיקיות, כדי לבנות ולהריץ את המדריך בעצמכם.

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

תוצאות למידה (C++)

• למדו כיצד ליצור מעגלים בעומק קטן יותר בהשוואה ל-Trotterization

• עברו על זרימת עבודה מקצה לקצה להערכת מצב יסוד באמצעות qDRIFT ו-SQD

• למדו כיצד להשתמש ב-qiskit-fermions יחד עם תוספי Qiskit אחרים כדי ליישם זרימת עבודה כזו

דרישות קדם (C++)

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

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

רקע (C++)

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

שני המימושים משתמשים ב-Qiskit Fermions כדי לייצג את ההמילטוניאן המולקולרי ולבנות מעגלים עבור אלגוריתם ה-qDRIFT. מימוש ה-Python משתמש במעברי (passes) transpiler פרמיוניים, ואילו מימוש ה-C++ דוגם קבוצות במפורש וממפה אותן לאופרטורי קיוביטים.

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

$$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$-י, המוגדר כ:

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

דרישות (C++)

המימוש הזה רץ על חומרת IBM Quantum® ומשתמש בקובץ הרצה חיצוני של SBD לאלכסון קלאסי. הוא אינו כולל את דוגמת הסימולטור ב-Python.

עקבו אחר הוראות הבנייה של C++ כדי להתקין את המהדר, CMake, Rust, Boost ואת התלויות של Qiskit. התקינו את Rust באמצעות rustup: התלויות של Qiskit נועלות גרסת Rust משלהן, ש-rustup מוריד אוטומטית, אבל Rust שהותקן מ-Homebrew או ממנהל חבילות של Linux עלול להיות ישן מדי. שלב ה-SBD האחרון דורש בנוסף MPI, OpenMP ו-BLAS/LAPACK. Windows אינו נתמך כרגע על ידי המימוש הזה.

השתמשו בקוד המקור המלא של C++ ובפרויקט ה-CMake מהמאגר. שמרו את התיקייה הסמוכה fcidump_files בעת הורדת הדוגמה. מתוך עותק של מאגר התיעוד, בנו עם הפקודות הבאות. הבנייה הראשונה מורידה ומקמפלת את כל התלויות ולוקחת כמה דקות.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

הבנייה כותבת את קובץ ההרצה SqDRIFT אל docs/tutorials/assets/sqdrift/cpp (שורש הפרויקט), ולא אל build/. הריצו ./SqDRIFT משורש הפרויקט כדי שיוכל למצוא את ../fcidump_files/N2_sto_3g. כל שאר הפקודות בלשונית הזו רצות גם הן משורש הפרויקט אלא אם נאמר אחרת. אם ./SqDRIFT נכשל בשגיאת ספרייה משותפת חסרה, עקבו אחר ההוראות לנתיב הספריות לפי פלטפורמה במדריך הבנייה; ב-macOS, קובץ ההרצה בדרך כלל רץ בלעדיהן.

הגדרה (C++)

הקטעים שלהלן הם מובאות מתוך תוכנית SqDRIFT.cpp המלאה. בנו והריצו את הקובץ הזה כתוכנית אחת.

לפני הרצת SqDRIFT, ודאו שאסימון ה-API של IBM Quantum שלכם מוגדר עבור IBM Cloud Platform.

אם מותקן אצלכם Python עם qiskit-ibm-runtime, הדרך הקלה ביותר ליצור קובץ תצורה בפורמט נכון היא:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

זה כותב את ~/.qiskit/qiskit-ibm.json (Linux/macOS) או %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) בצורה הנכונה. אם אתם מעדיפים ליצור את הקובץ ידנית:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

קבלו את מפתח ה-API שלכם מ-IBM Quantum Platform.

הערה: קובץ תצורה חסר או פגום גורם ל-panic לא מועיל של unwrap() ב-Rust בעת ההפעלה במקום שגיאה נקייה. בדקו שוב את שם המפתח (default-ibm-quantum-platform) ואת ערך הערוץ (ibm_quantum_platform) אם נתקלתם בבעיה.

מימוש (C++)

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

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

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

טוען את ההמילטוניאן של מולקולת $N_2$ מקובץ FCIDump. מספר הקיוביטים שווה לכפל מספר האורביטלים (ספין למעלה וספין למטה).

יצירת ההמילטוניאן וסידורו הנורמלי

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

ממיר את נתוני ה-FCIDump לאופרטור פרמיוני ומחיל סידור נורמלי כדי לפשט את מבנה האופרטור.

קיבוץ איברים לפי מבנה אלקטרוני

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

מקבץ איברי המילטוניאן קשורים לפי מבנה אלקטרוני לפני בניית מעגלי האבולוציה שלהם.

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

הסדר חשוב כאן. הנרמול משנה את קנה המידה של כל מקדם לגודל יחידה במקום, ולכן משקלי הדגימה חייבים להיחשב קודם. אם הנרמול רץ ראשון, כל $|c_j|$ נקרא בחזרה כ-1.0 בדיוק, והמשקל של כל קבוצה קורס לספירה גרידא של האיברים שלה, וההתפלגות הדוגמת הייתה מונעת אז על ידי גודל הקבוצה ולא על ידי המקדמים הפיזיקליים.

המשקל של כל קבוצה הוא הערך המוחלט הממוצע של מקדמי האיברים שלה. מעבר QDriftTrotterization בלשונית Python משתמש באותו משקל להמילטוניאן מקובץ. שימוש בסכום במקום זאת היה גורם ל-$\lambda$ לגדול עם גודל הקבוצה, מה שמעוות בתורו את זמן האבולוציה $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

המשקל $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ הוא הערך המוחלט הממוצע של המקדמים על פני איברי הקבוצה $G_i$, הנקרא מהמקדמים המקוריים. לאחר מכן הנרמול משנה את קנה המידה של כל איבר פרמיוני לגודל יחידה לפני התמרת Jordan-Wigner, ושומר באופרטור המתפתח רק את הפאזה/הסימן שלו. זה תואם את האבולוציה של איברים מנורמלים ב-qDRIFT: גדלי המקדמים קובעים את הסתברויות הדגימה דרך $\lambda$, ואילו אבולוציית המעגל משתמשת באיבר המנורמל כדי שמקדמים גדולים לא ייספרו פעמיים.

דגימת קבוצות אופרטורים

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

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

שלב 2: בניית מעגלי qDRIFT (C++)

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

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

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

עבור כל אופרטור שנדגם, יוצר מעגל אבולוציית Suzuki-Trotter עם זמן $t_k = \frac{\lambda \tau}{n}$, כאשר $n$ הוא מספר האופרטורים שנדגמים לכל מעגל ו-$\tau$ הוא צעד הזמן. מכיוון שכל קבוצה שנדגמה נורמלה איבר-איבר לפני המיפוי, האבולוציה הזו מפעילה רק את הפאזה/הסימן של כל איבר במהלך הסיבוב; הגדלים תורמים רק דרך התפלגות הדגימה. מצב ההתחלה של Hartree-Fock מוכן על ידי הפעלת שערי X על קיוביטים 0-6 ו-10-16 (שבעה אלקטרונים בכל מקטע ספין). ההוראות של כל מעגל אבולוציה מתווספות ידנית למעגל הראשי באמצעות ה-API של C++.

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

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

מתחבר ל-IBM Quantum, בוחר את ה-backend לפי שם, מבצע transpile לכל מעגל עבור ה-backend הזה, ושולח משימה אחת לכל מעגל: 100 משימות Sampler נפרדות של 100 shots כל אחת על אוגר של 20 קיוביטים. לאחר מכן הוא בודק כל משימה שוב ושוב עד שהיא מסתיימת ואוסף את מחרוזות הסיביות מ-result[0]. כל משימה משתמשת בזמן QPU וממתינה בתור בנפרד, ולכן צפו שהשלב הזה ייקח הרבה יותר זמן ממשימה אחת מקובצת.

תוצאות החומרה שמצוטטות בעמוד הזה מגיעות מהרצה על ibm_pittsburgh (100 מעגלים × 100 shots), כאשר backend_name שונה מברירת המחדל ibm_fez. הדגימות השמורות עובדו עם סדר הסיביות המתוקן של הדטרמיננטות שמוצג בשלב 4. רעש החומרה משתנה בין הרצות, ולכן הספירות והאנרגיה הסופית שלכם יהיו שונות.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

שינוי ה-backend: שנו את backend_name ב-SqDRIFT.cpp כדי לכוון להתקן אחר; service.backends() מציג את מה שהחשבון שלכם יכול להגיע אליו.

שלב 4: עיבוד סופי והחזרת התוצאה בפורמט הקלאסי הרצוי (C++)

Postselection של מחרוזות סיביות וכתיבת מחרוזות CI

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

מבצע postselection למחרוזות סיביות עם משקל המינג נכון (n_alpha אלקטרונים עם ספין למעלה, n_beta עם ספין למטה, הנגזרים משדה NELEC של FCIDump), ממיר אותן למחרוזות Configuration Interaction (CI) באמצעות סימטריזציית ספין, וכותב אותן ל-alphadets_from_sqd.txt עבור Selected Basis Diagonalization הבא.

ניקוי משאבים

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

משחרר כראוי את כל הזיכרון והמשאבים שהוקצו. ה-std::vector<QuantumCircuit> של C++ מנוקה אוטומטית באמצעות RAII. משחרר ידנית את אופרטורי הקיוביטים, האופרטורים הפרמיוניים ונתוני ה-FCIDump.

עכשיו, כשיצרנו את הבסיס שעליו נטיל את ההמילטוניאן, אפשר להמשיך בתהליך הדיאגונליזציה כדי לקבל את הערכת מצב היסוד. הריצו את הפקודות הבאות משורש הפרויקט, docs/tutorials/assets/sqdrift/cpp.

בניית הקובץ הבינארי diag של SBD

CMake מביא את SBD כקוד מקור בלבד אל deps/sbd/ (הבנייה הראשית מדלגת בכוונה על הקומפילציה שלו). יישום ה-SBD משתמש ב-Makefile פשוט, ו-SBD עצמו דורש MPI, OpenMP ו-BLAS/LAPACK — ולכן make פשוט אינו מספיק עד שאלה מותקנים וקובץ Configuration של היישום תואם לשרשרת הכלים שלכם.

ה-Makefile קורא את דגלי המהדר והקישור שלו מקובץ Configuration באותה תיקייה. הגרסה שמסופקת במקור מכוונת לפריסה אחת ספציפית של macOS/Homebrew ולא תיבנה כמו שהיא על מחשב רגיל — בעיקר -I/opt/homebrew/opt/llvm/include אינו מכיל את omp.h (הוא נמצא בתיקיית משאבים של clang שתלויה בגרסה), ו-clang++ של Apple דוחה את -fopenmp לחלוטין:

clang++: error: unsupported option '-fopenmp'

המתכון הבא אומת כעובד ב-macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — נבדק

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

לאחר מכן החליפו את deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration ב:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

בנייה מוצלחת מדפיסה את שתי שורות הקומפילציה/הקישור ומשאירה קובץ בינארי diag ליד ה-Makefile. אזהרה אחת של ld: warning: ignoring duplicate libraries: '-lomp' אינה מזיקה.

macOS עם Intel: הקידומת של Homebrew היא /usr/local ולא /opt/homebrew, ולכן החליפו אותה בשני הנתיבים הקודמים.

Linux (Ubuntu/Debian) — לא נבדק

דרישות הקדם הבאות מכסות MPI, OpenMP ו-BLAS/LAPACK. GCC תומך ב--fopenmp באופן מובנה, ולכן אין צורך בדריסת מהדר. זה לא אומת ב-Linux:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

הרצת SqDRIFT ליצירת מחרוזות CI

אם כבר הרצתם את ./SqDRIFT דרך שלבים 1–4, הקובץ alphadets_from_sqd.txt כבר קיים בשורש הפרויקט ואפשר לדלג על השלב הזה. אחרת, הריצו אותו משורש הפרויקט:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

הרצת הדיאגונליזציה

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

פרמטרים:

• --fcidump: נתיב לקובץ FCIDump שמכיל את ההמילטוניאן
• --adetfile: נתיב לקובץ דטרמיננטות אלפא (מחרוזות CI מ-SqDRIFT)
• --method 0: שיטת אלכסון Davidson
• --iteration 100: מספר איטרציות Davidson מרבי
• --block 10: גודל בלוק לאלגוריתם Davidson
• --tolerance 1e-8: סבילות התכנסות לאנרגיה

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

שורת Energy = -107.5108722590187 היא הערכת מצב היסוד על פני מרחב המשנה שנדגם.

במה זה שונה מדוגמת ה-Python

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

	לשונית C++	לשונית Python
מחרוזות ביטים עם סימטריה שגויה	מושלכות על ידי postselect_bitstrings לפי משקל המינג	מתוקנות על ידי שחזור קונפיגורציות איטרטיבי
מבנה הלולאה	מעבר יחיד: דגימה → postselection → דיאגונליזציה פעם אחת	לולאה חיצונית: דיאגונליזציה, קריאת אכלוסי אורביטלים, שחזור קונפיגורציות, דיאגונליזציה מחדש
מספר מעגלים	100 מעגלים (num_circuits = 100)	400 מעגלים (200 לכל זמן אבולוציה)
זמני אבולוציה	$\tau$ יחיד (time_step = 1)	שני זמנים, times = [1.0, 10.0]
איברים אלכסוניים	נשמרים באופרטור שנדגם	מוסרים על ידי filter_diagonal_terms, מה שמשנה את התפלגות הדגימה ואת המעגלים
סידור איברים	מקובצים בלבד	canonical_order מופעל לפני הקיבוץ לשחזוריות
פותר דיאגונליזציה	SBD (diag, קובץ בינארי חיצוני של MPI/OpenMP)	פותר qiskit-addon-sqd בתוך התהליך
Backend	קבוע לפי שם (backend_name = "ibm_fez"; ההרצה בעמוד הזה השתמשה ב-ibm_pittsburgh)	backend תפעולי מסוג least_busy עם מספיק קיוביטים
Transpilation	כל מעגל בנפרד, רמת אופטימיזציה 2	כל המעגלים יחד, רמת אופטימיזציה 3
שליחת משימות	100 משימות, מעגל אחד בכל אחת	משימת Sampler אחת שמכילה את כל 400 המעגלים
הרצת סימולטור	אין; חומרה בלבד	דוגמת סימולטור ללא רעש לפני הרצת החומרה

שתי הגרסאות משתמשות באותו המילטוניאן, מקבצות את איבריו לפי מבנה אלקטרוני, משקללות כל קבוצה לפי הערך המוחלט הממוצע של מקדמי איבריה, מתחילות ממצב Hartree-Fock, ודוגמות 10 קבוצות ו-100 shots לכל מעגל.

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

Postselection הוא מסנן קפדני, ו-shot עם מספר חלקיקים שגוי נזרק לחלוטין. בהרצה שלמעלה זה השליך את רוב הנתונים: 1137 מתוך 10000 shots שרדו את בדיקת משקל המינג (7,7) (כ-11%), ו-collapsו לכדי 49 דטרמיננטות $\alpha$ שונות בלבד. שחזור קונפיגורציות משתמש במקום זאת באכלוסי אורביטלים ממוצעים מאלכסון קודם כדי להפוך סיביות ולתקן shots כאלה, ומנצל מחדש את מה ש-postselection היה משליך. המשוב הזה הופך אותו מטבעו לאיטרטיבי: recover_configurations ב-qiskit-addon-sqd-hpc דורש ארגומנט avg_occupancies שקיים רק אחרי שאלכסון רץ.

הרצת C++ הזו מניבה מרחב משנה של 49 דטרמיננטות $\alpha$ ואנרגיה של -107.5109 Ha. מרחב דטרמיננטות ה-$\alpha$ המלא לשבעה אלקטרונים ב-10 אורביטלים הוא רק $\binom{10}{7} = 120$ מחרוזות, ולכן הזנת כל ה-120 לאותו קובץ בינארי diag נותנת את אנרגיית ה-FCI המדויקת -107.6482 Ha. לכן ההרצה נמצאת כ-0.1373 Ha מעל המדויק בעוד שהיא פורשת 49 מתוך 120 דטרמיננטות ה-$\alpha$. מכיוון ש-SQD הוא וריאציוני במרחב המשנה הנבחר, מרחב משנה קטן יותר יכול רק להעלות את האנרגיה, ולכן מרחב המשנה החלקי הוא ההסבר הישיר ביותר לפער.

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

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

כדי לסגור את הפער, הנה כמה צעדים אפשריים לנסות אחרי השלמת המדריך הזה: להגדיל את shots, num_circuits או ops_per_circuit; לדגום ביותר מ-time_step אחד; ולהוסיף לולאת שחזור מעל Qiskit::addon::sqd::recover_configurations.

צעדים הבאים

המלצות

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

• אלכסון קוונטי מבוסס דגימות של המילטוניאן כימי — מדריך על בניית מעגל local unitary cluster Jastrow (LUCJ) לסימולציה קוונטית של כימיה.

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