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

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

نتائج التعلّم

• تعلّم كيف تنشئ دارات بعمق أصغر مقارنةً بتقريب تروتر

• اتبع سير عمل متكاملًا لتقدير الحالة الأرضية باستخدام qDRIFT وSQD

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

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

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

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

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

الخلفية

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

يستخدم هذا الدرس 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$ ويدل على الدارة $k_{th}$ من SqDRIFT، ويُعرَّف بـ:

$$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$ هنا لإبقاء الاثنين متمايزين بوضوح.)

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

المتطلبات

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

• بيئة 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 — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# 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 = "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 = "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، الذي يوفر تمريرات محوّل (transpiler) وبوابات خاصة بالدارات الفرميونية. ستُستخدم هذه لاحقًا قبل تمريرات المحوّل التقليدية في Qiskit لسير العمل هذا.

تجميع الحدود

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

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

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

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

ترشيح الحدود القطرية

نزيل الحدود القطرية من الهاميلتوني المستخدم لتوليد الدارات، حتى تُنفَق خانات أخذ العيّنات $n$ في qDRIFT على حدود تنقل الإشغال بين الإعدادات. الأفضل ترشيح هذه الحدود من الهاميلتوني عند هذه النقطة، قبل بناء بوابة 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))

5060

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

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

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

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

تحضير الـ Ansatz

نحضّر حالة هارتري-فوك باستخدام الصنف 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 للحدود بشكل عشوائي، باحتمالات تتناسب مع معاملاتها في الهاملتوني. يقوم تمريرة الـ Transpiler الخاصة بـ qDRIFT بذلك نيابةً عنا. يمكننا الآن إنشاء دوائر أقل عمقًا يمكن تنفيذها على العتاد بكفاءة أكبر رغم محدودية اتصال الكيوبتات، حتى عندما يحتوي الهاملتوني على اقترانات بعيدة المدى وحدود من رتبة أعلى من التربيعية. بعد تجميع الحدود، يأخذ عيّنات من المعاملات (operators) بناءً على أوزانها. لكل معامل $h_i$ يُعرَّف الوزن $W_{h_i}$ كما يلي:

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

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

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

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

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

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

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

• F2QLayout: يطبّق مدير التمريرات المسبق الإعداد تمريرة TrivialF2QLayout، التي تعيّن $n$ من البتات الفرميونية إلى $n$ كيوبت بشكل بديهي.

• F2QSynth: تمريرة transpile لتحويل تعليمات الدائرة المبنية على الفرميونات إلى تعليمات مبنية على الكيوبتات.

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: التنفيذ باستخدام primitives الخاصة بـ Qiskit

بعد أن أصبحت لدينا الدوائر، يمكننا تشغيلها باستخدام primitives الخاصة بـ Qiskit على AerSimulator. سنجمع كل العدّات من الدوائر المختلفة، ونحوّلها إلى متجهات منطقية (boolean) قبل المعالجة اللاحقة النهائية باستخدام 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

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

أخيرًا، للحصول على تقدير الحالة الأساسية، نضيف 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$، يمكن تجاهل أي سلاسل بتات تحتوي على أكثر أو أقل من سبعة 1 في النصف الأول والنصف الثاني من المخرجات. نعرّف دالة تتحقق مما إذا كانت سلاسل البتات صالحة، وتتجاهلها إن لم تكن كذلك. بعد تصفية السلاسل الزائفة، تُرسل البقية إلى مخطط القطرنة. استخدم العلَم 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")

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

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

توصيات

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

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