خوارزمية SqDRIFT لتقدير الحالة الأساسية

Python

تنزيل دفتر Python

نزّل دفتر Python لتشغيل هذا الدرس التعليمي بنفسك. يقوم الدفتر بتنزيل مدخل FCIDump إذا لم يكن موجودًا بالفعل.

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

مخرجات التعلم

• تعلّم كيفية إنشاء دارات بعمق أصغر مقارنة بـ Trotterization

• اطّلع على سير عمل متكامل من البداية إلى النهاية لتقدير الحالة الأساسية باستخدام qDRIFT وSQD

• تعلّم كيفية استخدام qiskit-fermions جنبًا إلى جنب مع إضافات Qiskit الأخرى لتنفيذ سير عمل كهذا

المتطلبات الأساسية

• اقرأ نظرة عامة على التقطير الكمومي المعتمد على العينات (SQD)

• اقرأ درس التقطير الكمومي بطريقة Krylov المعتمدة على العينات (SKQD)

الخلفية

SqDRIFT هي نسخة متغيرة من SKQD تستبدل الحاجة إلى اختيار أنساتز (ansatz) تُؤخذ منه عينات من سلاسل البتات بمجموعة من دارات التطور الزمني المبنية مباشرة من الهاميلتوني المستهدف. يتحقق ذلك بأخذ عينات فرعية من مؤثرات تطور زمني أصغر من الهاميلتوني استنادًا إلى معاملاته، وهو ما يُعرف بطريقة qDRIFT للتروتّرة (Trotterization).

يستخدم كلا التنفيذين Qiskit Fermions لتمثيل الهاملتوني الجزيئي وبناء الدارات لخوارزمية qDRIFT. يستخدم تنفيذ Python مسارات المحوّل (transpiler) الفرميونية، بينما يأخذ تنفيذ C++ عينات المجموعات صراحةً ويربطها بمؤثرات الكيوبت.

ليكن الهاميلتوني بالصيغة التالية:

$$H = \sum_{i=1}^{N} c_i h_i$$

حيث، دون فقدان للعمومية، نشترط $c_i > 0$ وأن تكون القيمة الذاتية الكبرى لـ $h_i$ مساوية، بالقيمة المطلقة، لـ $1$. يُستوعَب أي معامل إشاري أو معقد داخل $h_i$، لذا فإن المعاملات $c_i$ هي أوزان موجبة تمامًا بينما تحمل $h_i$ اتجاه كل حد. هنا $N$ هو عدد الحدود (أو، بعد التجميع، عدد المجموعات) في الهاميلتوني؛ وهو خاصية للهاميلتوني ويختلف عن عدد المؤثرات المأخوذة بالعينة في دارة واحدة، والمكتوب $n$ أدناه.

ثم تُحقق خوارزمية qDRIFT، للوقت المستهدف $t$، مؤثرًا ما $V_k$، حيث يتراوح $k$ من $1 \cdots K$ ويدل على دارة SqDRIFT رقم $k$، معرّفًا كالتالي:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

هنا $n$ هو عدد المؤثرات المأخوذة بالعينة لكل دارة و$K$ هو عدد الدارات في المجموعة. يمتد الجداء عبر السحوبات الـ $n$، وليس عبر جميع حدود الهاميلتوني الـ $N$، ونظرًا لأن الحدود تُسحب مع الإرجاع، يمكن أن يظهر نفس $h_i$ أكثر من مرة في $V_k$ واحد.

الكمية:

$$\lambda = \sum_{i=1}^{N} c_i$$

هي معيار $L_1$ للمعاملات، لذا تتطور كل خطوة من الخطوات الـ $n$ لنفس المدة $\lambda t / n$ بغض النظر عن الحد الذي تم سحبه. تجانس زاوية الخطوة هو السمة المميزة لـ qDRIFT: يؤثر المعامل في النتيجة من خلال عدد مرات سحب حده، لا من خلال مقدار دوران ذلك الحد. تُؤخذ الفهارس بالعينة من التوزيع:

$$P[k_i] = \frac{c_i}{\lambda}$$

لذا فإن المتسلسلة $(k_1, \ldots, k_n)$ هي سلسلة عشوائية من فهارس الحدود مسحوبة من هذا التوزيع. وبما أن $c_i$ موجبة ومجموعها يساوي $\lambda$، فإن هذا توزيع احتمالي مُطبَّع، ويقرِّب توقع القناة الناتجة عبر السحوبات العشوائية التطور تحت $H$، بخطأ يتناقص مع تزايد $n$. لاحظ أن خطأ التقريب يعتمد على $\lambda$ لا على عدد الحدود $N$.

(تكتب ورقة SqDRIFT عدد الحدود بصيغة $\mathcal{N}$ وطول المتسلسلة بصيغة $N$؛ نستخدم هنا $N$ و$n$ للإبقاء على تمييز واضح بينهما.)

يوضح هذا الدرس التعليمي كيفية توليد مجموعة من هذه الدارات العشوائية. وبعد إنشاء هذه الدارات، على غرار كيفية إنشاء فضاء فرعي من نوع Krylov لمؤثرات مختلفة، نأخذ عينات من سلاسل البتات من عدة مؤثرات كهذه بمعاملات زمنية مختلفة. يضمن هذا تداخلاً أعلى بين متجهات الحالة الأساسية وسلاسل البتات المأخوذة بالعينة.

المتطلبات

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

• بيئة Python افتراضية (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (لاحظ أن الاسم بصيغة الجمع)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

يمكنك تثبيت جميع الحزم المطلوبة باستخدام:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

الإعداد

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

التنفيذ

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

تستخدم خطوات Python أدناه محاكي Aer. ويأتي مثال العتاد بعد ذلك.

قراءة وتحضير FCIDump

في هذا الدرس التعليمي، سنحمّل هاميلتوني البنية الإلكترونية للنيتروجين (N2). توجد طرق أخرى أيضًا لإنشاء مؤثرات فرميونية. راجع التوثيق في qiskit_fermions.operators.library.

عن ملف FCIDump هذا. يصف الملف N2_sto_3g جزيء نيتروجين ($N_2$) بالأساس الأدنى STO-3G عند مسافة بين ذرية قدرها 1.09 $\AA$، وهي طول رابطة التوازن التجريبي. يعلن ترويسته NORB=10 وNELEC=14 وMS2=0: 10 مدارات فضائية (وبالتالي 20 مدارًا مغزليًا، و20 كيوبت تحت تحويل جوردان-فيغنر)، و14 إلكترونًا في حالة مفردة مغزلية، أي سبعة إلكترونات $\alpha$ وسبعة إلكترونات $\beta$. تُعطى جميع المدارات رمز تناظر 1، أي أنه لا يُستغل أي تناظر زمرة نقطية. وبما أن هذا ملف STO-3G للفضاء الكامل، فلا تُجمَّد أي مدارات، وفضاء الترابط صغير بما يكفي لحساب طاقة مرجعية FCI دقيقة كلاسيكيًا للمقارنة، كما هو موضح في كتلة التعليمات البرمجية التالية.

يمكن إعادة توليد ملف مكافئ باستخدام PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

بما أن التكاملات تعتمد على مدارات SCF المتقاربة، فقد يختلف الملف المُعاد توليده عن الملف المرفق في طور المدار أو ترتيبه؛ لا تتأثر الطاقات الكلية بذلك.

الحصول على الملف. ابحث عن ملف FCIDump في مستودع GitHub هذا. يمكنك تشغيل كتلة التعليمات البرمجية أدناه لجلبه إلى الموقع الذي يتوقعه بقية الدرس.

أولًا نستخدم cisolver المقدَّم من pyscf للحصول على الطاقة المرجعية. هذه هي طاقة الحالة الأساسية الحقيقية للجزيء الذي نعمل عليه. لهذا سنعلن أولًا norb وnelec، وهما عدد المدارات وعدد الإلكترونات على التوالي. ثم نعلن h1e وh2e، وهما التكاملات أحادية وثنائية الإلكترون على التوالي. كل هذه ستُستخدم لاحقًا في SQD أيضًا.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

تحميل الهاملتوني

مع جهوزية البيانات اللازمة، نقرأ الهاملتوني من ملف FCI بصيغة متوافقة مع qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

سير عمل فرميوني باستخدام qiskit-fermions

سنقوم أولًا بتحويل الهاملتوني إلى نموذج دارة فرميونية باستخدام qiskit-fermions، الذي يوفر مسارات ترانسبايلر وبوابات خاصة بالدارات الفرميونية. ستُستخدم هذه لاحقًا قبل مسارات ترانسبايلر Qiskit التقليدية لسير العمل هذا.

تجميع الحدود

لضمان قابلية إعادة إنتاج النتائج، نستخدم أولًا canonical_order لترتيب الحدود بناءً على بنيتها فقط. لذلك يكون ترتيب المؤثرات في قائمة canon ثابتًا. يضمن هذا قابلية إعادة إنتاج المؤثرات المُنشأة لأن مسار QDriftTrotterization الذي سنستخدمه لاحقًا يأخذ عينات من فهارس عشوائية لإنشاء مؤثرات qDRIFT.

في هذه الخطوة، نستغل التناظرات الكثيرة الموجودة في هاملتوني البنية الإلكترونية عن طريق تجميع الحدود المترابطة ذات المعاملات المتطابقة. ومع أن القيام بذلك يغيّر توزيع معاملات المؤثر الذي يأخذ منه بروتوكول qDRIFT عيناته، فإن هذا لا يؤثر على ضمانات تقاربه. والأهم أن تجميع الحدود المرتبطة بالتناظر يؤدي إلى إلغاء مواتٍ لحدود باولي وإلى عمق دارة أقل إجمالًا عند تطوير حالة زمنيًا تحت تأثيرها.

يوفر qiskit-fermions الدالة group_terms_by_electronic_structure التي تقوم بهذا التجميع نيابة عنا.

لاحظ أن group_terms_by_electronic_structure تفترض حدودًا مرتبة ترتيبًا طبيعيًا.

تصفية الحدود القطرية

نزيل الحدود القطرية من الهاميلتوني المستخدم لتوليد الدارات، حتى تُخصَّص فتحات أخذ عينات qDRIFT الـ$n$ للحدود التي تنقل الإشغال بين التهيئات. من الأفضل تصفية هذه الحدود من الهاميلتوني في هذه المرحلة، قبل بناء بوابة Evolution في الخطوة التالية.

الحدود المقصودة هي تلك القطرية في أساس عدد الإشغال، أي نواتج مؤثرات العدد $a^\dagger_i a_i$. تندرج ثلاثة أنواع من الحدود تحت هذا الوصف:

• إزاحة الطاقة الثابتة، وهي ناتج صفر مؤثرات عدد، حيث يساهم تطورها الزمني بطور شامل فقط؛

• مؤثرات العدد الفردية $n_i$، التي يختزل تطورها الزمني إلى دورانات $Z$ أحادية الكيوبت؛

• النواتج من الرتبة الأعلى مثل $n_i n_j$.

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

تزيل الدالة filter_diagonal_terms() هذه الحدود من المؤثر في مكانه. تحددها من بنيتها المرتبة ترتيبًا طبيعيًا — التوافق بين مجموعة أنماط الإنشاء ومجموعة أنماط الإفناء — لذا فهي صالحة فقط على مؤثر مرتّب بالفعل ترتيبًا طبيعيًا. لا يُتحقق من هذا الافتراض وقت التشغيل.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

الآن بعد أن جمّعنا الحدود في الهاملتوني، سنحدد المعاملات التالية لتوليد مجموعة الدارات:

• عدد الدارات المراد توليدها: num_circuits
• طول كل دارة من حيث مجموعات الإثارة: num_exc
• عامل أزمنة التطور المختلفة: times

إنشاء دارات فرميونية

سننشئ الآن دارات فرميونية لكل خطوة زمنية. تتكون كل دارة من بوابة تطور واحدة، بزمن التطور الذي حددناه سابقًا. مؤثر التطور هو الهاملتوني. لاحقًا سنشغّل مسارات الترانسبايلر على هذه الدارات لإنشاء دارات qDRIFT.

تحضير الأنساتز

نحضّر حالة هارتري-فوك باستخدام الصنف InitializeModes. بالنسبة للنيتروجين، تتمثل العملية ببساطة في تطبيق بوابات X على أول num_elec_a كيوبت ثم على num_elec_b كيوبت، وكلاهما يساوي سبعة بالنسبة للنيتروجين. تمثل هذه الحالة إلكترونات النيتروجين السبعة $\alpha$ والسبعة $\beta$.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

الخطوة 2: تحسين المسألة لتنفيذها على عتاد الحوسبة الكمومية

الآن بعد أن أصبحت لدينا دارّاتنا، سنستخدم أولًا المسارات المتوفرة في qiskit-fermions لإجراء تحسينات على المستوى الفرميوني، ثم تحويل (transpile) دارّاتنا للـ Backend المختار. بما أن هذه تجربة محاكاة، فسنقوم بذلك أولًا لأجل AerSimulator.

حساب الوزن لكل مجموعة

في هذه الخطوة، نجري أخذ عينات qDRIFT للحدود بشكل عشوائي مع احتمالات تتناسب مع معاملاتها في الهاملتوني. يقوم مسار ترانسبايلر qDRIFT بهذا نيابة عنا. يمكننا الآن إنشاء دارات أقل عمقًا يمكن تنفيذها على العتاد بكفاءة أكبر رغم محدودية اتصال الكيوبتات، حتى عندما يحتوي الهاملتوني على اقترانات بعيدة المدى وحدود أعلى من التربيعية. بعد تجميع الحدود، يأخذ عينات من المؤثرات بناءً على أوزانها. لكل مؤثر $h_i$، يُعرَّف الوزن $W_{h_i}$ كما يلي:

$$W_{h_i} = |c_i| / \lambda$$

بما أن الحدود جُمّعت في الخطوة 1، فإن كل $h_i$ هنا هو مجموعة كاملة: $c_i$ هو متوسط القيمة المطلقة لمعاملات الحدود في المجموعة $i$، ويتطور كل حد في المجموعة بمعامل مُختزَل إلى إشارته.

تحسينات فرميونية وأصيلة للعتاد

تُرجع الدالة generate_preset_jw_pass_manager() كائن MultiStagePassManager يأخذ FermionicCircuit وينتج دارة نهائية محسّنة يمكننا ترجمتها للعمل على عتادنا. نستبدل مرحلة التحسين الافتراضية بـFermionicPassManager يحتوي مسار QDriftTrotterization الخاص بنا:

• يستخدم مسار QDriftTrotterization حساب الأوزان وأخذ العينات داخليًا لتوليد الدارات التي سنستخدمها لأخذ العينات

• مسار RelabelModes هو مسار تحسين آخر يمكن استخدامه لتبديل الأنماط الفرميونية لتحسين الاتصال بين الكيوبتات وتقليل عمق البوابات؛ اقرأ المزيد في مرجع الواجهة البرمجية

المراحل المتبقية من MultiStagePassManager تعمل تلقائيًا وتتولى التخطيط الكامل من الفرميون إلى الكيوبت:

• F2QLayout: يطبق مدير المسارات المُعد مسبقًا مسار TrivialF2QLayout، الذي يخطط بشكل مباشر $n$ بتًا فرميونيًا إلى $n$ كيوبت.

• F2QSynth: مسار تحويل (transpilation) لتخطيط تعليمات الدارة القائمة على الفرميونات إلى تعليمات قائمة على الكيوبتات.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

الآن بعد أن انتهينا من التحسينات على المستوى الفرميوني، يمكننا تحويل (transpile) الدارات للتنفيذ على المحاكي.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

الخطوة 3: التنفيذ باستخدام بدائيات Qiskit

الآن بعد أن أصبحت لدينا دارتنا، يمكننا تشغيلها باستخدام بدائيات Qiskit على AerSimulator. سنجمع كل العدّات من الدارات المختلفة. نحوّلها إلى متجهات منطقية قبل المعالجة النهائية باستخدام SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

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

استخدام السلاسل الثنائية لأجل SQD

يمكننا الآن تشغيل مخطط القطر على السلاسل الثنائية المختارة لإيجاد القيمة الذاتية الأدنى التي ستقابل طاقة الحالة الأساسية للجزيء. ننشئ دالة استدعاء، ونعلن الإشغالات الابتدائية، ونضبط المعاملات قبل تشغيل مخطط القطر أخيرًا. تُستخدم دالة الاستدعاء لطباعة التكرار الحالي وتقدير القيمة الذاتية الحالية عند كل تكرار.

أخيرًا، للحصول على تقدير الحالة الأساسية، نضيف nuclear_repulsion_energy إلى الطاقة الناتجة.

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

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

مثال على العتاد

يستخدم هذا المثال 20 كيوبت (10 مدارات فضائية). هذا الاختيار وسيلة عملية لدرس يُراد له أن يعمل بسرعة، وليس سقفًا صارمًا للطريقة.

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

عمليًا إذن، يعتمد حجم النظام القابل للتحقيق على بُعد الفضاء الفرعي الذي تحتاجه للحصول على الدقة المرغوبة وعلى الذاكرة والأنوية المتوفرة لحلّال القيم الذاتية. عادةً ما تتطلب الفضاءات المدارية الأكبر فضاءً فرعيًا أكبر للوصول إلى الدقة الكيميائية، وهذا ما يحفز في النهاية استخدام موارد موزعة — انظر qiskit-addon-sqd-hpc لتوسيع هذه الخطوة. بدلًا من افتراض حد أقصى ثابت، فإن النهج العملي هو مراقبة بُعد الفضاء الفرعي المُبلَّغ عنه وتقارب الطاقة عبر التكرارات وزيادة حجم الفضاء الفرعي حتى تتوقف الطاقة عن التحسن أو تستنفد الذاكرة المتوفرة.

ملاحظة: بسبب خطأ أخذ العينات الناتج عن الضجيج في العتاد، سيكون الفضاء الفرعي المُنشأ للقطر في تشغيل العتاد أكبر مما نحصل عليه عند استخدام المحاكي. وبينما يزيد ذلك من بُعد الفضاء الفرعي الذي نريد قطره، لا يزال سير العمل يعطينا إجابة دقيقة بفضل متانة SQD تجاه الضجيج.

تقليم السلاسل الزائفة

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

بما أن النيتروجين لا يمكن أن يحتوي إلا على سبعة إلكترونات $\alpha$ وسبعة إلكترونات $\beta$، فإن أي سلاسل ثنائية تحتوي على أكثر أو أقل من سبعة آحاد في النصف الأول والنصف الثاني من المخرجات يمكن التخلص منها. نعرّف دالة تتحقق من صلاحية السلاسل الثنائية، وإذا لم تكن صالحة، تتخلص منها. بمجرد تصفية السلاسل الثنائية الزائفة، تُرسل البقية إلى مخطط القطر. استخدم العلامة PRUNE أدناه للتبديل بين السلوكين.

ضع في اعتبارك أن التقليم ليس سوى أحد عدة خيارات تشكّل الفضاء الفرعي النهائي، إلى جانب عدد الدارات، ومجموعة أزمنة التطور، وتصفية الحدود القطرية. لا تكون مقارنة تشغيل مُقلَّم بآخر غير مُقلَّم مفيدة إلا إذا بقي كل شيء آخر ثابتًا؛ يناقش تبويب C++ هذا بتفصيل أكبر، إذ إنه ينتقي لاحقًا بدلًا من الاسترجاع ويختلف أيضًا في تلك المعاملات الأخرى.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT_SQDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

C++

احصل على شيفرة C++

توجد شيفرة C++ المصدرية ومشروع CMake وتعليمات البناء في مجلد docs/tutorials/assets/sqdrift/cpp في مستودع توثيق Qiskit. ويوجد مدخل FCIDump في المجلد المجاور fcidump_files. انسخ المستودع (clone)، أو نزّل المجلدين معًا، لبناء هذا الدرس التعليمي وتشغيله بنفسك.

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

مخرجات التعلم (C++)

• تعلّم كيفية إنشاء دارات بعمق أصغر مقارنة بـ Trotterization

• اطّلع على سير عمل متكامل من البداية إلى النهاية لتقدير الحالة الأساسية باستخدام qDRIFT وSQD

• تعلّم كيفية استخدام qiskit-fermions جنبًا إلى جنب مع إضافات Qiskit الأخرى لتنفيذ سير عمل كهذا

المتطلبات الأساسية (C++)

• اقرأ نظرة عامة على التقطير الكمومي المعتمد على العينات (SQD)

• اقرأ درس التقطير الكمومي بطريقة Krylov المعتمدة على العينات (SKQD)

الخلفية (C++)

SqDRIFT هي نسخة متغيرة من SKQD تستبدل الحاجة إلى اختيار أنساتز (ansatz) تُؤخذ منه عينات من سلاسل البتات بمجموعة من دارات التطور الزمني المبنية مباشرة من الهاميلتوني المستهدف. يتحقق ذلك بأخذ عينات فرعية من مؤثرات تطور زمني أصغر من الهاميلتوني استنادًا إلى معاملاته، وهو ما يُعرف بطريقة qDRIFT للتروتّرة (Trotterization).

يستخدم كلا التنفيذين Qiskit Fermions لتمثيل الهاملتوني الجزيئي وبناء الدارات لخوارزمية qDRIFT. يستخدم تنفيذ Python مسارات المحوّل (transpiler) الفرميونية، بينما يأخذ تنفيذ C++ عينات المجموعات صراحةً ويربطها بمؤثرات الكيوبت.

ليكن الهاميلتوني بالصيغة التالية:

$$H = \sum_{i=1}^{N} c_i h_i$$

حيث، دون فقدان للعمومية، نشترط $c_i > 0$ وأن تكون القيمة الذاتية الكبرى لـ $h_i$ مساوية، بالقيمة المطلقة، لـ $1$. يُستوعَب أي معامل إشاري أو معقد داخل $h_i$، لذا فإن المعاملات $c_i$ هي أوزان موجبة تمامًا بينما تحمل $h_i$ اتجاه كل حد. هنا $N$ هو عدد الحدود (أو، بعد التجميع، عدد المجموعات) في الهاميلتوني؛ وهو خاصية للهاميلتوني ويختلف عن عدد المؤثرات المأخوذة بالعينة في دارة واحدة، والمكتوب $n$ أدناه.

ثم تُحقق خوارزمية qDRIFT، للوقت المستهدف $t$، مؤثرًا ما $V_k$، حيث يتراوح $k$ من $1 \cdots K$ ويدل على دارة SqDRIFT رقم $k$، معرّفًا كالتالي:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

هنا $n$ هو عدد المؤثرات المأخوذة بالعينة لكل دارة و$K$ هو عدد الدارات في المجموعة. يمتد الجداء عبر السحوبات الـ $n$، وليس عبر جميع حدود الهاميلتوني الـ $N$، ونظرًا لأن الحدود تُسحب مع الإرجاع، يمكن أن يظهر نفس $h_i$ أكثر من مرة في $V_k$ واحد.

الكمية:

$$\lambda = \sum_{i=1}^{N} c_i$$

هي معيار $L_1$ للمعاملات، لذا تتطور كل خطوة من الخطوات الـ $n$ لنفس المدة $\lambda t / n$ بغض النظر عن الحد الذي تم سحبه. تجانس زاوية الخطوة هو السمة المميزة لـ qDRIFT: يؤثر المعامل في النتيجة من خلال عدد مرات سحب حده، لا من خلال مقدار دوران ذلك الحد. تُؤخذ الفهارس بالعينة من التوزيع:

$$P[k_i] = \frac{c_i}{\lambda}$$

لذا فإن المتسلسلة $(k_1, \ldots, k_n)$ هي سلسلة عشوائية من فهارس الحدود مسحوبة من هذا التوزيع. وبما أن $c_i$ موجبة ومجموعها يساوي $\lambda$، فإن هذا توزيع احتمالي مُطبَّع، ويقرِّب توقع القناة الناتجة عبر السحوبات العشوائية التطور تحت $H$، بخطأ يتناقص مع تزايد $n$. لاحظ أن خطأ التقريب يعتمد على $\lambda$ لا على عدد الحدود $N$.

(تكتب ورقة SqDRIFT عدد الحدود بصيغة $\mathcal{N}$ وطول المتسلسلة بصيغة $N$؛ نستخدم هنا $N$ و$n$ للإبقاء على تمييز واضح بينهما.)

يوضح هذا الدرس التعليمي كيفية توليد مجموعة من هذه الدارات العشوائية. وبعد إنشاء هذه الدارات، على غرار كيفية إنشاء فضاء فرعي من نوع Krylov لمؤثرات مختلفة، نأخذ عينات من سلاسل البتات من عدة مؤثرات كهذه بمعاملات زمنية مختلفة. يضمن هذا تداخلاً أعلى بين متجهات الحالة الأساسية وسلاسل البتات المأخوذة بالعينة.

المتطلبات (C++)

يعمل هذا التنفيذ على عتاد IBM Quantum® ويستخدم ملفًا تنفيذيًا خارجيًا لـ SBD في القطرنة الكلاسيكية. ولا يتضمن مثال محاكي Python.

اتبع تعليمات بناء C++ لتثبيت المترجم وCMake وRust وBoost واعتماديات Qiskit. ثبّت Rust باستخدام rustup: تثبّت اعتماديات Qiskit إصدار Rust الخاص بها، وينزّله rustup تلقائيًا، لكن Rust المثبّت من Homebrew أو من مدير حزم Linux قد يكون قديمًا جدًا. وتتطلب خطوة SBD الأخيرة أيضًا MPI وOpenMP وBLAS/LAPACK. نظام Windows غير مدعوم حاليًا في هذا التنفيذ.

استخدم شيفرة C++ المصدرية الكاملة ومشروع CMake من المستودع. احتفظ بمجلد fcidump_files المجاور عند تنزيل المثال. من نسخة مسحوبة (checkout) من مستودع التوثيق، ابنِ باستخدام الأوامر التالية. يقوم البناء الأول بتنزيل كل الاعتماديات وترجمتها ويستغرق عدة دقائق.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

يكتب البناء الملف التنفيذي SqDRIFT في docs/tutorials/assets/sqdrift/cpp (جذر المشروع)، وليس في build/. شغّل ./SqDRIFT من جذر المشروع حتى يتمكن من العثور على ../fcidump_files/N2_sto_3g. وتُشغَّل كل الأوامر الأخرى في هذا التبويب أيضًا من جذر المشروع ما لم يُذكر خلاف ذلك. إذا فشل ./SqDRIFT بخطأ مكتبة مشتركة مفقودة، فاتبع تعليمات مسار المكتبات الخاصة بمنصتك في دليل البناء؛ وعلى macOS، يعمل الملف التنفيذي عادةً بدونها.

الإعداد (C++)

المقتطفات أدناه مأخوذة من برنامج SqDRIFT.cpp الكامل. ابنِ ذلك الملف وشغّله كبرنامج واحد.

قبل تشغيل SqDRIFT، تأكد من أن رمز IBM Quantum API الخاص بك مهيّأ لـ IBM Cloud Platform.

إذا كان لديك Python مع تثبيت qiskit-ibm-runtime، فأسهل طريقة لإنشاء ملف إعداد بالتنسيق الصحيح هي:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

يكتب هذا الملف ~/.qiskit/qiskit-ibm.json (Linux/macOS) أو %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) بالشكل الصحيح. وإذا فضّلت إنشاء الملف يدويًا:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

احصل على مفتاح API الخاص بك من IBM Quantum Platform.

ملاحظة: يتسبب ملف إعداد مفقود أو مشوّه في انهيار Rust من نوع unwrap() غير واضح عند البدء بدلًا من خطأ صريح. تحقق جيدًا من اسم المفتاح (default-ibm-quantum-platform) وقيمة القناة (ibm_quantum_platform) إذا واجهت مشكلة.

التنفيذ (C++)

الخطوة 1: تحويل المدخلات الكلاسيكية إلى مسألة كمومية (C++)

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

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

يحمّل هاملتوني جزيء $N_2$ من ملف FCIDump. يساوي عدد الكيوبتات ضعف عدد المدارات (مغزل علوي ومغزل سفلي).

إنشاء الهاملتوني وترتيبه الطبيعي

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

يحوّل بيانات FCIDump إلى مؤثر فرميوني ويطبّق الترتيب الطبيعي لتبسيط بنية المؤثر.

تجميع الحدود بحسب البنية الإلكترونية

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

يجمّع حدود الهاملتوني المترابطة بحسب البنية الإلكترونية قبل بناء دارات تطورها.

حساب أوزان العينات، ثم التطبيع والربط بمؤثرات الكيوبت

الترتيب مهم هنا. يعيد التطبيع تحجيم كل معامل إلى مقدار الوحدة في مكانه، لذا يجب حساب أوزان العينات أولًا. وإذا جرى التطبيع أولًا، فإن كل $|c_j|$ يُقرأ على أنه 1.0 تمامًا، ويؤول وزن كل مجموعة إلى مجرد عدد حدودها — وعندها يتحدد توزيع العينات بحجم المجموعة بدلًا من المعاملات الفيزيائية.

وزن كل مجموعة هو متوسط القيمة المطلقة لمعاملات حدودها. يستخدم تمرير QDriftTrotterization في تبويب Python الوزن نفسه للهاملتوني المجمّع. أما استخدام المجموع بدلًا من ذلك فسيجعل $\lambda$ يتناسب مع حجم المجموعة، مما يشوّه بدوره زمن التطور $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

الوزن $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ هو متوسط القيمة المطلقة للمعامل على حدود المجموعة $G_i$، مقروءًا من المعاملات الأصلية. ثم يعيد التطبيع تحجيم كل حد فرميوني إلى مقدار الوحدة قبل تحويل Jordan-Wigner، محتفظًا فقط بطوره/إشارته في المؤثر المتطور. وهذا يطابق تطور الحدود المطبَّعة في qDRIFT: تحدد مقادير المعاملات احتمالات أخذ العينات عبر $\lambda$، بينما يستخدم تطور الدارة الحد المطبَّع حتى لا تُحسب المعاملات الكبيرة مرتين.

أخذ عينات من مجموعات المؤثرات

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

يسحب مجموعات المؤثرات لـ 100 دارة عشوائية، تحتوي كل منها على 10 مجموعات مأخوذة من التوزيع الموزون. تحوّل الخطوة التالية هذه المجموعات إلى مجموعة دارات qDRIFT.

الخطوة 2: بناء دارات qDRIFT (C++)

في نسخة C++، تبني هذه الخطوة الدارات من مجموعات المؤثرات المأخوذة. ويجري التحويل (transpilation) للـ Backend المستهدف في الخطوة 3، مباشرةً قبل إرسال كل دارة.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

يطبع البرنامج الكامل بعد ذلك إحصاءات الدارات ومؤشرات المؤثرات للدارة 0، وقد حُذفت هنا.

لكل مؤثر مأخوذ، ينشئ دارة تطور Suzuki-Trotter بزمن $t_k = \frac{\lambda \tau}{n}$، حيث $n$ هو عدد المؤثرات المأخوذة لكل دارة و$\tau$ هو الخطوة الزمنية. ولأن كل مجموعة مأخوذة جرى تطبيعها حدًا حدًا قبل الربط، فإن هذا التطور يطبّق فقط طور/إشارة كل حد أثناء الدوران؛ وتسهم المقادير عبر توزيع أخذ العينات فقط. تُحضَّر حالة Hartree-Fock الابتدائية بتطبيق بوابات X على الكيوبتات 0-6 و10-16 (سبعة إلكترونات في كل قطاع لف). وتُلحَق تعليمات كل دارة تطور يدويًا بالدارة الرئيسية باستخدام واجهة C++ البرمجية.

الخطوة 3: التنفيذ باستخدام بدائيات Qiskit (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

يتصل بـ IBM Quantum، ويختار الـ Backend بالاسم، ويحوّل (transpile) كل دارة لذلك الـ Backend، ويرسل مهمة واحدة لكل دارة: 100 مهمة Sampler منفصلة، كل منها بـ 100 قياس (shots) على سجل من 20 كيوبت. ثم يستعلم عن كل مهمة حتى تنتهي ويجمع السلاسل الثنائية من result[0]. تستهلك كل مهمة وقت QPU وتنتظر في الطابور بشكل منفصل، لذا توقّع أن تستغرق هذه الخطوة وقتًا أطول بكثير من مهمة واحدة مجمّعة.

النتائج على العتاد المقتبسة في هذه الصفحة مأخوذة من تشغيل على ibm_pittsburgh (100 دارة × 100 قياس)، مع تغيير backend_name من القيمة الافتراضية ibm_fez. وجرت المعالجة اللاحقة للعينات المحفوظة بترتيب بتات المحدِّدات المصحَّح الموضح في الخطوة 4. تتفاوت ضوضاء العتاد بين التشغيلات، لذا ستختلف أعدادك وطاقتك النهائية.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

تغيير الـ Backend: غيّر backend_name في SqDRIFT.cpp لاستهداف جهاز مختلف؛ وتسرد service.backends() ما يمكن لحسابك الوصول إليه.

الخطوة 4: المعالجة اللاحقة وإرجاع النتيجة بالصيغة الكلاسيكية المطلوبة (C++)

الانتقاء اللاحق للسلاسل الثنائية وكتابة سلاسل CI

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

ينتقي لاحقًا السلاسل الثنائية ذات وزن Hamming الصحيح (n_alpha إلكترونات مغزل علوي، n_beta إلكترونات مغزل سفلي، مشتقة من حقل NELEC في FCIDump)، ويحوّلها إلى سلاسل تفاعل التكوين (CI) باستخدام تنظير المغزل (spin symmetrization)، ويكتبها في alphadets_from_sqd.txt لاستخدامها في قطرنة الأساس المختار اللاحقة.

تنظيف الموارد

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

يحرر بشكل سليم كل الذاكرة والموارد المخصصة. يُنظَّف std::vector<QuantumCircuit> في C++ تلقائيًا عبر RAII. ويحرر يدويًا مؤثرات الكيوبت والمؤثرات الفرميونية وبيانات FCIDump.

بعد أن أنشأنا الأساس الذي سنسقط عليه الهاملتوني، يمكننا المتابعة إلى عملية القطرنة للحصول على تقدير الحالة الأساسية. شغّل الأوامر التالية من جذر المشروع، docs/tutorials/assets/sqdrift/cpp.

بناء الملف التنفيذي diag الخاص بـ SBD

يجلب CMake مكتبة SBD كمصدر فقط إلى deps/sbd/ (ويتخطى البناء الرئيسي ترجمتها عمدًا). يستخدم تطبيق SBD ملف Makefile عاديًا، وتتطلب SBD نفسها MPI وOpenMP وBLAS/LAPACK — لذا فإن make العادي لا يكفي حتى تُثبَّت هذه ويتطابق ملف Configuration الخاص بالتطبيق مع سلسلة أدواتك.

يقرأ Makefile رايات الترجمة والربط من ملف Configuration في المجلد نفسه. النسخة المرفقة في المصدر الأصلي مضبوطة لتخطيط macOS/Homebrew محدد ولن تُبنى كما هي على جهاز عادي — وأبرز ذلك أن -I/opt/homebrew/opt/llvm/include لا يحتوي على omp.h (فهو موجود في مجلد موارد clang خاص بالإصدار)، كما أن clang++ من Apple يرفض -fopenmp تمامًا:

clang++: error: unsupported option '-fopenmp'

الوصفة التالية جرى التحقق من عملها على macOS (Apple Silicon، Darwin 25.6، Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — مُختبَر

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

ثم استبدل deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration بما يلي:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

يطبع البناء الناجح سطري الترجمة/الربط ويترك ملف diag التنفيذي بجانب Makefile. وتحذير ld: warning: ignoring duplicate libraries: '-lomp' الواحد غير ضار.

Intel macOS: بادئة Homebrew هي /usr/local بدلًا من /opt/homebrew، لذا استبدلها في المسارين السابقين.

Linux (Ubuntu/Debian) — غير مُختبَر

تغطي المتطلبات المسبقة التالية MPI وOpenMP وBLAS/LAPACK. يدعم GCC الخيار -fopenmp أصلًا، لذا لا حاجة لتجاوز المترجم. لم يُتحقق من هذا على Linux:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

شغّل SqDRIFT لتوليد سلاسل CI

إذا كنت قد شغّلت ./SqDRIFT بالفعل عبر الخطوات 1–4، فإن alphadets_from_sqd.txt موجود بالفعل في جذر المشروع ويمكنك تخطي هذه الخطوة. وإلا، فشغّله من جذر المشروع:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

تشغيل القطرنة

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

المعاملات:

• --fcidump: مسار ملف FCIDump الذي يحتوي على الهاملتوني
• --adetfile: مسار ملف محدِّدات ألفا (سلاسل CI من SqDRIFT)
• --method 0: طريقة قطرنة Davidson
• --iteration 100: الحد الأقصى لعدد تكرارات Davidson
• --block 10: حجم الكتلة لخوارزمية Davidson
• --tolerance 1e-8: سماحية التقارب للطاقة

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

السطر Energy = -107.5108722590187 هو تقدير الحالة الأساسية على الفضاء الفرعي المُعايَن.

كيف يختلف هذا عن مثال Python

يختلف المثالان في عدة جوانب في آن واحد، وليس فقط في طريقة معالجة السلاسل الثنائية ذات التماثل الخاطئ، لذا لا يُتوقع أن ينتجا الطاقات نفسها:

	تبويب C++	تبويب Python
السلاسل الثنائية ذات التماثل الخاطئ	تُستبعد بواسطة postselect_bitstrings بحسب وزن Hamming	تُصلَح بالاسترجاع التكراري للتكوينات
بنية الحلقة	تمريرة واحدة: أخذ عينات ← انتقاء لاحق ← قطرنة مرة واحدة	حلقة خارجية: قطرنة، قراءة إشغالات المدارات، استرجاع التكوينات، إعادة القطرنة
عدد الدارات	100 دارة (num_circuits = 100)	400 دارة (200 لكل زمن تطور)
أزمنة التطور	$\tau$ واحد (time_step = 1)	زمنان، times = [1.0, 10.0]
الحدود القطرية	تبقى في المؤثر المأخوذ	تُزال بواسطة filter_diagonal_terms، ما يغيّر توزيع أخذ العينات والدارات
ترتيب الحدود	مجمّعة فقط	يُطبَّق canonical_order قبل التجميع لضمان إمكانية إعادة الإنتاج
حلّال القطرنة	SBD (diag، ملف تنفيذي خارجي MPI/OpenMP)	حلّال qiskit-addon-sqd داخل العملية
Backend	محدد بالاسم (backend_name = "ibm_fez"؛ واستُخدم ibm_pittsburgh في التشغيل المعروض في هذه الصفحة)	Backend تشغيلي من نوع least_busy به عدد كافٍ من الكيوبتات
التحويل (Transpilation)	كل دارة على حدة، مستوى التحسين 2	كل الدارات معًا، مستوى التحسين 3
إرسال المهام	100 مهمة، دارة واحدة لكل منها	مهمة Sampler واحدة تحتوي على كل الدارات الـ 400
تشغيل المحاكي	لا يوجد؛ عتاد فقط	مثال محاكي خالٍ من الضوضاء قبل التشغيل على العتاد

يستخدم كلا الإصدارين الهاملتوني نفسه، ويجمّعان حدوده بحسب البنية الإلكترونية، ويزنان كل مجموعة بمتوسط القيمة المطلقة لمعاملات حدودها، ويبدآن من حالة Hartree-Fock، ويأخذان عينات من 10 مجموعات و100 قياس لكل دارة.

يؤثر كل اختلاف على الفضاء الفرعي المُعايَن، وبالتالي على الطاقة، بشكل مستقل. فالمجموعة الأكبر من الدارات وزمن التطور الثاني الأطول كلاهما يوسّعان التكوينات التي يستكشفها مثال Python؛ ويمنح أخذ العينات عند أكثر من زمن واحد SqDRIFT تغطيته الشبيهة بـ Krylov للفضاء منخفض الطاقة. أما تصفية الحدود القطرية فتغيّر المؤثرات التي يمكن لأداة أخذ العينات سحبها وأوزانها، لذا يأخذ التشغيلان عينات من توزيعين مختلفين.

الانتقاء اللاحق مرشّح صارم، ويُسقَط القياس ذو عدد الجسيمات الخاطئ تمامًا. في التشغيل أعلاه أدى ذلك إلى التخلص من معظم البيانات، إذ نجا 1137 قياسًا من أصل 10000 من فحص وزن Hamming (7,7) (نحو 11%)، لتنهار في 49 محدِّد $\alpha$ متمايزًا فقط. أما استرجاع التكوينات فيستخدم متوسط إشغالات المدارات من قطرنة سابقة لقلب البتات وإصلاح هذه القياسات، معيدًا تدوير ما كان الانتقاء اللاحق سيتخلص منه. وتجعله هذه التغذية الراجعة تكراريًا بطبيعته: إذ تتطلب recover_configurations في qiskit-addon-sqd-hpc وسيطة avg_occupancies التي لا توجد إلا بعد تشغيل قطرنة.

ينتج تشغيل C++ هذا فضاءً فرعيًا من 49 محدِّد $\alpha$ وطاقة قدرها -107.5109 Ha. وفضاء محدِّدات $\alpha$ الكامل لسبعة إلكترونات في 10 مدارات هو فقط $\binom{10}{7} = 120$ سلسلة، لذا فإن تغذية الـ 120 كلها إلى الملف التنفيذي diag نفسه تعطي طاقة FCI الدقيقة -107.6482 Ha. ولذلك يقع التشغيل على بُعد نحو 0.1373 Ha فوق القيمة الدقيقة بينما يمتد على 49 من محدِّدات $\alpha$ الـ 120. وبما أن SQD تغايري في الفضاء الفرعي المختار، فإن الفضاء الفرعي الأصغر لا يمكنه إلا رفع الطاقة، لذا فإن الفضاء الفرعي غير المكتمل هو التفسير الأكثر مباشرة للفجوة.

لا يمكن إجراء فحص FCI الدقيق هذا إلا لأن المثال صغير؛ وهو ليس جزءًا من سير عمل SqDRIFT، وفي أحجام الأنظمة التي تحتاج فعلًا إلى SQD لا يوجد مرجع كهذا. ويفيد هنا فقط في قياس مقدار ما استرجعه هذا التشغيل من الفضاء.

المقارنة مع مثال Python غير مضبوطة، لذا لا ينبغي عزو أي فرق إلى غياب استرجاع التكوينات وحده؛ لاحظ أن عدد الدارات وزمن التطور الواحد وتصفية الحدود القطرية وضوضاء العتاد كلها تختلف في آن واحد وكلها تشكّل الفضاء الفرعي.

لسد الفجوة، إليك بعض الخطوات التالية الممكنة بعد إكمال هذا الدرس التعليمي: ارفع shots أو num_circuits أو ops_per_circuit؛ وخذ عينات عند أكثر من time_step واحد؛ وأضف حلقة استرجاع فوق Qiskit::addon::sqd::recover_configurations.

الخطوات التالية

توصيات

• القطرنة الكمومية Krylov المعتمدة على العينات لنموذج شبكي فرميوني — درس تعليمي ذو صلة يستخدم دارات التطور الزمني بدلًا من ansatz تغايري.

• القطرنة الكمومية المعتمدة على العينات لهاملتوني كيميائي — درس تعليمي حول كيفية بناء دارة LUCJ (local unitary cluster Jastrow) لمحاكاة الكيمياء الكمومية.

• ورقة SqDRIFT — المرجع العلمي الذي يستند إليه هذا الدرس التعليمي. (لاحظ أن بعض التحسينات المناقشة في هذه الورقة لا تزال قيد العمل حاليًا، وأن هذا الدرس التعليمي قابل للتغيير مستقبلًا بناءً على تطور المكتبات المستخدمة.)