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

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

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

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

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

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

Цей туторіал подано як ноутбук Python з навчальною метою.

Передумови

• Прочитай огляд Sample-based quantum diagonalization (SQD)

• Прочитай урок Sample-based Krylov Quantum Diagonalization (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_{th}$ схему SqDRIFT, визначену як:

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

Тут $n$ — кількість вибраних операторів на схему, а $K$ — кількість схем в ансамблі. Добуток іде за $n$ виборами, а не за всіма $N$ доданками гамільтоніана, і оскільки доданки вибираються з поверненням, той самий $h_i$ може з'явитися більше одного разу в одному $V_k$.

Величина:

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

є $L_1$-нормою коефіцієнтів, тож кожен із $n$ кроків еволюціонує протягом однакової тривалості $\lambda t / n$ незалежно від того, який доданок вибрано. Однорідність кута кроку є характерною рисою qDRIFT: коефіцієнт впливає на результат через те, як часто вибирається його доданок, а не через те, наскільки сильно цей доданок повернуто. Індекси вибираються з розподілу:

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

тож послідовність $(k_1, \ldots, k_n)$ є випадковою послідовністю індексів доданків, вибраних із цього розподілу. Оскільки $c_i$ додатні та в сумі дають $\lambda$, це нормований розподіл імовірностей, а математичне сподівання отриманого каналу за випадковими виборами наближає еволюцію під дією $H$ з похибкою, що зменшується зі зростанням $n$. Зауваж, що похибка наближення залежить від $\lambda$, а не від кількості доданків $N$.

(У статті про SqDRIFT кількість доданків позначено $\mathcal{N}$, а довжину послідовності — $N$; ми використовуємо тут $N$ і $n$, щоб чітко розрізняти їх.)

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

Вимоги

Перш ніж починати цей туторіал, переконайся, що встановлено

• Віртуальне середовище Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (зауваж, що назва у множині)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Встановити всі потрібні пакети можна так:

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

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

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

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

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

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

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

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

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

from pyscf import gto, scf, tools

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

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

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

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

import os
from urllib.request import urlopen

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

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

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

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

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

Спочатку ми відобразимо гамільтоніан на ферміонну модель схем за допомогою qiskit-fermions, яка надає проходи транспайлера та вентилі, специфічні для ферміонних схем. Їх буде використано перед традиційними проходами транспайлера 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")

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