דלג לתוכן הראשי

דיאגונליזציה קוונטית מבוססת דגימות של המילטוניאן כימי

אומדן זמן שימוש: פחות מדקה על מעבד Heron r2 (הערה: זהו אומדן בלבד. זמן הריצה שלך עשוי להשתנות.)

רקע​

במדריך זה, נראה כיצד לעבד במעבד-הבא דגימות קוונטיות רועשות כדי למצוא קירוב למצב היסוד של מולקולת החנקן N2\text{N}_2 באורך קשר שיווי משקל, תוך שימוש באלגוריתם דיאגונליזציה קוונטית מבוססת דגימות (SQD), עבור דגימות שנלקחו ממעגל קוונטי בן 59 קיוביטים (52 קיוביטי מערכת + 7 קיוביטי עזר). מימוש מבוסס Qiskit זמין בהרחבות Qiskit של SQD, פרטים נוספים ניתן למצוא בתיעוד המתאים עם דוגמה פשוטה כדי להתחיל. SQD היא טכניקה למציאת ערכים עצמיים וווקטורים עצמיים של אופרטורים קוונטיים, כמו המילטוניאן של מערכת קוונטית, תוך שימוש במחשוב קוונטי וקלאסי מבוזר ביחד. מחשוב קלאסי מבוזר משמש לעיבוד דגימות שהתקבלו ממעבד קוונטי, ולהקרנה ודיאגונליזציה של המילטוניאן יעד בתת-מרחב שנפרש על ידן. זה מאפשר ל-SQD להיות עמיד לדגימות שנפגמו על ידי רעש קוונטי ולהתמודד עם המילטוניאנים גדולים, כמו המילטוניאנים כימיים עם מיליוני איברי אינטראקציה, מעבר להישג ידם של כל שיטות הדיאגונליזציה המדויקות. לזרימת עבודה מבוססת SQD יש את השלבים הבאים:

  1. בחר ansatz מעגלי והחל אותו על מחשב קוונטי על מצב התייחסות (במקרה זה, מצב Hartree-Fock).
  2. דגום מחרוזות ביטים ממצב קוונטי זה.
  3. הרץ את תהליך שחזור התצורה העצמי-עקבי על מחרוזות הביטים כדי לקבל את הקירוב למצב היסוד.

ידוע ש-SQD עובד היטב כאשר מצב העצמי היעד הוא דליל: פונקציית הגל נתמכת בקבוצה של מצבי בסיס S={∣x⟩}\mathcal{S} = \{|x\rangle \} שגודלה אינו גדל באופן אקספוננציאלי עם גודל הבעיה.

כימיה קוונטית​

תכונות של מולקולות נקבעות במידה רבה על ידי מבנה האלקטרונים בתוכן. כחלקיקים פרמיוניים, ניתן לתאר אלקטרונים באמצעות פורמליזם מתמטי הנקרא קוונטיזציה שניה. הרעיון הוא שישנם מספר אורביטלים, שכל אחד מהם יכול להיות ריק או תפוס על ידי פרמיון. מערכת של NN אורביטלים מתוארת על ידי קבוצה של אופרטורי השמדה פרמיוניים {a^p}p=1N\{\hat{a}_p\}_{p=1}^N שמקיימים את יחסי האנטי-קומוטציה הפרמיוניים,

a^pa^q+a^qa^p=0,a^pa^q†+a^q†a^p=δpq.\begin{align*} \hat{a}_p \hat{a}_q + \hat{a}_q \hat{a}_p &= 0, \\ \hat{a}_p \hat{a}_q^\dagger + \hat{a}_q^\dagger \hat{a}_p &= \delta_{pq}. \end{align*}

הצמוד a^p†\hat{a}_p^\dagger נקרא אופרטור יצירה.

עד כה, ההצגה שלנו לא לקחה בחשבון ספין, שהוא תכונה יסודית של פרמיונים. כאשר לוקחים בחשבון ספין, האורביטלים באים בזוגות הנקראים אורביטלים מרחביים. כל אורביטל מרחבי מורכב משני אורביטלי ספין, אחד שמסומן "ספין-α\alpha" ואחד שמסומן "ספין-β\beta". לאחר מכן אנו כותבים a^pσ\hat{a}_{p\sigma} עבור אופרטור ההשמדה המשויך לאורביטל הספין עם ספין σ\sigma (σ∈{α,β}\sigma \in \{\alpha, \beta\}) באורביטל המרחבי pp. אם ניקח NN להיות מספר האורביטלים המרחביים, אז יש סך של 2N2N אורביטלי ספין. מרחב הילברט של מערכת זו נפרש על ידי 22N2^{2N} וקטורי בסיס אורתונורמליים המסומנים במחרוזות ביטים דו-חלקיות ∣z⟩=∣zβzα⟩=∣zβ,N⋯zβ,1zα,N⋯zα,1⟩\lvert z \rangle = \lvert z_\beta z_\alpha \rangle = \lvert z_{\beta, N} \cdots z_{\beta, 1} z_{\alpha, N} \cdots z_{\alpha, 1} \rangle.

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

H^=∑prσhpr a^pσ†a^rσ+12∑prqsστhprqs a^pσ†a^qτ†a^sτa^rσ,\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},

כאשר ה-hprh_{pr} וה-hprqsh_{prqs} הם מספרים מרוכבים הנקראים אינטגרלים מולקולריים שניתן לחשב מהמפרט של המולקולה באמצעות תוכנית מחשב. במדריך זה, אנו מחשבים את האינטגרלים באמצעות חבילת התוכנה PySCF.

לפרטים על האופן שבו ההמילטוניאן המולקולרי נגזר, התייעץ עם ספר לימוד על כימיה קוונטית (לדוגמה, Modern Quantum Chemistry מאת Szabo ו-Ostlund). להסבר ברמה גבוהה על האופן שבו בעיות כימיה קוונטית ממופות למחשבים קוונטיים, בדוק את ההרצאה Mapping Problems to Qubits מ-Qiskit Global Summer School 2024.

אנסץ local unitary cluster Jastrow (LUCJ)​

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

אנסץ ה-LUCJ הוא צורה מיוחדת של אנסץ ה-unitary cluster Jastrow (UCJ) הכללי, שיש לו את הצורה

∣Ψ⟩=∏μ=1LeK^μeiJ^μe−K^μ∣Φ0⟩ \lvert \Psi \rangle = \prod_{\mu=1}^L e^{\hat{K}_\mu} e^{i \hat{J}_\mu} e^{-\hat{K}_\mu} | \Phi_0 \rangle

כאשר ∣Φ0⟩\lvert \Phi_0 \rangle הוא מצב התייחסות, לעתים קרובות נלקח להיות מצב Hartree-Fock, וה-K^μ\hat{K}_\mu וה-J^μ\hat{J}_\mu יש להם את הצורה

K^μ=∑pq,σKpqμ a^pσ†a^qσ†  ,  J^μ=∑pq,στJpq,στμ n^pσn^qτ  ,\hat{K}_\mu = \sum_{pq, \sigma} K_{pq}^\mu \, \hat{a}^\dagger_{p \sigma} \hat{a}^{\phantom{\dagger}}_{q \sigma} \;,\; \hat{J}_\mu = \sum_{pq, \sigma\tau} J_{pq,\sigma\tau}^\mu \, \hat{n}_{p \sigma} \hat{n}_{q \tau} \;,

כאשר הגדרנו את אופרטור המספר

n^pσ=a^pσ†a^pσ†.\hat{n}_{p \sigma} = \hat{a}^\dagger_{p \sigma} \hat{a}^{\phantom{\dagger}}_{p \sigma}.

האופרטור eK^μe^{\hat{K}_\mu} הוא סיבוב אורביטלי, שניתן ליישם באמצעות אלגוריתמים ידועים בעומק לינארי ותוך שימוש בקישוריות לינארית. יישום איבר ה-eiJke^{i \mathcal{J}_k} של אנסץ ה-UCJ דורש קישוריות מכולם לכולם או שימוש ברשת החלפה פרמיונית, מה שהופך אותו למאתגר עבור מעבדים קוונטיים טרום-סובלניים לשגיאות רועשים שיש להם קישוריות מוגבלת. הרעיון של אנסץ ה-UCJ המקומי הוא להטיל אילוצי דלילות על המטריצות Jαα\mathbf{J}^{\alpha\alpha} ו-Jαβ\mathbf{J}^{\alpha\beta} המאפשרים להן להיות מיושמות בעומק קבוע על טופולוגיות קיוביט עם קישוריות מוגבלת. האילוצים מוגדרים על ידי רשימה של אינדקסים המציינים אילו ערכי מטריצה במשולש העליון מורשים להיות לא-אפס (מאחר שהמטריצות סימטריות, רק המשולש העליון צריך להיות מוגדר). אינדקסים אלה ניתן לפרש כזוגות של אורביטלים המורשים לקיים אינטראקציה.

כדוגמה, שקול טופולוגיית קיוביט של סריג מרובע. אנו יכולים למקם את האורביטלים α\alpha ו-β\beta בקווים מקבילים על הסריג, עם חיבורים בין קווים אלה היוצרים "שלבי סולם" של צורת סולם, כזה:

דיאגרמת מיפוי קיוביטים עבור אנסץ LUCJ על סריג מרובע

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

Jαα:{ (p,p+1)  ,  p=0,…,N−2 }Jαβ:{ (p,p)  ,  p=0,…,N−1 }\begin{align*} \mathbf{J}^{\alpha\alpha} &: \set{(p, p+1) \; , \; p = 0, \ldots, N-2} \\ \mathbf{J}^{\alpha\beta} &: \set{(p, p) \;, \; p = 0, \ldots, N-1} \end{align*}

במילים אחרות, אם מטריצות ה-J\mathbf{J} הן לא-אפס רק באינדקסים המוגדרים במשולש העליון, אז איבר ה-eiJke^{i \mathcal{J}_k} ניתן ליישום על טופולוגיה מרובעת בלי להשתמש בשערי החלפה, בעומק קבוע. כמובן, הטלת אילוצים כאלה על האנסץ הופכת אותו פחות אקספרסיבי, כך שייתכן שיידרשו יותר חזרות אנסץ.

החומרה של IBM® יש לה טופולוגיית קיוביט של סריג heavy-hex, ובמקרה זה אנו יכולים לאמץ תבנית "זיג-זג", המתוארת למטה. בתבנית זו, אורביטלים עם אותו ספין ממופים לקיוביטים עם טופולוגיית קו (עיגולים אדומים וכחולים), וחיבור בין אורביטלים של ספינים שונים קיים בכל אורביטל מרחבי רביעי, כאשר החיבור מתבצע על ידי קיוביט עזר (עיגולים סגולים). במקרה זה, אילוצי האינדקס הם

Jαα:{ (p,p+1)  ,  p=0,…,N−2 }Jαβ:{ (p,p)  ,  p=0,4,8,…(p≤N−1) }\begin{align*} \mathbf{J}^{\alpha\alpha} &: \set{(p, p+1) \; , \; p = 0, \ldots, N-2} \\ \mathbf{J}^{\alpha\beta} &: \set{(p, p) \;, \; p = 0, 4, 8, \ldots (p \leq N-1)} \end{align*}

דיאגרמת מיפוי קיוביטים עבור אנסץ LUCJ על סריג heavy-hex

שחזור תצורה עצמי-עקבי​

תהליך שחזור התצורה העצמי-עקבי מתוכנן לחלץ כמה שיותר אות מדגימות קוונטיות רועשות. מכיוון שההמילטוניאן המולקולרי שומר על מספר החלקיקים ו-spin Z, הגיוני לבחור ansatz מעגלי ששומר גם על הסימטריות הללו. כאשר מוחלים על מצב Hartree-Fock, המצב המתקבל יש לו מספר חלקיקים קבוע ו-spin Z קבוע בהגדרה נטולת רעש. לכן, חצאי ה-spin-α\alpha וה-spin-β\beta של כל מחרוזת ביטים הנדגמת ממצב זה צריכים להיות בעלי אותו משקל Hamming כמו במצב Hartree-Fock. בשל נוכחות הרעש במעבדים קוונטיים נוכחיים, חלק ממחרוזות הביטים הנמדדות יפרו תכונה זו. צורה פשוטה של בחירה לאחר-מעשה תשליך את מחרוזות הביטים הללו, אבל זה בזבזני מכיוון שמחרוזות ביטים אלה עדיין עשויות להכיל אות כלשהו. תהליך השחזור העצמי-עקבי מנסה לשחזר חלק מהאות הזה בעיבוד-לאחר. התהליך איטרטיבי ודורש כקלט אומדן של התפוסות הממוצעות של כל אורביטל במצב היסוד, שמחושב תחילה מהדגימות הגולמיות. התהליך מורץ בלולאה, וכל איטרציה כוללת את השלבים הבאים:

  1. עבור כל מחרוזת ביטים המפרה את הסימטריות המוגדרות, הפוך את הביטים שלה עם תהליך הסתברותי המיועד להביא את מחרוזת הביטים קרוב יותר לאומדן הנוכחי של התפוסות האורביטליות הממוצעות, כדי לקבל מחרוזת ביטים חדשה.
  2. אסוף את כל מחרוזות הביטים הישנות והחדשות שמקיימות את הסימטריות, ודגום תת-קבוצות בגודל קבוע, שנבחר מראש.
  3. עבור כל תת-קבוצה של מחרוזות ביטים, הקרן את ההמילטוניאן לתוך תת-המרחב שנפרש על ידי וקטורי הבסיס המתאימים (ראה את הסעיף הקודם לתיאור וקטורי הבסיס הללו), וחשב אומדן מצב יסוד של ההמילטוניאן המוקרן על מחשב קלאסי.
  4. עדכן את האומדן של התפוסות האורביטליות הממוצעות עם אומדן מצב היסוד עם האנרגיה הנמוכה ביותר.

דיאגרמת זרימת עבודה של SQD​

זרימת העבודה של SQD מתוארת בדיאגרמה הבאה:

דיאגרמת זרימת עבודה של אלגוריתם SQD

דרישות​

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

  • Qiskit SDK v1.0 ומעלה, עם תמיכה בויזואליזציה
  • Qiskit Runtime v0.22 ומעלה (pip install qiskit-ibm-runtime)
  • הרחבת Qiskit של SQD v0.11 ומעלה (pip install qiskit-addon-sqd)
  • ffsim v0.0.58 ומעלה (pip install ffsim)

הגדרות​

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy pyscf qiskit qiskit-addon-sqd qiskit-ibm-runtime rustworkx
import pyscf
import pyscf.cc
import pyscf.mcscf
import ffsim
import numpy as np
import matplotlib.pyplot as plt

from qiskit import QuantumCircuit, QuantumRegister
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

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

במדריך זה, נמצא קירוב למצב היסוד של המולקולה בשיווי משקל בקבוצת הבסיס cc-pVDZ. ראשית, אנו מציינים את המולקולה ואת תכונותיה.

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="cc-pvdz",
    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)

# Store reference energy from SCI calculation performed separately
exact_energy = -109.22690201485733
converged SCF energy = -108.929838385609

לפני בניית מעגל האנסץ LUCJ, אנו מבצעים תחילה חישוב CCSD בתא הקוד הבא. אמפליטודות t1t_1 ו-t2t_2 מחישוב זה ישמשו לאתחול הפרמטרים של האנסץ.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2
E(CCSD) = -109.2177884185543  E_corr = -0.2879500329450045

כעת, אנו משתמשים ב-ffsim כדי ליצור את מעגל האנסץ. מכיוון שלמולקולה שלנו יש מצב Hartree-Fock של קליפה סגורה, אנו משתמשים בוריאנט מאוזן-ספין של אנסץ ה-UCJ, UCJOpSpinBalanced. אנו מעבירים זוגות אינטראקציה המתאימים לטופולוגיית קיוביט של סריג heavy-hex (ראה את סעיף הרקע על אנסץ LUCJ). אנו מגדירים optimize=True בשיטת from_t_amplitudes כדי לאפשר את הפקטוריזציה הכפולה ה"דחוסה" של אמפליטודות ה-t2t_2 (ראה The local unitary cluster Jastrow (LUCJ) ansatz בתיעוד של ffsim לפרטים).

n_reps = 1
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell from running
    # too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

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

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

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)

print(f"Using backend {backend.name}")
Using backend ibm_fez

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

  • בחר קיוביטים פיזיים (initial_layout) מחומרת היעד התואמת לתבנית ה"זיג-זג" שתוארה קודם לכן. פריסת קיוביטים בתבנית זו מובילה למעגל תואם-חומרה יעיל עם פחות שערים. כאן, אנו כוללים קוד לאוטומציה של בחירת תבנית ה"זיג-זג" המשתמשת באבחון הערכה כדי למזער את השגיאות הקשורות לפריסה הנבחרת.
  • ייצר מנהל מעבר מדורג באמצעות פונקציית generate_preset_pass_manager מ-qiskit עם הבחירה שלך של backend ו-initial_layout.
  • הגדר את שלב ה-pre_init של מנהל המעבר המדורג שלך ל-ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT כולל מעברי טרנספיילר של qiskit שמפרקים שערים לסיבובי אורביטלים ואז ממזגים את סיבובי האורביטלים, וכתוצאה מכך פחות שערים במעגל הסופי.
  • הרץ את מנהל המעבר על המעגל שלך.
קוד לבחירה אוטומטית של פריסת "זיג-זג"
from typing import Sequence

import rustworkx
from qiskit.providers import BackendV2
from rustworkx import NoEdgeBetweenNodes, PyGraph

IBM_TWO_Q_GATES = {"cx", "ecr", "cz"}

def create_linear_chains(num_orbitals: int) -> PyGraph:
    """In zig-zag layout, there are two linear chains (with connecting qubits between
    the chains). This function creates those two linear chains: a rustworkx PyGraph
    with two disconnected linear chains. Each chain contains `num_orbitals` number
    of nodes, that is, in the final graph there are `2 * num_orbitals` number of nodes.

    Args:
        num_orbitals (int): Number orbitals or nodes in each linear chain. They are
            also known as alpha-alpha interaction qubits.

    Returns:
        A rustworkx.PyGraph with two disconnected linear chains each with `num_orbitals`
            number of nodes.
    """
    G = rustworkx.PyGraph()

    for n in range(num_orbitals):
        G.add_node(n)

    for n in range(num_orbitals - 1):
        G.add_edge(n, n + 1, None)

    for n in range(num_orbitals, 2 * num_orbitals):
        G.add_node(n)

    for n in range(num_orbitals, 2 * num_orbitals - 1):
        G.add_edge(n, n + 1, None)

    return G

def create_lucj_zigzag_layout(
    num_orbitals: int, backend_coupling_graph: PyGraph
) -> tuple[PyGraph, int]:
    """This function creates the complete zigzag graph that 'can be mapped' to an IBM QPU with
    heavy-hex connectivity (the zigzag must be an isomorphic sub-graph to the QPU/backend
    coupling graph for it to be mapped).
    The zigzag pattern includes both linear chains (alpha-alpha interactions) and connecting
    qubits between the linear chains (alpha-beta interactions).

    Args:
        num_orbitals (int): Number of orbitals, that is, number of nodes in each alpha-alpha linear chain.
        backend_coupling_graph (PyGraph): The coupling graph of the backend on which the LUCJ ansatz
            will be mapped and run. This function takes the coupling graph as a undirected
            `rustworkx.PyGraph` where there is only one 'undirected' edge between two nodes,
            that is, qubits. Usually, the coupling graph of a IBM backend is directed (for example, Eagle devices
            such as ibm_brisbane) or may have two edges between two nodes (for example, Heron `ibm_torino`).
            A user needs to be make such graphs undirected and/or remove duplicate edges to make them
            compatible with this function.

    Returns:
        G_new (PyGraph): The graph with IBM backend compliant zigzag pattern.
        num_alpha_beta_qubits (int): Number of connecting qubits between the linear chains
            in the zigzag pattern. While we want as many connecting (alpha-beta) qubits between
            the linear (alpha-alpha) chains, we cannot accommodate all due to qubit and connectivity
            constraints of backends. This is the maximum number of connecting qubits the zigzag pattern
            can have while being backend compliant (that is, isomorphic to backend coupling graph).
    """
    isomorphic = False
    G = create_linear_chains(num_orbitals=num_orbitals)

    num_iters = num_orbitals
    while not isomorphic:
        G_new = G.copy()
        num_alpha_beta_qubits = 0
        for n in range(num_iters):
            if n % 4 == 0:
                new_node = 2 * num_orbitals + num_alpha_beta_qubits
                G_new.add_node(new_node)
                G_new.add_edge(n, new_node, None)
                G_new.add_edge(new_node, n + num_orbitals, None)
                num_alpha_beta_qubits = num_alpha_beta_qubits + 1
        isomorphic = rustworkx.is_subgraph_isomorphic(
            backend_coupling_graph, G_new
        )
        num_iters -= 1

    return G_new, num_alpha_beta_qubits

def lightweight_layout_error_scoring(
    backend: BackendV2,
    virtual_edges: Sequence[Sequence[int]],
    physical_layouts: Sequence[int],
    two_q_gate_name: str,
) -> list[list[list[int], float]]:
    """Lightweight and heuristic function to score isomorphic layouts. There can be many zigzag patterns,
    each with different set of physical qubits, that can be mapped to a backend. Some of them may
    include less noise qubits and couplings than others. This function computes a simple error score
    for each such layout. It sums up 2Q gate error for all couplings in the zigzag pattern (layout) and
    measurement of errors of physical qubits in the layout to compute the error score.

    Note:
        This lightweight scoring can be refined using concepts such as mapomatic.

    Args:
        backend (BackendV2): A backend.
        virtual_edges (Sequence[Sequence[int]]): Edges in the device compliant zigzag pattern where
            nodes are numbered from 0 to (2 * num_orbitals + num_alpha_beta_qubits).
        physical_layouts (Sequence[int]): All physical layouts of the zigzag pattern that are isomorphic
            to each other and to the larger backend coupling map.
        two_q_gate_name (str): The name of the two-qubit gate of the backend. The name is used for fetching
            two-qubit gate error from backend properties.

    Returns:
        scores (list): A list of lists where each sublist contains two items. First item is the layout, and
            second item is a float representing error score of the layout. The layouts in the `scores` are
            sorted in the ascending order of error score.
    """
    props = backend.properties()
    scores = []
    for layout in physical_layouts:
        total_2q_error = 0
        for edge in virtual_edges:
            physical_edge = (layout[edge[0]], layout[edge[1]])
            try:
                ge = props.gate_error(two_q_gate_name, physical_edge)
            except Exception:
                ge = props.gate_error(two_q_gate_name, physical_edge[::-1])
            total_2q_error += ge
        total_measurement_error = 0
        for qubit in layout:
            meas_error = props.readout_error(qubit)
            total_measurement_error += meas_error
        scores.append([layout, total_2q_error + total_measurement_error])

    return sorted(scores, key=lambda x: x[1])

def _make_backend_cmap_pygraph(backend: BackendV2) -> PyGraph:
    graph = backend.coupling_map.graph
    if not graph.is_symmetric():
        graph.make_symmetric()
    backend_coupling_graph = graph.to_undirected()

    edge_list = backend_coupling_graph.edge_list()
    removed_edge = []
    for edge in edge_list:
        if set(edge) in removed_edge:
            continue
        try:
            backend_coupling_graph.remove_edge(edge[0], edge[1])
            removed_edge.append(set(edge))
        except NoEdgeBetweenNodes:
            pass

    return backend_coupling_graph

def get_zigzag_physical_layout(
    num_orbitals: int, backend: BackendV2, score_layouts: bool = True
) -> tuple[list[int], int]:
    """The main function that generates the zigzag pattern with physical qubits that can be used
    as an `intial_layout` in a preset passmanager/transpiler.

    Args:
        num_orbitals (int): Number of orbitals.
        backend (BackendV2): A backend.
        score_layouts (bool): Optional. If `True`, it uses the `lightweight_layout_error_scoring`
            function to score the isomorphic layouts and returns the layout with less erroneous qubits.
            If `False`, returns the first isomorphic subgraph.

    Returns:
        A tuple of device compliant layout (list[int]) with zigzag pattern and an int representing
            number of alpha-beta-interactions.
    """
    backend_coupling_graph = _make_backend_cmap_pygraph(backend=backend)

    G, num_alpha_beta_qubits = create_lucj_zigzag_layout(
        num_orbitals=num_orbitals,
        backend_coupling_graph=backend_coupling_graph,
    )

    isomorphic_mappings = rustworkx.vf2_mapping(
        backend_coupling_graph, G, subgraph=True
    )
    isomorphic_mappings = list(isomorphic_mappings)

    edges = list(G.edge_list())

    layouts = []
    for mapping in isomorphic_mappings:
        initial_layout = [None] * (2 * num_orbitals + num_alpha_beta_qubits)
        for key, value in mapping.items():
            initial_layout[value] = key
        layouts.append(initial_layout)

    two_q_gate_name = IBM_TWO_Q_GATES.intersection(
        backend.configuration().basis_gates
    ).pop()

    if score_layouts:
        scores = lightweight_layout_error_scoring(
            backend=backend,
            virtual_edges=edges,
            physical_layouts=layouts,
            two_q_gate_name=two_q_gate_name,
        )

        return scores[0][0][:-num_alpha_beta_qubits], num_alpha_beta_qubits

    return layouts[0][:-num_alpha_beta_qubits], num_alpha_beta_qubits
initial_layout, _ = get_zigzag_physical_layout(num_orbitals, backend=backend)

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")
Gate counts (w/o pre-init passes): OrderedDict({'rz': 12438, 'sx': 12169, 'cz': 3986, 'x': 573, 'measure': 52, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'sx': 7059, 'rz': 6962, 'cz': 1876, 'measure': 52, 'x': 35, 'barrier': 1})

שלב 3: הרצה באמצעות Primitives של Qiskit​

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

sampler = Sampler(mode=backend)
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

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

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

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
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)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    pub_result.data.meas,
    samples_per_batch=samples_per_batch,
    norb=num_orbitals,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=12345,
)
Iteration 1
	Subsample 0
		Energy: -108.97781410104506
		Subspace dimension: 28561
	Subsample 1
		Energy: -108.97781410104506
		Subspace dimension: 28561
	Subsample 2
		Energy: -108.97781410104506
		Subspace dimension: 28561
Iteration 2
	Subsample 0
		Energy: -109.05150860754739
		Subspace dimension: 287296
	Subsample 1
		Energy: -109.08152283263908
		Subspace dimension: 264196
	Subsample 2
		Energy: -109.11636893067873
		Subspace dimension: 284089
Iteration 3
	Subsample 0
		Energy: -109.15878555367885
		Subspace dimension: 429025
	Subsample 1
		Energy: -109.16462831275786
		Subspace dimension: 473344
	Subsample 2
		Energy: -109.15895026995382
		Subspace dimension: 435600
Iteration 4
	Subsample 0
		Energy: -109.17784051223317
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.1774651326829
		Subspace dimension: 657721
	Subsample 2
		Energy: -109.18085212360191
		Subspace dimension: 617796
Iteration 5
	Subsample 0
		Energy: -109.18636242518915
		Subspace dimension: 815409
	Subsample 1
		Energy: -109.18451014767518
		Subspace dimension: 837225
	Subsample 2
		Energy: -109.18333728638888
		Subspace dimension: 857476

הצגה חזותית של התוצאות​

התרשים הראשון מראה שאחרי כמה איטרציות אנו מעריכים את אנרגיית מצב היסוד בטווח של ~41 mH (דיוק כימי מקובל בדרך כלל להיות 1 kcal/mol ≈\approx 1.6 mH). ניתן לשפר את האנרגיה על ידי אפשור יותר איטרציות של שחזור תצורה או הגדלת מספר הדגימות לכל batch.

התרשים השני מראה את התפוסה הממוצעת של כל אורביטל מרחבי לאחר האיטרציה הסופית. אנו יכולים לראות שגם האלקטרונים של spin-up וגם של spin-down תופסים את חמשת האורביטלים הראשונים עם הסתברות גבוהה בפתרונות שלנו.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

פלט של תא הקוד הקודם

סקר המדריך​

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

קישור לסקר