SQD لتقدير طاقة هاميلتوني كيميائي

في هذا الدرس، سنطبّق SQD لتقدير طاقة الحالة الأساسية لجزيء.

تحديدًا، سنناقش المواضيع التالية باستخدام نهج نمط Qiskit المكوّن من $4$ خطوات:

1. الخطوة 1: تعيين المسألة إلى دوائر كمية وعوامل
   • إعداد الهاميلتوني الجزيئي لـ $N_2$.
   • شرح نموذج Jastrow للعنقود الأحادي الموحّد المحلي (LUCJ) المستوحى من الكيمياء والمتوافق مع العتاد [1]
2. الخطوة 2: تحسين الأداء للعتاد المستهدف
   • تحسين عدد البوابات وتخطيط النموذج الأولي للتنفيذ على العتاد
3. الخطوة 3: التنفيذ على العتاد المستهدف
   • تشغيل الدائرة المحسَّنة على وحدة معالجة كمية حقيقية (QPU) لتوليد عينات الفضاء الجزئي.
4. الخطوة 4: معالجة النتائج
   • تقديم حلقة الاسترداد الذاتي للتكوينات [2]
     • معالجة المجموعة الكاملة من عينات السلاسل الثنائية، مع الاستفادة من المعرفة المسبقة بعدد الجسيمات ومتوسط إشغال المدارات المحسوب في آخر تكرار.
     • إنشاء دُفعات احتمالية من العينات الفرعية المستعادة من السلاسل الثنائية.
     • إسقاط الهاميلتوني الجزيئي وتقطيره على كل فضاء جزئي من العينات.
     • حفظ أدنى طاقة للحالة الأساسية عبر جميع الدُفعات وتحديث متوسط إشغال المدارات.

سنستخدم عدة حزم برمجية خلال هذا الدرس.

• PySCF لتعريف الجزيء وإعداد الهاميلتوني.
• حزمة ffsim لبناء النموذج الأولي LUCJ.
• Qiskit لتهيئة النموذج الأولي للتنفيذ على العتاد.
• Qiskit IBM Runtime لتنفيذ الدائرة على وحدة QPU وجمع العينات.
• إضافة Qiskit addon SQD لاسترداد التكوينات وتقدير طاقة الحالة الأساسية باستخدام إسقاط الفضاء الجزئي وتقطير المصفوفات.

1. تعيين المسألة إلى دوائر كمية وعوامل

الهاميلتوني الجزيئي

يأخذ الهاميلتوني الجزيئي الصورة العامة التالية:

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$ هي عوامل إنشاء/إفناء الفيرميونات المرتبطة بعنصر مجموعة الأساس من الرتبة $p$ والدوران $\sigma$. أما $h_{pr}$ و$(pr|qs)$ فهما تكاملات الإلكترون أحادية وثنائية الجسم. باستخدام pySCF، سنعرّف الجزيء ونحسب التكاملات أحادية وثنائية الجسم للهاميلتوني بمجموعة الأساس 6-31g.

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy pyscf qiskit qiskit-addon-sqd qiskit-ibm-runtime

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# 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)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

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

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

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

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

في هذا الدرس، سنستخدم تحويل جوردان-ويغنر (JW) لتعيين دالة موجة الفيرميونات إلى دالة موجة كيوبتية حتى يمكن إعدادها باستخدام دائرة كمية. يُعيّن تحويل JW فضاء فوك للفيرميونات في M مدار مكاني إلى فضاء هيلبرت لـ 2M كيوبت، أي أن المدار المكاني يُقسَّم إلى مدارَي دوران، أحدهما مرتبط بإلكترون دوران للأعلى ($\alpha$) والآخر بدوران للأسفل ($\beta$). يمكن أن يكون مدار الدوران محتلًا أو شاغرًا. عادةً، حين نشير إلى عدد المدارات، نعني عدد المدارات المكانية. أما عدد مدارات الدوران فيكون ضعف ذلك. في الدوائر الكمية، سنمثّل كل مدار دوران بكيوبت واحد. وبذلك، ستمثّل مجموعة من الكيوبتات المدارات ذات الدوران للأعلى أو $\alpha$، وتمثّل مجموعة أخرى المدارات ذات الدوران للأسفل أو $\beta$. على سبيل المثال، جزيء $N_2$ بمجموعة الأساس 6-31g يمتلك $16$ مدارًا مكانيًا (أي $16$ مدار $\alpha$ + $16$ مدار $\beta$ = $32$ مدار دوران). وبالتالي، نحتاج إلى دائرة كمية من $32$ كيوبت (قد نحتاج إلى كيوبتات مساعدة إضافية كما سيُناقش لاحقًا). تُقاس الكيوبتات في الأساس الحسابي لتوليد سلاسل ثنائية تمثّل التكوينات الإلكترونية أو المحددات (Slater). سنستخدم طوال هذا الدرس مصطلحات السلاسل الثنائية والتكوينات والمحددات بشكل متبادل. تُخبرنا السلاسل الثنائية عن إشغال الإلكترونات في مدارات الدوران: الرقم $1$ في موضع بت ما يعني أن مدار الدوران المقابل مشغول، بينما $0$ يعني أن المدار شاغر. ولأن مسائل البنية الإلكترونية تحافظ على عدد الجسيمات، يجب أن يكون عدد محدد من مدارات الدوران مشغولًا. يمتلك جزيء $N_2$ $5$ إلكترونات ذات دوران للأعلى ($\alpha$) و$5$ إلكترونات ذات دوران للأسفل ($\beta$). وبالتالي، يجب أن تحتوي أي سلسلة ثنائية تمثّل مدارات $\alpha$ و$\beta$ على خمسة أرقام $1$ لكل منها في جزيء $N_2$.

1.1 الدائرة الكمية لتوليد العينات: النموذج الأولي LUCJ

في هذا الدرس، سنستخدم النموذج الأولي Jastrow للعنقود الأحادي الموحّد المحلي (LUCJ) \[1\] لإعداد الحالة الكمية وأخذ العينات منها لاحقًا. أولًا، سنشرح مختلف المكوّنات الأساسية لنموذج UCJ الكامل والتقريبات المُجراة في نسخته المحلية. بعد ذلك، باستخدام حزمة ffsim، سنبني النموذج الأولي LUCJ ونحسّنه باستخدام محوّل Qiskit للتنفيذ على العتاد.

يأخذ النموذج الأولي UCJ الشكل التالي (لحاصل ضرب $L$ طبقة أو تكرار من عامل UCJ):

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

حيث $\vert \Phi_{0} \rangle$ هي الحالة المرجعية، وعادةً تُؤخذ كحالة هارتري-فوك (HF). بما أن حالة هارتري-فوك تُعرَّف بأن تكون المدارات الأقل رقمًا هي المشغولة، فإن إعداد حالة HF سيتضمّن تطبيق بوابات X لضبط الكيوبتات المقابلة للمدارات المشغولة على القيمة واحد. على سبيل المثال، قد يبدو كتلة إعداد حالة HF لـ 4 مدارات مكانية و2 إلكترون ذو دوران للأعلى و2 ذو دوران للأسفل كالتالي: يتكوّن تكرار واحد من عامل UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ من تطور كولوم قطري ($e^{iJ^{(\mu)}}$) محاط بدورانات مدارية ($e^{K^{(\mu)}}$ و$e^{-K^{(\mu)}}$). تعمل كتل الدوران المداري على نوع واحد من الدوران ($\alpha$ (للأعلى)/$\beta$ (للأسفل)). لكل نوع إلكتروني، يتكوّن الدوران المداري من طبقة من بوابات $R_{z}$ أحادية الكيوبت تليها سلسلة من بوابات دوران Givens ثنائية الكيوبت (بوابات $XX + YY$).

تعمل بوابات الكيوبتَين على مدارات الدوران المتجاورة (الكيوبتات الأقرب جوارًا)، وبالتالي يمكن تنفيذها على وحدات IBM® QPU دون الحاجة إلى بوابات SWAP. يتكوّن $e^{iJ^{(\mu)}}$، المعروف أيضًا بعامل كولوم القطري، من ثلاثة كتل. اثنتان منها تعملان على قطاعات الدوران نفسها ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ و$e^{iJ_{\beta \beta}^{(\mu)}}$)، وواحدة تعمل بين قطاعَي الدوران ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

تتكوّن جميع الكتل في $e^{iJ^{(\mu)}}$ من بوابات عدد-عدد $U_{nn}(\phi)$ [1]. يمكن تحليل بوابة $U_{nn}(\phi)$ بشكل أعمق إلى بوابة $R_{ZZ}(\frac{\phi}{2})$ تليها بوابتا $Rz(-\frac{\phi}{2})$ أحاديتا الكيوبت تعملان على كيوبتَين منفصلَين.

تحتوي المكوّنات ذات الدوران المتماثل ($J_{\alpha \alpha}$ و$J_{\beta \beta}$) على بوابات $U_{nn}$ بين جميع أزواج الكيوبتات المحتملة. ومع ذلك، نظرًا لأن وحدات QPU فائقة التوصيل ذات اتصال محدود، يجب تبديل الكيوبتات لتحقيق البوابات بين الكيوبتات غير المتجاورة.

على سبيل المثال، تأمّل الكتلة $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (أو $e^{iJ_{\beta \beta}^{(\mu)}}$) التالية لـ $N = 4$ مدارات مكانية. بالنسبة لاتصال خطي بين الكيوبتات، لا يمكن تنفيذ البوابات الثلاثة الأخيرة مباشرةً لأنها تعمل بين كيوبتات غير متجاورة (مثلًا، Q0 وQ2 غير متصلَين مباشرةً). لذلك نحتاج إلى بوابات SWAP لجعلها متجاورة (يُظهر الشكل التالي مثالًا بـ $3$ بوابات SWAP). بعد ذلك، تُنفّذ $J_{\alpha \beta}$ بوابات بين المدارات ذات نفس الفهرس من قطاعات الدوران المختلفة (مثلًا، بين $0\alpha$ و$0\beta$). وبالمثل، إذا لم تكن الكيوبتات متجاورة فيزيائيًا على وحدة QPU، فستحتاج هذه البوابات أيضًا إلى SWAPs. من النقاش السابق، يواجه النموذج الأولي UCJ بعض العقبات للتنفيذ على العتاد لأنه يحتاج إلى بوابات SWAP بسبب التفاعلات بين الكيوبتات غير المتجاورة. يعالج النموذج المحلي للنموذج الأولي UCJ، وهو LUCJ، هذا التحدي عن طريق إزالة بعض $U_{nn}$ من عامل كولوم القطري.

في كتل نفس النوع الإلكتروني ($J_{\alpha \alpha}$ و$J_{\beta \beta}$)، نحتفظ فقط ببوابات $U_{nn}$ المتوافقة مع الاتصال بين أقرب الجيران ونحذف البوابات بين الكيوبتات غير المتجاورة في نسخة LUCJ. يُظهر الشكل التالي كتلة LUCJ بعد حذف البوابات غير المتجاورة. بعد ذلك، يمكن أن تأخذ نسخة LUCJ من كتلة $J_{\alpha \beta}$ التي تعمل بين أنواع إلكترونية مختلفة أشكالًا متعددة بناءً على طوبولوجيا الجهاز.

هنا أيضًا، تتخلص نسخة LUCJ من البوابات غير المتوافقة. يُظهر الشكل أدناه متغيرات كتلة $J_{\alpha \beta}$ لطوبولوجيا كيوبت مختلفة تشمل الشبكة المربعة والسداسية وهيفي-هكس والخطية.

• مربعة: يمكن وضع بوابات $U_{nn}$ بين جميع مدارات $\alpha$ و$\beta$ دون أي SWAPs، وبالتالي لا نحتاج إلى حذف أي بوابات $U_{nn}$.
• هيفي-هكس: تُحفظ تفاعلات $\alpha$-$\beta$ بين المدارات ذات الفهرس المضاعف للـ $4$-الأول (مثل الفهرس 0 و4 و8) وتكون مُوسَّطة بكيوبتات مساعدة، أي نحتاج إلى كيوبتات مساعدة بين السلاسل الخطية التي تمثّل مدارات $\alpha$ و$\beta$. يحتاج هذا الترتيب إلى عدد محدود من SWAPs.
• سداسية: كل مدار آخر، مثل المدارات ذات الفهرس 0 و2 و4، يصبح أقرب جار عند ترتيب $\alpha$ و$\beta$ في سلسلتَين خطيتَين متجاورتَين.
• خطية: يكون مدار $\alpha$ واحد فقط ومدار $\beta$ واحد متصلَين، مما يعني أن كتلة $J_{\alpha \beta}$ ستحتوي على بوابة واحدة فقط. بينما يجعل حذف البوابات من نموذج UCJ لبناء نسخة LUCJ التنفيذ على العتاد أكثر ملاءمة، إلا أن النموذج الأولي يفقد بعض القدرة التعبيرية. وبالتالي، قد تكون هناك حاجة إلى تكرارات ($L$) أكثر من عامل UCJ المعدَّل عند استخدام نموذج LUCJ.

1.2 تهيئة النموذج الأولي LUCJ

LUCJ نموذج أولي مُعامَل، ونحتاج إلى تهيئة المعاملات قبل التنفيذ على العتاد. إحدى طرق تهيئة النموذج الأولي هي استخدام سَعات t1 وt2 من طريقة العنقود الأحادي المزدوج الكلاسيكية (CCSD)، حيث سَعات t1 هي معاملات عوامل الإثارة المفردة وسَعات t2 هي لعوامل الإثارة المزدوجة.

لاحظ أنه رغم أن تهيئة نموذج LUCJ بسَعات t1 وt2 تُعطي نتائج جيدة، إلا أن معاملات النموذج الأولي قد تحتاج إلى مزيد من التحسين.

# 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]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 بناء النموذج الأولي LUCJ باستخدام ffsim

سنستخدم حزمة ffsim لإنشاء النموذج الأولي وتهيئته بسَعات t1 وt2 المحسوبة أعلاه. بما أن جزيئنا يمتلك حالة هارتري-فوك ذات غلاف مغلق، سنستخدم النسخة متوازنة الدوران من نموذج UCJ، وهي UCJOpSpinBalanced.

نظرًا لأن عتاد IBM يمتلك طوبولوجيا هيفي-هكس، سنعتمد نمط zig-zag المستخدم في [1] والموضَّح أعلاه لتفاعلات الكيوبتات. في هذا النمط، المدارات (الكيوبتات) ذات الدوران نفسه متصلة بطوبولوجيا خطية (الدوائر الحمراء والزرقاء). بسبب طوبولوجيا هيفي-هكس، تمتلك المدارات ذات دوران مختلف اتصالات بين كل مدار من الرتبة الرابعة، أي المدارات ذات الفهرس 0 و4 و8 وما إلى ذلك (الدوائر البنفسجية).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
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),
)

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()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

يمكن تحسين نموذج LUCJ ذي الطبقات المتكررة عن طريق دمج بعض الكتل المتجاورة. تأمّل حالة n_reps=2. يمكن دمج كتلتَي الدوران المداري في الوسط في كتلة دوران مداري واحدة. تمتلك حزمة ffsim مدير تمريرات بالاسم ffsim.qiskit.PRE_INIT لتحسين الدائرة عن طريق دمج هذه الكتل المتجاورة.

2. تحسين الأداء للعتاد المستهدف

أولًا، نُحضر الواجهة الخلفية التي اخترناها. سنحسّن دائرتنا للواجهة الخلفية، ثم ننفّذ الدائرة المحسَّنة على نفس الواجهة الخلفية لتوليد عينات الفضاء الجزئي.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

بعد ذلك، نوصي باتباع الخطوات التالية لتحسين النموذج الأولي وجعله متوافقًا مع العتاد.

• اختر كيوبتات فيزيائية (initial_layout) من العتاد المستهدف تتوافق مع نمط zig-zag (سلسلتَين خطيتَين مع كيوبت مساعد بينهما) الموضَّح أعلاه. يؤدي ترتيب الكيوبتات وفق هذا النمط إلى دائرة متوافقة مع العتاد بكفاءة وعدد بوابات أقل.
• أنشئ مدير تمريرات مُرحَّل باستخدام دالة generate_preset_pass_manager من Qiskit مع backend وinitial_layout اللذَين اخترتهما.
• اضبط مرحلة pre_init في مدير التمريرات المُرحَّل على ffsim.qiskit.PRE_INIT. يتضمّن ffsim.qiskit.PRE_INIT تمريرات محوّل Qiskit التي تحلّل البوابات إلى دورانات مدارية ثم تدمج هذه الدورانات، مما يؤدي إلى عدد بوابات أقل في الدائرة النهائية.
• شغّل مدير التمريرات على دائرتك.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

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': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. التنفيذ على العتاد المستهدف

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

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. المعالجة اللاحقة للنتائج

يمكن تلخيص جزء المعالجة اللاحقة في سير عمل SQD باستخدام المخطط التالي.

أخذ عينات من دالة LUCJ الأنساتز في الأساس الحسابي يولّد مجموعة من التهيئات الضوضائية $\tilde{\mathcal{\chi}}$، والتي تُستخدم في روتين المعالجة اللاحقة. يتضمن هذا الروتين طريقةً تُسمى استرداد التهيئة (سيأتي شرحها لاحقاً) لتصحيح التهيئات ذات الأعداد الخاطئة من الإلكترونات باحتمالية معينة. ثم تُأخذ عينات فرعية من التهيئات التي تمتلك أعداداً صحيحة من الإلكترونات $\tilde{\mathcal{\chi}}_{R}$ فقط، وتُوزَّع على دُفعات متعددة بناءً على تكرار ظهور كل تهيئة فريدة. كل دفعة من العينات تُعرّف فضاءً جزئياً ($\mathcal{S^{(k)}}$). بعد ذلك، يُسقَط هاميلتوني الجزيء $H$ على الفضاءات الجزئية:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

ثم يُغذَّى كل هاميلتوني مُسقَط $H_{\mathcal{S}^{(k)}}$ إلى حلّال القيم الذاتية، حيث يُقطَّر لحساب القيم الذاتية والمتجهات الذاتية لإعادة بناء الحالة الذاتية. في هذا الدرس، نُسقط الهاميلتوني ونُقطّره باستخدام حزمة qiskit-addon-sqd، التي تستخدم طريقة دافيدسون من PySCF للتقطير.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

نجمع بعد ذلك أدنى قيمة ذاتية (طاقة) من الدُفعات، ونحسب أيضاً متوسط إشغال المدارات $\text{n}$. تُستخدم معلومات متوسط الإشغال في خطوة استرداد التهيئة لتصحيح التهيئات الضوضائية باحتمالية معينة.

نشرح فيما يلي حلقة استرداد التهيئة ذاتية الاتساق بالتفصيل، ونعرض أمثلة برمجية ملموسة لتطبيق الخطوات المذكورة أعلاه لتقدير طاقة الحالة الأساسية لهاميلتوني $N_2$.

4.1 استرداد التهيئة: نظرة عامة

كل بت في السلسلة الثنائية (محدد سلاتر) يمثل مداراً دورانياً. النصف الأيمن من السلسلة يمثل المدارات ذات الدوران الصاعد، والنصف الأيسر يمثل المدارات ذات الدوران النازل. القيمة 1 تعني أن المدار مشغول بإلكترون، والقيمة 0 تعني أن المدار فارغ. نعرف مسبقاً العدد الصحيح للجسيمات (إلكترونات الدوران الصاعد والنازل على حدٍّ سواء). لنفترض أن لدينا محدداً $x$ يحتوي على $N_x$ إلكتروناً (أي أن هناك $N_x$ من القيم 1 في السلسلة). العدد الصحيح للجسيمات هو $N$. إذا كان $N_x \neq N$، فنعلم أن السلسلة قد أُفسدت بالضوضاء. يحاول روتين التهيئة ذاتية الاتساق تصحيح السلسلة بقلب $|N_x - N|$ بت باحتمالية معينة، مستعيناً بمعلومات متوسط إشغال المدارات. يخبرنا متوسط الإشغال ($n$) بمدى احتمال أن يكون المدار مشغولاً بإلكترون. إذا كان $N_x < N$، فهذا يعني أن لدينا إلكترونات أقل مما ينبغي، ونحتاج إلى قلب بعض الأصفار إلى آحاد، والعكس بالعكس.

احتمال القلب يمكن أن يكون $|x[i] - avg\_occupancy[i]|$ للمدار الدوراني i. في [2]، استخدم المؤلفون احتمالاً موزوناً للقلب باستخدام دالة ReLU المعدَّلة.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

هنا $h$ يحدد موضع "الزاوية" في دالة ReLU، والمعامل $\delta$ يحدد قيمة دالة ReLU عند الزاوية. عند $\delta = 0$، تصبح $w$ دالة ReLU الحقيقية، وعند $\delta > 0$، تصبح ReLU المعدَّلة. في الورقة البحثية، استخدم المؤلفون $\delta = 0.01$ و$h =$ عدد جسيمات ألفا (أو بيتا) / عدد مدارات ألفا (أو بيتا) الدورانية $= N/M$ (عامل الملء).

متوسط إشغال المدارات ($n$) غير معروف مسبقاً. تبدأ الدورة الأولى لتقدير الحالة الأساسية بالتهيئات التي تمتلك فقط أعداداً صحيحة من الجسيمات في كلا النوعين الدورانيين. بعد الدورة الأولى، يكون لدينا تقدير أولي للحالة الأساسية، ونستطيع باستخدام هذا التقدير بناء التخمين الأول لـ$n$. يُستخدم هذا التخمين لـ$n$ لاسترداد التهيئات، وتشغيل دورة التقدير التالية للحالة الأساسية، وتحسين التخمين لـ$n$ بشكل ذاتي الاتساق. تتكرر العملية حتى يُستوفى معيار التوقف.

خذ المثال التالي لـ$N = 2$ و$x = |1000\rangle$ ($N_x = 1$). نحتاج إلى قلب أحد الأصفار إلى واحد لتصحيح أعداد الجسيمات، والخيارات هي 1100 و1010 و1001. بناءً على احتمال القلب، سيُختار أحد الخيارات بوصفه التهيئة المُستردة (أي السلسلة ذات العدد الصحيح من الجسيمات).

لنفترض أننا نُشغّل دُفعتين في الدورة الأولى، والحالتان الأساسيتان المقدَّرتان منهما هما:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

باستخدام حالات الأساس الحسابي وسعاتها، يمكننا حساب احتمال إشغال الإلكترونات (باختصار الإشغالات) لكل مدار دوراني (كيوبت) (علماً بأن الاحتمال = |السعة|$^2$). في الجداول أدناه نُدوّن الإشغالات لكل كيوبت لكل سلسلة تظهر في الحالة الأساسية المقدَّرة، ونحسب إجمالي إشغال المدارات لكل دُفعة. لاحظ أنه وفق اتفاقية ترتيب Qiskit، يمثل البت الأيمن الكيوبت 0 (Q0)، ويمثل البت الأيسر Q3.

الإشغال (الدُفعة 0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

الإشغال (الدُفعة 1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

الإشغال (المتوسط عبر الدُفعات)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (average)	0.49	0.51	0.35	0.65

باستخدام متوسط إشغال المدارات المحسوب أعلاه، يمكننا إيجاد احتمالات القلب للمدارات المختلفة في التهيئة $x = \vert 1000 \rangle$. بما أن المدار الذي يمثله Q3 مشغول بالفعل ولا يحتاج إلى قلب، نضع p(flip) له يساوي $0$. أما بقية المدارات غير المشغولة، فاحتمال قلب كل منها هو $\vert x[i] - \text{n}[i] \vert$. إلى جانب p(flip)، نحسب أيضاً الوزن الاحتمالي المرتبط بالقلب باستخدام دالة ReLU المعدَّلة الموضحة أعلاه.

احتمال القلب ($x = \vert 1000 \rangle$، $\delta = 0.01$، $h = N/M = 2/4 = 0.50$)

	Q3	Q2	Q1	Q0
p(flip) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(flip))	0	0.03	0.007	0.31

أخيراً، باستخدام الاحتمالات الموزونة أعلاه، يمكننا قلب أحد المدارات غير المشغولة Q2 أو Q1 أو Q0. بناءً على القيم أعلاه، Q0 هو الأرجح للقلب، والتهيئة المحتملة المُستردة هي $\vert \text{1001} \rangle$. يمكن تلخيص عملية استرداد التهيئة ذاتية الاتساق الكاملة على النحو التالي:

الدورة الأولى: لنفترض أن السلاسل الثنائية (التهيئات أو محددات سلاتر) التي يولّدها الحاسوب الكمومي تشكّل مجموعة $\widetilde{\chi}$، والتي تضم التهيئات ذات الأعداد الصحيحة ($\widetilde{\chi}_{correct}$) والأعداد الخاطئة ($\widetilde{\chi}_{incorrect}$) من الجسيمات في كل قطاع دوراني.

1. تُأخذ عينات عشوائية من التهيئات في ($\widetilde{\chi}_{correct}$) لإنشاء دُفعات $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ من المتجهات للإسقاط على الفضاء الجزئي. عدد الدُفعات وعدد العينات في كل دُفعة معاملان يحددهما المستخدم. كلما زاد عدد العينات في كل دُفعة، زادت أبعاد الفضاء الجزئي وأصبح التقطير أكثر تطلباً من الناحية الحسابية. في المقابل، قد يُفضي عدد العينات الصغير جداً إلى إغفال متجهات دعم الحالة الأساسية والتوصل إلى تقديرات خاطئة.
2. تشغيل حلّال الحالة الذاتية (أي الإسقاط على الفضاء الجزئي والتقطير) على الدُفعات والحصول على الحالات الذاتية التقريبية $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. بناء التخمين الأول لـ$n$ من الحالات الذاتية التقريبية.

الدورات اللاحقة:

1. استخدام $n$ لتصحيح التهيئات ذات الأعداد الخاطئة من الجسيمات في $\widetilde{\chi}_{incorrect}$. لنسمِّها $\widetilde{\chi}_{correct\_new}$. عندئذٍ، $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ تشكّل المجموعة الجديدة من التهيئات ذات الأعداد الصحيحة من الجسيمات.
2. أخذ عينات من $\widetilde{\chi}_{R}$ لإنشاء دُفعات $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. تشغيل حلّال الحالة الذاتية مع الدُفعات الجديدة وتوليد تقديرات جديدة للحالات الأساسية $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. بناء التخمين المحسَّن لـ$n$ من الحالات الذاتية التقريبية.
5. إذا لم يُستوفَ معيار التوقف، العودة إلى الخطوة 2.1.

4.2 تقدير الحالة الأساسية

أولاً، سنحوّل العدّات إلى مصفوفة سلاسل ثنائية ومصفوفة احتمالات للمعالجة اللاحقة.

كل صف في المصفوفة يمثل سلسلة ثنائية فريدة واحدة. بما أن الكيوبتات مُرقَّمة من اليمين في السلسلة الثنائية في Qiskit، يمثل العمود 0 الكيوبت N-1، ويمثل العمود N-1 الكيوبت 0، حيث N هو عدد الكيوبتات.

المدارات ألفا ممثَّلة في نطاق فهرس الأعمدة (N, N/2] (النصف الأيمن)، والمدارات بيتا ممثَّلة في النطاق (N/2, 0] (النصف الأيسر).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

ثمة خيارات عديدة يتحكم فيها المستخدم وهي مهمة لهذه التقنية:

• iterations: عدد دورات استرداد التهيئة ذاتية الاتساق
• n_batches: عدد دُفعات التهيئات المستخدمة في الاستدعاءات المختلفة لحلّال الحالة الذاتية
• samples_per_batch: عدد التهيئات الفريدة المُدرجة في كل دُفعة
• max_davidson_cycles: الحد الأقصى لعدد دورات دافيدسون التي يُشغّلها كل حلّال قيم ذاتية

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We post-select here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_davidson=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 مناقشة النتائج

يُظهر الرسم الأول أننا بعد بضع دورات نقدّر طاقة الحالة الأساسية بفارق ~24 mH (الدقة الكيميائية مقبولة عموماً عند 1 kcal/mol $\approx$ 1.6 mH). يُظهر الرسم الثاني متوسط إشغال كل مدار فراغي بعد الدورة الأخيرة. يتضح أن إلكترونات الدوران الصاعد والنازل على حدٍّ سواء تشغل المدارات الخمسة الأولى باحتمالية عالية في حلولنا.

على الرغم من أن طاقة الحالة الأساسية المقدَّرة معقولة، إلا أنها لا تقع ضمن حد الدقة الكيميائية ($\pm \approx 1.6$ mH). يمكن إرجاع هذا الفارق إلى صغر أبعاد الفضاء الجزئي المستخدم أعلاه للإسقاط والتقطير. بما أننا استخدمنا samples_per_batch=500، فإن الفضاء الجزئي يمتد على أقصى تقدير على $500$ متجه، وهو ما يُغفل متجهات دعم الحالة الأساسية. زيادة قيمة المعامل samples_per_batch ينبغي أن تُحسّن الدقة على حساب موارد الحوسبة الكلاسيكية ووقت التشغيل.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

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

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

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-6)
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})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

تمرين للقارئ

زِد قيمة المعامل samples_per_batch تدريجياً (مثلاً، من $1000$ إلى $10000$ بخطوة $1000$؛ بحسب ما تسمح به ذاكرة حاسوبك) وقارن طاقات الحالة الأساسية المقدَّرة.

المراجع

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.