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

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

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

الخلفية​

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

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

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

نموذج المجموعة العنقودية الوحدوية المحلية Jastrow (LUCJ)​

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

نموذج LUCJ هو شكل متخصص من النموذج العام للمجموعة العنقودية الوحدوية 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 هي الحالة المرجعية، وغالبًا ما تكون حالة هارتري-فوك، و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

استعادة التكوين ذاتي التناسق​

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

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

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

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

ناتج خلية الكود السابقة

استطلاع البرنامج التعليمي​

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

رابط الاستطلاع