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

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

تبحث عن نسخة ++C؟

يستخدم هذا الدرس التعليمي لغة Python. للحصول على التنفيذ بلغة ++C، بما في ذلك الكود المصدري وتعليمات البناء، راجع درس ++SqDRIFT C التعليمي.

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

تعلّم كيفية إنشاء دارات بعمق أصغر مقارنة بـ Trotterization
اطّلع على سير عمل متكامل من البداية إلى النهاية لتقدير الحالة الأساسية باستخدام qDRIFT وSQD
تعلّم كيفية استخدام qiskit-fermions جنبًا إلى جنب مع إضافات Qiskit الأخرى لتنفيذ سير عمل كهذا

يُقدَّم هذا الدرس التعليمي كدفتر Python لأغراض تعليمية.

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

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

الخلفية

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

يستخدم هذا الدرس التعليمي Qiskit Fermions لإنشاء الدارات الفرميونية الأكثر طبيعية لخوارزمية qDRIFT، يتبعها استخدام تمريرات التخطيط والتوليف الفرميونية قبل إدخال الدارات في خط أنابيب Qiskit التقليدي لتنفيذها على العتاد.

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

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

الإعداد

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

# 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: تحويل المدخلات الكلاسيكية إلى مسألة كمومية

قراءة وتحضير 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 = "assets/sqdrift/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}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/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")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
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 على المؤثر المصفّى بدلًا من الأصلي. هذا مقبول هنا لأن الدارات ليست سوى مرشد ابتكاري لأخذ العينات يُستخدم لاقتراح التهيئات: لا يُفقد أي حد من تقدير الطاقة نفسه، إذ تنطبق التصفية فقط على الهاملتوني المستخدم لبناء الدارات، بينما يستخدم القطر الكلاسيكي لاحقًا الهاملتوني الكامل، بما في ذلك الحدود القطرية. تعتمد دقة 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))

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

عدد الدارات المراد توليدها: 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 لإجراء تحسينات على المستوى الفرميوني، ثم ترجمة دارتنا للنظام الخلفي المختار. بما أن هذه تجربة محاكاة، فسنقوم بذلك أولًا لأجل AerSimulator. حساب الوزن لكل مجموعة

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

W_{h_i} = |c_i| / \lambda

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

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

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

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

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

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

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

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 = "assets/sqdrift/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")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
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

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

توصيات

إذا وجدت هذا العمل مثيرًا للاهتمام، فقد تهتم بالمواد التالية:

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

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

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

الخلفية​

المتطلبات​

الإعداد​

مثال على المحاكي​

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

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

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

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

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

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