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

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

Шукаєш версію C++?

Цей туторіал використовує Python. Для реалізації на C++, включно з вихідним кодом та інструкціями зі збірки, дивись туторіал C++ SqDRIFT.

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

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

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

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

Цей туторіал представлений у вигляді блокнота Python для навчальних цілей.

Передумови

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

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

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

SqDRIFT — це варіант SKQD, який усуває потребу вибирати анзац, із якого семплюються бітові рядки, замінюючи його ансамблем схем часової еволюції, побудованих безпосередньо з цільового гамільтоніана. Це досягається підвибіркою менших операторів часової еволюції з гамільтоніана на основі його коефіцієнтів, що відомо як метод тротеризації 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$-ту схему 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 — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Приклад із симулятором

Крок 1: Відображення класичних вхідних даних у квантову задачу

Читання та підготовка FCIDump

Для цього туторіалу ми завантажимо гамільтоніан електронної структури для азоту (N2). Існують і інші способи створення ферміонних операторів. Дивись документацію на qiskit_fermions.operators.library.

Про цей FCIDump. Файл N2_sto_3g описує молекулу азоту ($N_2$) у мінімальному базисі STO-3G при міжатомній відстані 1.09 $\AA$ — експериментальній рівноважній довжині зв'язку. Його заголовок оголошує NORB=10, NELEC=14 та MS2=0: 10 просторових орбіталей (отже, 20 спінових орбіталей і 20 кубітів за перетворенням Джордана-Вігнера), 14 електронів у спіновому синглеті, тобто сім $\alpha$- та сім $\beta$-електронів. Усі орбіталі мають мітку симетрії 1, тобто симетрія точкової групи не використовується. Оскільки це повнопросторовий дамп STO-3G, жодні орбіталі не заморожені, а простір кореляції достатньо малий, щоб точну еталонну енергію FCI можна було обчислити класично для порівняння, як показано в наступній комірці.

Еквівалентний файл можна перегенерувати за допомогою PySCF:

from pyscf import gto, scf, tools

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

Оскільки інтеграли залежать від збіжних орбіталей SCF, перегенерований файл може відрізнятися від наданого фазою чи порядком орбіталей; на повні енергії це не впливає.

Отримання файлу. Знайдіть FCIDump у цьому репозиторії GitHub. Ви можете запустити комірку нижче, щоб завантажити його в місце, яке очікує решта частини навчального посібника.

Спочатку ми використовуємо cisolver, наданий pyscf, щоб отримати еталонну енергію. Це справжня енергія основного стану молекули, з якою ми працюємо. Для цього ми спершу оголосимо norb та nelec — кількість орбіталей та кількість електронів відповідно. Потім ми оголосимо h1e та h2e — одно- та двоелектронні інтеграли відповідно. Усе це пізніше також буде використано для SQD.

import os
from urllib.request import urlopen

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

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Завантаження гамільтоніана

Маючи необхідні дані, ми зчитуємо гамільтоніан з файлу FCI у форматі, сумісному з qiskit-fermions

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

Ферміонні робочі процеси з qiskit-fermions

Спочатку ми відобразимо гамільтоніан у модель ферміонного схемного представлення за допомогою qiskit-fermions, яка надає транспілерні паси та гейти, специфічні для ферміонних схем. Вони будуть використані пізніше, перед традиційними транспілерними пасами Qiskit у цьому робочому процесі.

Групування термів

Щоб забезпечити відтворюваність результатів, ми спочатку використовуємо canonical_order для сортування термів лише на основі їхньої структури. Порядок операторів у списку canon є, таким чином, фіксованим. Це забезпечує відтворюваність створених операторів, оскільки пас QDriftTrotterization, який ми використаємо пізніше, вибирає випадкові індекси для створення операторів qDRIFT.

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

qiskit-fermions надає функцію group_terms_by_electronic_structure, яка виконує це групування за нас.

Зауважте, що group_terms_by_electronic_structure передбачає терми в нормальному впорядкуванні.

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

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

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

Функція 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

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

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

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

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

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

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

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

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

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

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

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

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

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

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

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

Якщо ця робота вас зацікавила, вас можуть зацікавити такі матеріали:

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