שימוש ב-trim SQD לבחירה טובה יותר של קונפיגורציות

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

המימוש של trim SQD מבוסס על מחקר משותף של ה-Cleveland Clinic Foundation, RIKEN ו-IBM®. אם אתם משתמשים בתכונה הזו במחקר שלכם, צטטו גם את המאמר. את הציטוט אפשר למצוא בקובץ CITATION.bib, עם מפתח הציטוט merz2026crossing.

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

הגדרת המולקולה

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

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import pyscf
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

WARN: Unable to to identify input symmetry using original axes.
Different symmetry axes will be used.

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

הרצת הלולאה

איטרציה אחת של trim SQD עושה את הדברים הבאים:

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

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

ה-sci_solver של הלולאה מבצע את שני הסבבים, וברירת המחדל שלו היא פותר PySCF של החבילה הזו.

מפעילים את trim SQD על ידי העברת TrimPolicy כך:

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    max_iterations=5,
    policy=TrimPolicy(trim_ratio=0.5, max_strings_per_trim=250),
    symmetrize_spin=True,
    max_dim=250,
    seed=rng,
)

sci_state = result.sci_state
print(
    f"Merged subspace: {len(sci_state.ci_strs_a)} x {len(sci_state.ci_strs_b)} CI strings"
)
print(f"Exact energy:     {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Merged subspace: 148 x 148 CI strings
Exact energy:     -109.04667177808024
Estimated energy: -109.034369147063

בדיקה אם הסינון עוזר

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

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

הגבלת תת-המרחב

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

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

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

ארבעה ארגומנטים מגבילים את זה, ושימושי להגדיר לפחות אחד מהם:

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

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

• carryover_threshold, carryover_ratio ו-max_carryover מגבילים את המחרוזות שמזינות את האיטרציה הבאה, ומוחלים על התת-מרחב הממוזג כולו ולא לכל אצווה. carryover_threshold הוא סף האמפליטודה ש-StandardPolicy משתמש בו, עם אותה ברירת מחדל; שני האחרים מגבילים את הכמות, מה שהסף לבדו לא עושה.

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

ההרצה שלמעלה מגדירה יחד max_strings_per_trim=250 ו-max_dim=250, כך שלא מה שכל אצווה תורמת ולא התת-מרחב שנבנה ממנו יכולים להתפרץ.

שימוש בפותר ערכים עצמיים חיצוני

אפשר להשתמש ב-trim SQD עם פותר ערכים עצמיים חיצוני על ידי העברת הארגומנט sci_solver:

from sbd.sbd_solver import create_sbd_solver

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    policy=TrimPolicy(trim_ratio=0.1),
    sci_solver=create_sbd_solver(sbd_config={"method": 0, "eps": 1e-10}),
    seed=rng,
)

הקטע הקודם דורש שהחבילה sbd-eigensolver תותקן מ-PyPI.

חלוקת העבודה

עם TrimPolicy, הלולאה קוראת ל-sci_solver פעמיים בכל איטרציה: פעם אחת עם כל האצוות, ופעם אחת עם התת-מרחב הממוזג. שתי הקריאות עוברות דרך הממשק הרגיל של sci_solver, שהתמיכה הרב-תהליכית של החבילה הזו מפעילה באופן קולקטיבי בכל דרג MPI.

שתי רמות של מקביליות זמינות פוטנציאלית, והן בלתי תלויות:

• בתוך אלכסון. פותר כמו SBD יכול לפזר אלכסון בודד על פני כמה דרגים.

• על פני האצוות. פותר יכול, באופן עקרוני, לאלכסן את האצוות במקביל.

בהינתן קבוצה של אצוות, הפותר קובע איך להאציל עבודה לתהליכי העבודה שלו.