انتقل إلى المحتوى الرئيسي

التقطير الكمي القائم على العينات لهاميلتوني كيميائي

تقدير الاستخدام: أقل من دقيقة واحدة على معالج Heron r2 (ملاحظة: هذا تقدير فقط. قد يتفاوت وقت التشغيل الفعلي.)

الخلفية النظرية​

في هذا البرنامج التعليمي، نوضح كيفية المعالجة اللاحقة للعينات الكمية الضوضائية للحصول على تقريب لحالة الأساس لجزيئة النيتروجين N2\text{N}_2 عند طول الرابطة التوازنية، باستخدام خوارزمية التقطير الكمي القائم على العينات (SQD)، لعينات مأخوذة من دائرة كمية مكونة من 59 كيوبت (52 كيوبت نظام + 7 كيوبت مساعد). تتوفر تطبيقات مبنية على Qiskit في إضافات SQD لـ Qiskit، ويمكن الاطلاع على مزيد من التفاصيل في التوثيق المقابل مع مثال بسيط للبدء.

SQD هي تقنية لإيجاد القيم الذاتية والمتجهات الذاتية للمؤثرات الكمية، كالهاميلتوني لنظام كمي، باستخدام الحوسبة الكمية والحوسبة الكلاسيكية الموزعة معًا. تُستخدم الحوسبة الكلاسيكية الموزعة لمعالجة العينات المستخرجة من معالج كمي، وإسقاط هاميلتوني مستهدف وتقطيره قطريًا في الفضاء الجزئي الذي تمتد عليه تلك العينات. يجعل هذا النهجَ SQD مقاومًا للعينات المتأثرة بالضوضاء الكمية، وقادرًا على التعامل مع الهاميلتونيات الكبيرة، كهاميلتونيات الكيمياء التي تحتوي على ملايين حدود التفاعل، والتي تتجاوز نطاق أي طرق تقطير قطري دقيق. يتضمن سير عمل SQD الخطوات التالية:

  1. اختيار نموذج الدائرة (circuit ansatz) وتطبيقه على حاسوب كمي انطلاقًا من حالة مرجعية (في هذه الحالة، حالة هارتري-فوك).
  2. أخذ عينات من السلاسل الثنائية (bitstrings) من الحالة الكمية الناتجة.
  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.

نموذج التجمع الأحادي الموضعي (LUCJ)​

يستلزم SQD وجود نموذج دائرة كمية لأخذ العينات منه. في هذا البرنامج التعليمي، سنستخدم نموذج التجمع الأحادي الموضعي (LUCJ) نظرًا لجمعه بين الدوافع الفيزيائية وملاءمته للأجهزة الكمية.

نموذج LUCJ هو شكل متخصص من نموذج التجمع الأحادي العام (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 هي الحالة المرجعية، وغالبًا ما تكون حالة هارتري-فوك، و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 إما ترابطًا كليًا بين الكيوبتات أو استخدام شبكة مبادلة فيرميونية، مما يجعله تحديًا لمعالجات كمية قبل تحقيق الخطأ المتسامح به وذات الترابط المحدود. الفكرة الجوهرية لنموذج LUCJ هي فرض قيود تفرق على مصفوفتي 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)، ويمكن في هذه الحالة اعتماد نمط "متعرج" (zig-zag) كما هو موضح أدناه. في هذا النمط، تُربط المدارات ذات السبين المتماثل بكيوبتات ذات هيكل خطي (الدوائر الحمراء والزرقاء)، وتتوفر وصلة بين مدارات السبين المختلف عند كل رابع مدار مكاني، حيث تُتيح هذه الوصلة كيوبت مساعد (الدوائر البنفسجية). في هذه الحالة، قيود الفهارس هي:

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 على شبكة السداسي الثقيل

استرداد التهيئات بصورة ذاتية الاتساق​

صُمِّم إجراء استرداد التهيئات بصورة ذاتية الاتساق لاستخلاص أكبر قدر ممكن من الإشارة من العينات الكمية الضوضائية. نظرًا لأن الهاميلتوني الجزيئي يحفظ عدد الجسيمات والسبين Z، فمن المنطقي اختيار نموذج دائرة يحفظ هذين التماثلين أيضًا. عند تطبيق النموذج على حالة هارتري-فوك، تكون الحالة الناتجة ذات عدد جسيمات وسبين Z ثابتين في الحالة الخالية من الضوضاء. لذلك، يجب أن يكون لنصفَي السبين-α\alpha والسبين-β\beta في أي سلسلة ثنائية مأخوذة من هذه الحالة نفس الوزن الهاميلي كما في حالة هارتري-فوك. بسبب وجود الضوضاء في المعالجات الكمية الحالية، ستنتهك بعض السلاسل الثنائية المقيسة هذه الخاصية. أبسط صورة للانتقاء اللاحق هي استبعاد هذه السلاسل الثنائية، غير أن ذلك مُضيِّع لأن تلك السلاسل قد تحتوي بعض الإشارة. يسعى إجراء الاسترداد ذاتي الاتساق إلى استعادة بعض تلك الإشارة في مرحلة ما بعد المعالجة. الإجراء تكراري ويتطلب كمدخلات تقديرًا لمتوسط الإشغال لكل مدار في حالة الأساس، والذي يُحسب أولًا من العينات الخام. يُشغَّل الإجراء في حلقة تكرارية، وتتضمن كل دورة الخطوات التالية:

  1. لكل سلسلة ثنائية تنتهك التماثلات المحددة، يتم قلب بتاتها باستخدام إجراء احتمالي مصمم لتقريب السلسلة الثنائية من التقدير الحالي لمتوسط إشغال المدارات، للحصول على سلسلة ثنائية جديدة.
  2. تجميع كل السلاسل الثنائية القديمة والجديدة التي تحقق التماثلات، وأخذ عينات فرعية من مجموعات بحجم ثابت محدد مسبقًا.
  3. لكل مجموعة فرعية من السلاسل الثنائية، يُسقَط الهاميلتوني في الفضاء الجزئي الذي تمتد عليه المتجهات الأساسية المقابلة (انظر القسم السابق لوصف هذه المتجهات الأساسية)، وتُحسب تقدير حالة الأساس للهاميلتوني المسقوط على حاسوب كلاسيكي.
  4. تحديث تقدير متوسط إشغال المدارات بتقدير حالة الأساس ذي الطاقة الأدنى.

مخطط سير عمل SQD​

يُوضَّح سير عمل SQD في المخطط التالي:

مخطط سير عمل خوارزمية SQD

المتطلبات​

قبل البدء بهذا البرنامج التعليمي، تأكد من تثبيت ما يلي:

  • Qiskit SDK الإصدار 1.0 أو أحدث، مع دعم الرسوميات
  • Qiskit Runtime الإصدار 0.22 أو أحدث (pip install qiskit-ibm-runtime)
  • إضافة SQD لـ Qiskit الإصدار 0.11 أو أحدث (pip install qiskit-addon-sqd)
  • ffsim الإصدار 0.0.58 أو أحدث (pip install ffsim)

الإعداد​

# Added by doQumentation — installs packages not in the Binder environment
%pip install -q ffsim qiskit-addon-sqd
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

الخطوة الأولى: تحويل المدخلات الكلاسيكية إلى مسألة كمية​

في هذا البرنامج التعليمي، سنجد تقريبًا لحالة الأساس للجزيئة عند حالة التوازن في مجموعة الأساس 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 لإنشاء دائرة النموذج. بما أن جزيئتنا تمتلك حالة هارتري-فوك ذات قشرة مغلقة، نستخدم المتغير المتوازن للسبين من نموذج UCJ، وهو UCJOpSpinBalanced. نمرر أزواج التفاعل المناسبة لهيكل كيوبت شبكة السداسي الثقيل (انظر قسم الخلفية النظرية حول نموذج LUCJ). نضبط optimize=True في دالة from_t_amplitudes لتفعيل التحليل المضغوط المزدوج لسعات t2t_2 (انظر نموذج التجمع الأحادي الموضعي (LUCJ) في توثيق 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()

الخطوة الثانية: تحسين المسألة لتنفيذها على الأجهزة الكمية​

بعد ذلك، نُحسِّن الدائرة لجهاز كمي مستهدف.

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) من الجهاز المستهدف وفق النمط المتعرج (zig-zag) الموصوف سابقًا. يؤدي ترتيب الكيوبتات بهذا النمط إلى دائرة متوافقة مع الجهاز وذات بوابات أقل. نُدرج هنا بعض التعليمات البرمجية لأتمتة اختيار نمط "المتعرج" باستخدام تقييم إرشادي يُقلل أخطاء التخطيط المختارة.
  • إنشاء مدير ممرات متعدد المراحل باستخدام دالة 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)

```python
            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: التنفيذ باستخدام Qiskit primitives​

بعد تحسين الدائرة لتنفيذها على العتاد، أصبحنا جاهزين لتشغيلها على العتاد المستهدف وجمع العينات لتقدير طاقة الحالة الأساسية. بما أن لدينا دائرة واحدة فقط، سنستخدم وضع تنفيذ Job في Qiskit Runtime وننفّذ دائرتنا.

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

الخطوة 4: المعالجة اللاحقة وإرجاع النتيجة بالتنسيق الكلاسيكي المطلوب​

الآن، نقدّر طاقة الحالة الأساسية للـ Hamiltonian باستخدام الدالة 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). يمكن تحسين الطاقة بالسماح بمزيد من تكرارات استرداد التهيئة أو بزيادة عدد العينات لكل دُفعة.

يُظهر الرسم البياني الثاني متوسط إشغال كل مدار مكاني بعد التكرار الأخير. يتضح لنا أن الإلكترونات ذات الدوران الصاعد وذات الدوران الهابط تحتل المدارات الخمسة الأولى باحتمالية عالية في حلولنا.

# 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()

Output of the previous code cell

استبيان الدرس التعليمي​

يُرجى ملء هذا الاستبيان القصير لتقديم ملاحظاتك حول هذا الدرس التعليمي. ستساعدنا آراؤك في تحسين محتوانا وتجربة المستخدم.