יצירת מעגלי התפתחות בזמן של קרילוב (SKQD)

חשוב

המושגים במדריך זה זמינים כרגע רק ב-API של Python. פונקציונליות מקבילה תהיה זמינה ב-API של C בגרסה עתידית.

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

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

שמתקבלים מהתפתחות של מצב ייחוס $\lvert \psi_0 \rangle$ למכפלות גדלות של צעד זמן קבוע $\Delta t$. מחרוזות ביטים שנדגמות מ-$D$ המעגלים האלה ממלאות את תת-מרחב קרילוב, שבו ההמילטוניאן מאולכסן לאחר מכן באופן קלאסי.

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

הערה

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

1. בניית ההמילטוניאן של SIAM כאופרטור פרמיוני

מודל אנדרסון עם מוקש בודד ממקם אורביטל מוקש אחד שמקיים אינטראקציה (אינדקס $0$) במגע עם אמבט בלתי מקיים אינטראקציה של $L$ אתרים המסודרים כשרשרת. ההמילטוניאן שלו מתפצל לחלק חד-גופי $H_1$ ולחלק דו-גופי $H_2$,

$$H = H_1 + H_2,$$

כאשר החלק החד-גופי אוסף את אנרגיית האתר של המוקש (הפוטנציאל הכימי $\mu$), את ההיברידיזציה בין המוקש לאמבט $V$, ואת קפיצת האמבט בין שכנים קרובים $t$,

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

והחלק הדו-גופי נושא את הדחייה הקולומבית באתר $U$ על המוקש בלבד,

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

מדריך זה עובד בנקודה הסימטרית חלקיק-חור $\mu = -U/2$. שני החלקים האלה ממופים ישירות לאינטגרלים החד-גופיים h1e (קפיצת האמבט, היברידיזציית המוקש, אנרגיית האתר של המוקש) ולאינטגרל הדו-גופי היחיד h2e (איבר הקולון של המוקש). קודם כול מרכיבים אותם בבסיס המיקום, שבו המודל מוגדר באופן טבעי.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. מעבר לבסיס התנע

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

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

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

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

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

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

המקומיות הזו היא מה שמשאיר את האינטראקציה זולה לסינתזה בהמשך.

מרכיבים את האינטגרלים החד-גופיים ל-FermionOperator באמצעות בנאי האינטגרלים האלקטרוניים from_1body_tril_spin_sym(), שמצפה שהם יהיו ארוזים בסדר משולש תחתון. החלק הדו-גופי הוא איבר יחיד של מספר-מספר, ולכן במקום להעביר אותו דרך מנגנון האינטגרלים הדו-גופיים, בונים אותו ישירות עם cre()/ann(): אופרטור הקולון באתר $U \, n_{p\uparrow} n_{p\downarrow}$ עבור מוד המוקש $p$ (מוד $\uparrow$ שלו באינדקס p, מוד $\downarrow$ שלו באינדקס p + norb בפריסת בלוקי ספין). שני הבנאים מכפילים בספין את norb האורביטלים המרחביים ל-2 * norb מודים פרמיוניים.

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

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

שינוי הבסיס הוא טרנספורמציית דמיון יוניטרית, ולכן הוא משאיר את הספקטרום ללא שינוי. אשרו שהאופרטור נכון על ידי השוואת הערך העצמי הנמוך ביותר שלו בסקטור $(\text{norb}, \text{nelec})$ מול אלכסון מדויק. ה-FermionOperator חושף LinearOperator של SciPy (שנתמך על ידי גרעין מטריצה-וקטור מקורי של אינטראקציית קונפיגורציות מלאה (FCI)), ולכן אפשר להעביר אותו ישירות ל-scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

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

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

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

3. הכנת מצב הייחוס

מצב הייחוס של SKQD $\lvert \psi_0 \rangle$ הוא סופרפוזיציה של כל העירורים של האלקטרונים הקרובים ביותר לרמת פרמי אל מודי התנע הריקים הקרובים ביותר. הפרסום של SKQD בונה אותו עם רשת שנכתבה ביד של XXPlusYYGates. הוא נראה בעל קורלציה חזקה (יש לו מאות משרעות שאינן אפס בבסיס החישובי), אבל רשת של סיבובים דו-מודיים ששומרים על מספר החלקיקים היא פעולה חד-חלקיקית (גאוסיאנית פרמיונית), ולכן המצב שהיא מפיקה הוא דטרמיננטת סלייטר יחידה. כלומר אפשר להכין אותו עם שער PrepareSlaterDeterminant אחד, בהינתן הסיבוב החד-חלקיקי שהרשת מממשת.

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

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

הערה

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

הכנת דטרמיננטת סלייטר נעשית לכל סקטור ספין בנפרד, כי לכל סקטור יש צורה משלו (אלקטרונים, אורביטלים). שער PrepareSlaterDeterminant אחד (איכלוס, שבו nocc המודים הנמוכים ביותר מלאים, יחד עם הסיבוב החד-חלקיקי) מכין סקטור אחד; אותו שער מוחל פעם אחת על מודי האלפא 0..norb ופעם אחת על מודי הבטא norb..2*norb, כי שני סקטורי הספין חולקים כאן את אותו איכלוס ייחוס ואותו סיבוב:

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

השער ההצהרתי הזה מחליף את רשת XXPlusYYGate שנכתבה ביד; ה-Transpiler מייצר את סינתזת סיבובי גיבנס.

4. הרכבת מעגלי קרילוב

כל מעגל קרילוב הוא הכנת מצב הייחוס ואחריה התפתחות תחת $e^{-i k \, \Delta t \, H}$. מפצלים את ההמילטוניאן לחלק החד-גופי ולחלק הדו-גופי, $H = H_1 + H_2$, ומבצעים טרוטריזציה של ההתפתחות לפי הפיצול הזה. זו הבחירה המכריעה לעומק המעגל. ההתפתחות החד-גופית $e^{-i t H_1}$ היא פעולה של פרמיונים חופשיים (גאוסיאנית פרמיונית), ולכן היא סיבוב בסיס חד-חלקיקי (OrbitalRotation עם יוניטרית $e^{-i t \, \mathbf{h}_1}$) ולא משהו שצריך לבצע לו טרוטריזציה פאולית איבר אחר איבר. ההתפתחות הדו-גופית $e^{-i t H_2}$ פועלת רק על איבר האינטראקציה היחיד באתר. בניית הצעד כך מאפשרת ל-Transpiler לפלוט מבנה לבנים רדוד של סיבובי גיבנס ועוד שער אינטראקציה דו-קיוביטי יחיד, המבנה הרדוד שSKQD בונה ביד בפרסום שלו.

כבר יש לנו את החלק הדו-גופי. זה האופרטור two_body_momentum שנבנה בשלב 2 (האיבר היחיד באתר). ההתפתחות החד-גופית $e^{-i t H_1}$ אינה תלויה בספין, ולכן היא OrbitalRotation יחיד לכל סקטור עם יוניטרית $e^{-i t \, \mathbf{h}_1}$ שמוחל פעם אחת על כל סקטור ספין, באותו מבנה מקביל כמו בהכנת מצב הייחוס.

כל מעגל קרילוב מתפתח לזמן כולל $k \, \Delta t$. זה ממומש באמצעות מכפלת טרוטר מסדר שני של $k$ צעדים בגודל $\Delta t$, כך ששגיאת הצעד נשארת קבועה כשמימד קרילוב גדל. כל צעד מכניס סיבוב חד-גופי מלא (אחד לכל סקטור ספין) בין שני חצאי צעדים של ההתפתחות הדו-גופית (הזולה, בעלת האיבר היחיד). מימד קרילוב $D$ הוא מספר המעגלים כאלה, כש-dim נע מ-0 (מצב ייחוס בלבד) עד $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

הערה

המעגלים שנבנו כאן עדיין אינם כוללים מדידות, כי אותם מעגלים פרמיוניים מניעים גם את הסימולציה המדויקת בשלב 6, ומדידה אינה פעולה יוניטרית. שלב 5 מוסיף אותן באמצעות FermionicCircuit.measure_all(), שקורא את מוד i לתוך הביט הקלאסי i.

עדיף למדוד את המעגל הפרמיוני ולא את המעגל שעבר הידור. מדידה פרמיונית נקבת בשם המוד שהיא קוראת, ולכן הצבה וניתוב נושאים אותה יחד עם ה-Qubit שהמוד הגיע אליו, ומחרוזת הביטים שנדגמת מאונדקסת לפי מוד. מדידת המעגל שעבר הידור משאירה לכם לבטל בעצמכם את final_index_layout(), ומחזירה בשקט תוצאות מעורבבות אם שוכחים. זה נשען על קידוד ג'ורדן-ויגנר בבסיס האיכלוס שבו משתמשת ההגדרה המוכנה, שבו כל מוד נישא על ידי Qubit אחד.

5. הידור למעגלי Qubit

הידור של FermionicCircuit ממפה את השערים הפרמיוניים שלו לשערי Qubit. זה משתמש במיפוי ג'ורדן-ויגנר jordan_wigner(). שערי OrbitalRotation עוברים סינתזה למבנה לבנים של סיבובי גיבנס, וה-Evolution בעל האיבר היחיד מורד לשער אינטראקציה דו-קיוביטי יחיד.

מכיוון שלכל שער פרמיוני כאן יש סינתזת ברירת מחדל, אפשר להשתמש ישירות ב-generate_preset_jw_pass_manager() המוכן במקום להרכיב ידנית את השלבים. כל ארגומנט מילת מפתח מועבר לשלב ה-Qubit; העבירו optimization_level=0 לתמונה נאמנה ובלתי מותאמת של העומק המסונתז. יצירת משפחת מעגלי קרילוב המלאה היא אז לולאה על מימד קרילוב. שמרו את ה-FermionicCircuits שלא עברו הידור לצד מעגלי ה-Qubit שעברו הידור. הפרמיוניים מניעים את הסימולציה המדויקת (וקטור מצב) בשלב 6, ולכן הם חייבים להישאר ללא מדידות; measure_all(inplace=False) מחזיר עותק עם מדידות להידור ומשאיר את המקור ללא שינוי, והעותק הזה הוא מה ש-Backend היה מריץ.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

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

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

הערה

פיצול ההמילטוניאן שומר על העומקים האלה צנועים. ההתפתחות החד-גופית היא סיבוב אורביטלים אמיתי, ולכן היא עוברת סינתזה למבנה לבנים של XXPlusYYGate על Qubits סמוכים. כל חצי צעד של האינטראקציה באתר הוא RZZGate יחיד שמצמד את שני מודי הספין של המוקש; האיבר הזול שהמכפלה מסדר שני מכניסה סביב הסיבוב. העברת ההמילטוניאן כולו ל-Evolution אחד הייתה מבצעת במקום זאת טרוטריזציה פאולית של החלק החד-גופי איבר אחר איבר, מפיצה מחרוזות Z ארוכות של ג'ורדן-ויגנר על פני צימודי המוקש-אמבט (ארוכי הטווח) ומנפחת את העומק בסדר גודל.

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

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. דגימת מחרוזות ביטים ממעגלי קרילוב

בחומרה, הייתם מריצים את המעגלים שעברו הידור וממדדים. זו הגרסה נטולת הרעש (וקטור מצב) של השלב הזה, שמראה את מחרוזות הביטים ש-SKQD היה אוסף. מדמים כל FermionicCircuit שלא עבר הידור עם ffsim.apply_unitary() ודוגמים את וקטור המצב המתקבל עם ffsim.sample_state_vector().

שער PrepareSlaterDeterminant פועל בסימולציה בשיטת אימות-ואז-סיבוב (הוא בודק שהמצב הנכנס תופס את דטרמיננטת הייחוס שלו, ואז מפעיל את הסיבוב). לכן הקלט הוא דטרמיננטת האיכלוס הפשוטה: nocc המודים הנמוכים ביותר מלאים לכל ספין, בלי סיבוב. השער מפעיל את reference_rotation בעצמו.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

כשמאגדים את הדגימות מכל חמשת מעגלי קרילוב, התומך קטן מאוד: כמה מאות מחרוזות ביטים שונות מתוך $\binom{8}{4}^2 = 4900$ הדטרמיננטות בסקטור (4, 4). הריכוז הזה, תוצאה ישירה של דלילות בסיס התנע שהוכחה בשלב 2, הופך את האלכסון הקלאסי הבא לבר-ביצוע.

הצגת הספירות באמצעות qiskit.visualization.plot_histogram() הופכת את הדלילות לגלויה. כמה קונפיגורציות שולטות, ובראשן הדטרמיננטה של הייחוס.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. אלכסון במרחב המשנה שנדגם

מחרוזות הביטים שנדגמו הן הקלט לחצי הקלאסי של SKQD. ההמילטוניאן מוקרן על מרחב המשנה שנפרש על ידי הקונפיגורציות שנדגמו ומאולכסן שם. זה מועבר לחבילה qiskit-addon-sqd, שהפונקציה diagonalize_fermionic_hamiltonian() שלה מריצה את לולאת האלכסון הקוונטי מבוסס הדגימה המלאה. היא בונה מרחב משנה מאצוות של הקונפיגורציות שנדגמו, מאלכסנת בו את ההמילטוניאן, ומשפרת את מרחב המשנה באופן איטרטיבי באמצעות שחזור קונפיגורציות: היפוך של תפוסות כדי לתקן קונפיגורציות שמפרות את סימטריית מספר החלקיקים הידועה. שחזור קונפיגורציות נועד לבטל את הנזק של רעש החומרה, שהופך מחרוזות ביטים שנדגמו לסקטור שגוי של מספר חלקיקים; הדגימה כאן נטולת רעש, ולכן אין מה לתקן מבחינה זו.

הפותר צורך את הספירות שנדגמו בתור BitArray, שנבנה ישירות מה-counts המאוחדים של שלב 6. העבירו את num_bits במפורש: בלעדיו, הרוחב נקבע לפי מחרוזת הביטים הרחבה ביותר שנדגמה, והוא צר מ-2 * norb בכל פעם שאף דגימה לא תופסת את המצב האחרון, והפותר אז דוחה את המערך. הפותר מאלכסן באותו בסיס תנע שבו משתמש שאר הבנייה, ולכן הוא מקבל ישירות את האינטגרלים h1e_momentum ו-h2e_momentum משלב 2, בלי אריזה מחדש. האנרגיה ששוחזרה תואמת את הייחוס המדויק בפער קטן בהרבה ממילי-הארטרי, מכמה מאות קונפיגורציות שנדגמו בלבד. מכיוון שאנרגיית SQD היא חסם עליון וריאציוני לאנרגיית מצב היסוד האמיתית, היא תמיד נמצאת ב-reference_energy או מעליו (שעדיין בתוקף משלב 2):

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

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

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

בכך נסגרת לולאת SKQD. המעגלים בבסיס התנע בשלבים 1–5 מפיקים דגימה דלילה (שלב 6), והאלכסון מבוסס הדגימה שלמעלה הופך את הדגימה הזו בחזרה לאנרגיית מצב יסוד.

השלבים הבאים

במדריך הזה הרצנו את כל צינור SKQD על סימולטור וקטור מצב נטול רעש. בחומרה, ה-circuits שעברו טרנספילציה היו מורצים ונמדדים במקום דגימת וקטור המצב בשלב 6, והספירות שהתקבלו היו מוזנות לאלכסון של שלב 7 ללא שינוי. עיינו בתיעוד של Qiskit לעזרה בהרצת מעגלים, ובמדריכי ה-SQD addon למידע על העיבוד שלאחר האלכסון של מרחב המשנה, כולל התנהגותו על דגימות רועשות, שם שחזור הקונפיגורציות עושה את רוב העבודה.

קראו את מדריך הקשר ל-ffsim כדי להבין מדוע ffsim.apply_unitary() ו-ffsim.sample_state_vector() עובדים באופן טבעי על המעגלים הפרמיוניים של החבילה הזו, ומדוע סקטור מספר החלקיקים הקבוע שבו נעשה שימוש בשלב 6 הופך את הדלילות בבסיס התנע, שנקבעה בשלב 2, ליעילה מבחינת דגימות.