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

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

מחפשים את גרסת ++C?

מדריך זה משתמש ב-Python. עבור היישום ב-++C, כולל קוד מקור והוראות בנייה, ראו את מדריך ה-SqDRIFT ב-++C.

תוצאות למידה

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

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

דרישות קדם

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

רקע

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

מדריך זה משתמש ב-Qiskit Fermions כדי ליצור את המעגלים הפרמיוניים הטבעיים יותר עבור אלגוריתם ה-qDRIFT, ולאחר מכן שימוש במעברי (passes) פריסה וסינתזה פרמיוניים לפני חיבור המעגלים לצינור העבודה המסורתי של 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$ -י, המוגדר כ:

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 — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

דוגמת סימולטור

שלב 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

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

עם הנתונים הדרושים מוכנים, אנחנו קוראים את ההמילטוניאן מקובץ ה-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))

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

מספר המעגלים ליצירה: 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

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

הפונקציה 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))

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

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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

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

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

צעדים הבאים

המלצות

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

דיאגונליזציה קוונטית מבוססת-קריילוב מבוססת-דגימה של מודל סריג פרמיוני - מדריך קשור המשתמש במעגלי אבולוציית זמן במקום באנזץ וריאציוני.
דיאגונליזציה קוונטית מבוססת-דגימה של המילטוניאן כימי - מדריך על אופן בניית מעגל local unitary cluster Jastrow (LUCJ) לסימולציית כימיה קוונטית.
המאמר SqDRIFT - הספרות שעליה מבוסס המדריך הזה. (שימו לב שחלק מהאופטימיזציות הנדונות במאמר זה נמצאות כרגע בעבודה מתמשכת, והמדריך הזה עשוי להשתנות בעתיד בהתאם להתפתחות הספריות בהן נעשה שימוש.)

תוצאות למידה​

דרישות קדם​

רקע​

דרישות​

הגדרה​

דוגמת סימולטור​

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

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

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

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

דוגמת חומרה​

צעדים הבאים​