Алгоритм SqDRIFT для оцінки основного стану

Python

Завантаж ноутбук Python

Завантаж ноутбук Python, щоб самостійно виконати цей туторіал. Ноутбук завантажує вхідні дані FCIDump, якщо їх ще немає.

Оцінка часу: 180 секунд на процесорі Heron r3 (ПРИМІТКА: це лише оцінка. Твій час виконання може відрізнятися.)

Результати навчання

• Дізнайся, як створювати схеми з меншою глибиною порівняно з тротеризацією

• Пройди наскрізний робочий процес для оцінки основного стану, використовуючи qDRIFT та SQD

• Дізнайся, як використовувати qiskit-fermions разом з іншими аддонами Qiskit для реалізації такого робочого процесу

Передумови

• Прочитай огляд Семпл-орієнтованої квантової діагоналізації (SQD)

• Прочитай урок Семпл-орієнтованої крилівської квантової діагоналізації (SKQD)

Довідкова інформація

SqDRIFT — це варіант SKQD, який усуває потребу вибирати анзац, із якого семплюються бітові рядки, замінюючи його ансамблем схем часової еволюції, побудованих безпосередньо з цільового гамільтоніана. Це досягається підвибіркою менших операторів часової еволюції з гамільтоніана на основі його коефіцієнтів, що відомо як метод тротеризації qDRIFT.

Обидві реалізації використовують Qiskit Fermions, щоб представити молекулярний гамільтоніан і побудувати схеми для алгоритму qDRIFT. Реалізація на Python використовує ферміонні проходи транспілятора, тоді як реалізація на 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$ і позначає $k$-ту схему 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

Налаштування

# 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 передбачає члени в нормальному впорядкуванні.

Фільтрація діагональних термів

Ми видаляємо діагональні члени з гамільтоніана, що використовується для генерування схем, щоб $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.

Підготовка анзацу

Ми готуємо стан Хартрі-Фока за допомогою класу 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$$

Оскільки доданки було згруповано на кроці 1, кожен $h_i$ тут є цілою групою: $c_i$ — це середній модуль коефіцієнтів доданків групи $i$, а кожен доданок групи еволюціонує з коефіцієнтом, зведеним до його знака.

Ферміонні та апаратно-специфічні оптимізації

Функція generate_preset_jw_pass_manager() повертає MultiStagePassManager, який приймає FermionicCircuit і генерує оптимізовану кінцеву схему, яку ми можемо транспілювати для виконання на нашому обладнанні. Ми замінюємо його стандартний етап оптимізації на FermionicPassManager, що містить наш пас QDriftTrotterization:

• Пас QDriftTrotterization внутрішньо використовує обчислення ваги та вибірку для генерування схем, які ми використаємо для вибірки

• Пас RelabelModes — це ще один пас оптимізації, який можна використати для перестановки ферміонних мод, щоб оптимізувати зв'язність між кубітами та зменшити глибину гейтів; докладніше читай у довіднику API

Решта етапів 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))

400

Тепер, коли ми завершили оптимізації на ферміонному рівні, ми можемо транспілювати схеми для виконання на симуляторі.

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. Клонуй репозиторій або завантаж обидва каталоги, щоб самостійно зібрати й запустити цей туторіал.

Оцінка часу: 180 секунд на процесорі Heron r3 (ПРИМІТКА: це лише оцінка. Твій час виконання може відрізнятися.)

Результати навчання (C++)

• Дізнайся, як створювати схеми з меншою глибиною порівняно з тротеризацією

• Пройди наскрізний робочий процес для оцінки основного стану, використовуючи qDRIFT та SQD

• Дізнайся, як використовувати qiskit-fermions разом з іншими аддонами Qiskit для реалізації такого робочого процесу

Передумови (C++)

• Прочитай огляд Семпл-орієнтованої квантової діагоналізації (SQD)

• Прочитай урок Семпл-орієнтованої крилівської квантової діагоналізації (SKQD)

Довідкова інформація (C++)

SqDRIFT — це варіант SKQD, який усуває потребу вибирати анзац, із якого семплюються бітові рядки, замінюючи його ансамблем схем часової еволюції, побудованих безпосередньо з цільового гамільтоніана. Це досягається підвибіркою менших операторів часової еволюції з гамільтоніана на основі його коефіцієнтів, що відомо як метод тротеризації qDRIFT.

Обидві реалізації використовують Qiskit Fermions, щоб представити молекулярний гамільтоніан і побудувати схеми для алгоритму qDRIFT. Реалізація на Python використовує ферміонні проходи транспілятора, тоді як реалізація на 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$ і позначає $k$-ту схему 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$, щоб чітко розрізнити ці два поняття.)

Цей туторіал показує, як згенерувати ансамбль таких рандомізованих схем. Після того як ми створили ці схеми, подібно до того, як ми створюємо крилівський підпростір для різних операторів, ми семплюємо бітові рядки з кількох таких операторів із різними параметрами часу. Це забезпечує більше перекриття між векторами основного стану та семпльованими бітовими рядками.

Вимоги (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 під час завантаження прикладу. З клону репозиторію документації збери проєкт такими командами. Перше збирання завантажує й компілює всі залежності та триває кілька хвилин.

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, переконайся, що твій токен API IBM Quantum налаштовано для 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$, зчитаним з оригінальних коефіцієнтів. Потім нормалізація перемасштабовує кожен ферміонний доданок до одиничного модуля перед перетворенням Джордана-Вігнера, зберігаючи в еволюціонованому операторі лише його фазу/знак. Це відповідає еволюції нормалізованих доданків у 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++ цей крок будує схеми з відібраних наборів операторів. Транспіляція для цільового бекенда відбувається на кроці 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, які тут опущено.

Для кожного відібраного оператора створює схему еволюції Сузукі-Троттера з часом $t_k = \frac{\lambda \tau}{n}$, де $n$ — кількість операторів, відібраних на схему, а $\tau$ — крок часу. Оскільки кожну відібрану групу було нормалізовано по доданках перед відображенням, ця еволюція застосовує під час обертання лише фазу/знак кожного доданка; величини вносять внесок лише через розподіл вибірки. Початковий стан Гартрі-Фока готується застосуванням вентилів X до кубітів 0-6 і 10-16 (сім електронів у кожному спіновому секторі). Інструкції кожної схеми еволюції вручну додаються до головної схеми за допомогою API 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, вибирає бекенд за назвою, транспілює кожну схему для цього бекенда та надсилає одне завдання на схему: 100 окремих завдань Sampler по 100 шотів кожне на 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_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!

Виконує постселекцію бітових рядків із правильною вагою Геммінга (n_alpha електронів зі спіном угору, n_beta зі спіном униз, виведених із поля NELEC у FCIDump), перетворює їх на рядки конфігураційної взаємодії (CI) за допомогою спінової симетризації та записує їх у 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;

Належним чином звільняє всю виділену пам'ять і ресурси. C++ std::vector<QuantumCircuit> автоматично очищується через RAII. Вручну звільняє кубітні оператори, ферміонні оператори та дані FCIDump.

Тепер, коли ми створили базис, на який проєктуватимемо наш гамільтоніан, можна перейти до діагоналізації, щоб отримати оцінку основного стану. Виконуй наведені нижче команди з кореня проєкту, docs/tutorials/assets/sqdrift/cpp.

Збирання бінарного файлу SBD diag

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: метод діагоналізації Девідсона
• --iteration 100: максимальна кількість ітерацій Девідсона
• --block 10: розмір блоку для алгоритму Девідсона
• --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 за вагою Геммінга	Виправляються ітеративним відновленням конфігурацій
Структура циклу	Один прохід: вибірка → постселекція → одна діагоналізація	Зовнішній цикл: діагоналізація, зчитування орбітальних заселеностей, відновлення конфігурацій, повторна діагоналізація
Кількість схем	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_name = "ibm_fez"; запуск на цій сторінці використовував ibm_pittsburgh)	Робочий бекенд least_busy із достатньою кількістю кубітів
Транспіляція	Кожна схема окремо, рівень оптимізації 2	Усі схеми разом, рівень оптимізації 3
Надсилання завдань	100 завдань, по одній схемі	Одне завдання Sampler з усіма 400 схемами
Запуск на симуляторі	Немає; лише апаратне забезпечення	Приклад на безшумовому симуляторі перед апаратним запуском

Обидві версії використовують той самий гамільтоніан, групують його доданки за електронною структурою, зважують кожну групу за середнім модулем коефіцієнтів її доданків, починають зі стану Гартрі-Фока та відбирають 10 груп і 100 шотів на схему.

Кожна відмінність впливає на відібраний підпростір, а отже й на енергію, незалежно. Більший ансамбль і другий, довший час еволюції обидва розширюють конфігурації, які досліджує приклад на Python; вибірка за більш ніж одного часу дає SqDRIFT його схоже на Крилова покриття низькоенергетичного простору. Фільтрація діагональних доданків змінює, які оператори може добирати семплер, і їхні ваги, тож два запуски відбирають із різних розподілів.

Постселекція — суворий фільтр, і шот із неправильною кількістю частинок відкидається повністю. У наведеному вище запуску це відкинуло більшість даних: 1137 із 10000 шотів пройшли перевірку ваги Геммінга (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 зі 120 $\alpha$-детермінантів. Оскільки SQD є варіаційним у вибраному підпросторі, менший підпростір може лише підвищити енергію, тож неповний підпростір є найпрямішим поясненням цього розриву.

Ця перевірка точним FCI можлива лише тому, що приклад малий; вона не є частиною робочого процесу SqDRIFT, а для розмірів систем, де SQD справді потрібен, такого еталона не існує. Тут вона слугує лише для оцінки того, яку частину простору охопив цей запуск.

Порівняння з прикладом на Python не є контрольованим, тому різницю не слід приписувати лише відсутності відновлення конфігурацій; зауваж, що кількість схем, єдиний час еволюції, фільтрація діагональних доданків і апаратний шум відрізняються одночасно, і всі вони формують підпростір.

Щоб закрити цей розрив, ось кілька можливих наступних кроків, які можна спробувати після завершення цього туторіалу: збільш shots, num_circuits або ops_per_circuit; роби вибірку за більш ніж одного time_step; і додай цикл відновлення поверх Qiskit::addon::sqd::recover_configurations.

Наступні кроки

Рекомендації

• Квантова діагоналізація Крилова на основі вибірки для ферміонної ґраткової моделі — пов'язаний туторіал, що використовує схеми еволюції в часі замість варіаційного ансатцу.

• Квантова діагоналізація на основі вибірки для гамільтоніана хімічної системи — туторіал про те, як побудувати схему локального унітарного кластерного Джастроу (LUCJ) для симуляції квантової хімії.

• Стаття про SqDRIFT — література, на якій ґрунтується цей туторіал. (Зауваж, що деякі оптимізації, обговорені в цій статті, наразі перебувають у розробці, і цей туторіал може змінитися в майбутньому залежно від розвитку використаних бібліотек.)