SQD для оцінки енергії хімічного гамільтоніана

У цьому уроці ми застосуємо SQD для оцінки енергії основного стану молекули.

Зокрема, ми розглянемо такі теми, використовуючи $4$-кроковий підхід шаблону Qiskit:

1. Крок 1: Відображення задачі на квантові схеми та оператори
   • Налаштування молекулярного гамільтоніана для $N_2$.
   • Пояснення натхненного хімією й апаратно сумісного локального унітарного кластерного анзацу Жастрова (LUCJ) [1]
2. Крок 2: Оптимізація під цільове залізо
   • Оптимізація кількості вентилів і розкладки анзацу для виконання на залізі
3. Крок 3: Виконання на цільовому залізі
   • Запуск оптимізованої схеми на реальному QPU для генерування зразків підпростору.
4. Крок 4: Постобробка результатів
   • Ознайомлення з самоузгодженим циклом відновлення конфігурацій [2]
     • Постобробка повного набору бітрядкових зразків з використанням апріорного знання про кількість частинок і середню орбітальну заповненість, обчислену на останній ітерації.
     • Ймовірнісне створення пакетів підзразків із відновлених бітрядків.
     • Проекція та діагоналізація молекулярного гамільтоніана для кожного зразкового підпростору.
     • Збереження мінімальної енергії основного стану серед усіх пакетів та оновлення середньої орбітальної заповненості.

Протягом уроку ми використовуватимемо кілька програмних пакетів.

• PySCF — для визначення молекули та налаштування гамільтоніана.
• Пакет ffsim — для побудови анзацу LUCJ.
• Qiskit — для транспіляції анзацу під виконання на залізі.
• Qiskit IBM Runtime — для виконання схеми на QPU та збору зразків.
• Qiskit addon SQD — відновлення конфігурацій та оцінка енергії основного стану за допомогою проекції підпростору та матричної діагоналізації.

1. Відображення задачі на квантові схеми та оператори

Молекулярний гамільтоніан

Молекулярний гамільтоніан має загальний вигляд:

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$ — феміонні оператори народження/знищення, пов'язані з $p$-м елементом базисного набору та спіном $\sigma$. $h_{pr}$ і $(pr|qs)$ — одно- та двочастинкові електронні інтеграли. За допомогою pySCF ми визначимо молекулу та обчислимо одно- й двочастинкові інтеграли гамільтоніана для базисного набору 6-31g.

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

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

У цьому уроці ми використаємо перетворення Йордана–Вігнера (JW) для відображення феміонної хвильової функції на кубітову, щоб її можна було підготувати за допомогою квантової схеми. Перетворення JW відображає простір Фока феміонів у M просторових орбіталях на гільбертів простір 2M кубітів, тобто просторова орбіталь розбивається на дві спін-орбіталі: одна пов'язана з електроном зі спіном вгору ($\alpha$), інша — зі спіном вниз ($\beta$). Спін-орбіталь може бути зайнятою або вільною. Зазвичай, коли ми говоримо про кількість орбіталей, мається на увазі кількість просторових орбіталей. Кількість спін-орбіталей буде вдвічі більшою. У квантових схемах кожна спін-орбіталь представлена одним кубітом. Отже, один набір кубітів представлятиме спін-вгору або $\alpha$-орбіталі, а інший — спін-вниз або $\beta$-орбіталі. Наприклад, молекула $N_2$ для базисного набору 6-31g має $16$ просторових орбіталей (тобто $16$ $\alpha$ + $16$ $\beta$ = $32$ спін-орбіталі). Отже, нам знадобиться квантова схема на $32$ кубіти (можливо, з додатковими допоміжними кубітами, як обговорюватиметься пізніше). Кубіти вимірюються у обчислювальному базисі для отримання бітрядків, які представляють електронні конфігурації або детермінанти (Слейтера). Протягом цього уроку ми будемо використовувати терміни «бітрядки», «конфігурації» та «детермінанти» як взаємозамінні. Бітрядки показують заповненість електронами спін-орбіталей: $1$ на певній позиції означає, що відповідна спін-орбіталь зайнята, а $0$ означає, що вона вільна. Оскільки задачі електронної структури зберігають кількість частинок, лише фіксована кількість спін-орбіталей має бути зайнята. Молекула $N_2$ має $5$ спін-вгору ($\alpha$) та $5$ спін-вниз ($\beta$) електронів. Отже, будь-який бітрядок, що представляє $\alpha$- та $\beta$-орбіталі, повинен мати по п'ять $1$ для молекули $N_2$.

1.1 Квантова схема для генерування зразків: анзац LUCJ

У цьому уроці ми використаємо анзац локального унітарного зв'язаного кластера Жастрова (LUCJ) \[1\] для підготовки квантового стану та подальшого зразкування. Спочатку ми пояснимо різні будівельні блоки повного анзацу UCJ та наближення, використані в його локальній версії. Потім, за допомогою пакету ffsim, ми побудуємо анзац LUCJ та оптимізуємо його за допомогою транспілятора Qiskit для виконання на залізі.

Анзац UCJ має такий вигляд (для добутку $L$ шарів або повторень оператора UCJ):

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

де $\vert \Phi_{0} \rangle$ — опорний стан, який зазвичай береться як стан Хартрі–Фока (HF). Оскільки стан Хартрі–Фока визначається як такий, у якому зайняті орбіталі з найменшими номерами, підготовка стану HF потребує застосування вентилів X для встановлення кубітів, що відповідають зайнятим орбіталям, у одиницю. Наприклад, блок підготовки стану HF для 4 просторових орбіталей і 2 електронів зі спіном вгору та 2 зі спіном вниз може виглядати так: Одне повторення оператора UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ складається з діагональної кулонівської еволюції ($e^{iJ^{(\mu)}}$), обрамленої орбітальними обертаннями ($e^{K^{(\mu)}}$ та $e^{-K^{(\mu)}}$). Блоки орбітальних обертань працюють на одному спіновому виді ($\alpha$ (спін вгору)/$\beta$ (спін вниз)). Для кожного виду електронів орбітальне обертання складається з шару однокубітних вентилів $R_{z}$, за яким іде послідовність двокубітних вентилів обертань Ґівенса ($XX + YY$ вентилі).

Двокубітні вентилі діють на сусідніх спін-орбіталях (сусідніх кубітах) і тому можуть бути реалізовані на QPU IBM® без потреби у вентилях SWAP. $e^{iJ^{(\mu)}}$, також відомий як діагональний кулонівський оператор, складається з трьох блоків. Два з них працюють в одному спіновому секторі ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ та $e^{iJ_{\beta \beta}^{(\mu)}}$), а один — між двома спіновими секторами ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Усі блоки в $e^{iJ^{(\mu)}}$ складаються з вентилів числа-числа $U_{nn}(\phi)$ [1]. Вентиль $U_{nn}(\phi)$ може бути далі розкладений на вентиль $R_{ZZ}(\frac{\phi}{2})$, за яким ідуть два однокубітних вентилі $Rz(-\frac{\phi}{2})$, що діють на двох окремих кубітах.

Однаково-спінові компоненти ($J_{\alpha \alpha}$ і $J_{\beta \beta}$) мають вентилі $U_{nn}$ між усіма можливими парами кубітів. Проте, оскільки надпровідні QPU мають обмежену зв'язність, кубіти необхідно переставляти для реалізації вентилів між несусідніми кубітами.

Наприклад, розглянемо такий блок $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (або $e^{iJ_{\beta \beta}^{(\mu)}}$) для $N = 4$ просторових орбіталей. При лінійній зв'язності кубітів останні три вентилі не можна реалізувати безпосередньо, оскільки вони діють між несусідніми кубітами (наприклад, Q0 і Q2 не з'єднані напряму). Тому потрібні вентилі SWAP для того, щоб зробити їх сусідніми (на наступному рисунку показано приклад із $3$ вентилями SWAP). Далі, $J_{\alpha \beta}$ реалізує вентилі між однаково індексованими орбіталями з різних спінових секторів (наприклад, між $0\alpha$ і $0\beta$). Аналогічно, якщо кубіти фізично не суміжні на QPU, ці вентилі також вимагатимуть SWAP. З наведеного вище обговорення випливає, що анзац UCJ стикається з певними труднощами при виконанні на залізі, оскільки потребує вентилів SWAP через взаємодії несусідніх кубітів. Локальний варіант анзацу UCJ — LUCJ — вирішує цю проблему шляхом видалення деяких $U_{nn}$ з діагонального кулонівського оператора.

У блоках одного виду електронів ($J_{\alpha \alpha}$ і $J_{\beta \beta}$) у версії LUCJ ми зберігаємо лише ті вентилі $U_{nn}$, що сумісні зі зв'язністю найближчих сусідів, і видаляємо вентилі між несусідніми кубітами. На наступному рисунку показано блок LUCJ після видалення несусідніх вентилів. Далі, версія LUCJ блоку $J_{\alpha \beta}$, що працює між різними видами електронів, може мати різний вигляд залежно від топології пристрою.

Тут також версія LUCJ позбавляється несумісних вентилів. На рисунку нижче показані варіанти блоку $J_{\alpha \beta}$ для різних топологій кубітів: квадратна решітка, шестикутна, heavy-hex та лінійна.

• Квадратна: ми можемо мати вентилі $U_{nn}$ між усіма $\alpha$- і $\beta$-орбіталями без жодних SWAP, тому видалення будь-яких вентилів $U_{nn}$ не потрібне.
• Heavy-hex: взаємодії $\alpha$-$\beta$ зберігаються між кожною $4$-ю індексованою (наприклад, 0-ю, 4-ю і 8-ю) спін-орбіталлю і реалізуються через допоміжні кубіти, тобто між лінійними ланцюжками, що представляють $\alpha$- і $\beta$-орбіталі, потрібні допоміжні кубіти. Такий порядок потребує обмеженої кількості SWAP.
• Шестикутна: кожна друга орбіталь, наприклад 0-та, 2-га і 4-та індексовані, стає найближчими сусідами, коли $\alpha$ і $\beta$ розташовані в двох суміжних лінійних ланцюжках.
• Лінійна: з'єднані лише одна $\alpha$- і одна $\beta$-орбіталь, тобто блок $J_{\alpha \beta}$ матиме лише один вентиль. Хоча видалення вентилів з анзацу UCJ для побудови версії LUCJ робить його більш сумісним із залізом, анзац втрачає частину виразності. Тому при використанні анзацу LUCJ може знадобитися більше повторень ($L$) модифікованого оператора UCJ.

1.2 Ініціалізація анзацу LUCJ

LUCJ — це параметризований анзац, і перед виконанням на залізі потрібно ініціалізувати параметри. Один зі способів ініціалізації анзацу — використання амплітуд t1 і t2 класичного методу зв'язаних кластерів із одиничними та подвійними збудженнями (CCSD), де амплітуди t1 є коефіцієнтами операторів одиничних збуджень, а амплітуди t2 — для подвійних збуджень.

Зверни увагу: хоча ініціалізація анзацу LUCJ з амплітудами t1 і t2 дає непогані результати, параметри анзацу можуть потребувати подальшої оптимізації.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Побудова анзацу LUCJ за допомогою ffsim

Ми використаємо пакет ffsim для створення та ініціалізації анзацу з амплітудами t1 і t2, обчисленими вище. Оскільки наша молекула має замкнутооболонковий стан Хартрі–Фока, ми використаємо спін-збалансований варіант анзацу UCJ — UCJOpSpinBalanced.

Оскільки апаратне забезпечення IBM має топологію heavy-hex, ми застосуємо зиґзаґоподібний патерн, використаний у [1] і пояснений вище, для взаємодій між кубітами. У цьому патерні орбіталі (кубіти) з однаковим спіном з'єднані лінійною топологією (червоні та сині кола). Завдяки топології heavy-hex орбіталі різних спінів мають зв'язки між кожною 4-ю орбіталлю, тобто 0-ю, 4-ю, 8-ю і так далі (фіолетові кола).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

Анзац LUCJ з повторюваними шарами можна оптимізувати шляхом злиття деяких суміжних блоків. Розглянемо випадок для n_reps=2. Два блоки орбітальних обертань посередині можна об'єднати в один блок орбітального обертання. Пакет ffsim має менеджер проходів під назвою ffsim.qiskit.PRE_INIT, що оптимізує схему шляхом злиття таких суміжних блоків.

2. Оптимізація для цільового обладнання

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

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Далі рекомендуємо виконати такі кроки для оптимізації ансацу та забезпечення його сумісності з обладнанням.

• Вибери фізичні кубіти (initial_layout) на цільовому обладнанні, які відповідають зигзагоподібному шаблону (два лінійних ланцюжки з допоміжним кубітом між ними), описаному вище. Розміщення кубітів за таким шаблоном дає ефективну схему, сумісну з обладнанням, з меншою кількістю гейтів.
• Згенеруй поетапний менеджер проходів за допомогою функції generate_preset_pass_manager з Qiskit, передавши обраний backend і initial_layout.
• Встанови етап pre_init свого поетапного менеджера проходів рівним ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT містить проходи транспілятора Qiskit, які розкладають гейти на орбітальні обертання, а потім об'єднують їх, що призводить до меншої кількості гейтів у фінальній схемі.
• Запусти менеджер проходів на своїй схемі.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Виконання на цільовому обладнанні

Після оптимізації схеми для виконання на обладнанні ми готові запустити її на цільовому обладнанні та зібрати зразки для оцінки енергії основного стану. Оскільки у нас лише одна схема, ми скористаємось режимом виконання завдань Qiskit Runtime і виконаємо нашу схему.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Постобробка результатів

Частину постобробки в робочому процесі SQD можна підсумувати за допомогою такої діаграми.

Семплювання ансацу LUCJ у обчислювальному базисі генерує пул зашумлених конфігурацій $\tilde{\mathcal{\chi}}$, які використовуються в процедурі постобробки. Вона включає метод (деталі розглянуто нижче) відновлення конфігурацій для імовірнісного виправлення конфігурацій з неправильною кількістю електронів. Конфігурації лише з правильною кількістю електронів $\tilde{\mathcal{\chi}}_{R}$ потім підвибірково розподіляються на кілька батчів відповідно до частоти появи кожної унікальної конфігурації. Кожен батч зразків визначає підпростір ($\mathcal{S^{(k)}}$). Далі молекулярний гамільтоніан $H$ проєктується на підпростори:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

Кожен проєктований гамільтоніан $H_{\mathcal{S}^{(k)}}$ передається до власного розв'язувача, де він діагоналізується для обчислення власних значень і власних векторів з метою реконструкції власного стану. У цьому уроці ми проєктуємо та діагоналізуємо гамільтоніан за допомогою пакета qiskit-addon-sqd, який використовує метод Девідсона з PySCF для діагоналізації.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

Потім ми збираємо найменше власне значення (енергію) з батчів і обчислюємо середню орбітальну заповненість $\text{n}$. Інформація про середню заповненість використовується на кроці відновлення конфігурацій для імовірнісного виправлення зашумлених конфігурацій.

Далі ми детально пояснюємо самоузгоджений цикл відновлення конфігурацій і наводимо конкретні приклади коду для реалізації вищезазначених кроків з метою оцінки енергії основного стану гамільтоніана $N_2$.

4.1 Відновлення конфігурацій: огляд

Кожен біт у бітрядку (детермінанті Слейтера) представляє спін-орбіталь. Права половина бітрядка відповідає спін-вгору орбіталям, а ліва — спін-вниз. 1 означає, що орбіталь зайнята електроном, а 0 — що вона порожня. Нам заздалегідь відома правильна кількість частинок (як електронів зі спіном вгору, так і зі спіном вниз). Припустимо, що є детермінант $x$ з $N_x$ електронами (тобто у бітрядку є $N_x$ одиниць). Правильна кількість частинок — $N$. Якщо $N_x \neq N$, то ми знаємо, що бітрядок пошкоджений шумом. Самоузгоджена процедура конфігурацій намагається виправити бітрядок, імовірнісно перевертаючи $|N_x - N|$ бітів, використовуючи інформацію про середню орбітальну заповненість. Середня орбітальна заповненість ($n$) показує, наскільки ймовірно, що орбіталь зайнята електроном. Якщо $N_x < N$, то електронів замало і потрібно перевернути деякі $0$ на $1$, і навпаки.

Імовірність перевертання може бути $|x[i] - avg\_occupancy[i]|$ для $i$-ї спін-орбіталі. У [2] автори використовували зважену ймовірність перевертання за допомогою модифікованої функції ReLU.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

Тут $h$ визначає місце «кута» функції ReLU, а параметр $\delta$ — значення функції ReLU в цьому куті. При $\delta = 0$ $w$ стає справжньою функцією ReLU, а при $\delta >0$ — модифікованою ReLU. У статті автори використовували $\delta = 0.01$ і $h =$ кількість альфа- (або бета-) частинок / кількість альфа- (або бета-) спін-орбіталей $= N/M$ (коефіцієнт заповнення).

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

Розглянемо такий приклад для $N = 2$ і $x = |1000\rangle$ ($N_x = 1$). Потрібно перевернути один із нулів на одиницю для виправлення кількості частинок, і можливі варіанти — 1100, 1010 та 1001. На основі ймовірності перевертання буде вибрано один із варіантів як відновлену конфігурацію (або бітрядок з правильною кількістю частинок).

Припустимо, що на першій ітерації ми запускаємо два батчі, і оцінені основні стани з них такі:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

Використовуючи стани обчислювального базису та їхні амплітуди, ми можемо обчислити ймовірності електронної заповненості (скорочено заповненості) для кожної спін-орбіталі (кубіта) (зауваж, що ймовірність = |амплітуда|$^2$). Нижче ми табулюємо кубіт-по-кубіт заповненості для кожного бітрядка, що з'являється в оцінці основного стану, і обчислюємо загальну орбітальну заповненість для батча. Зауваж, що відповідно до конвенції впорядкування Qiskit, крайній правий біт відповідає кубіту-0 (Q0), а крайній лівий — Q3.

Заповненість (Batch0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

Заповненість (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

Заповненість (середнє по батчах)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (середнє)	0.49	0.51	0.35	0.65

Використовуючи обчислену вище середню орбітальну заповненість, ми можемо знайти ймовірності перевертання для різних орбіталей у конфігурації $x = \vert 1000 \rangle$. Оскільки орбіталь, представлена Q3, вже зайнята і її перевертати не потрібно, встановимо її p(flip) рівним $0$. Для решти орбіталей, які незайняті, ймовірність перевертання кожної дорівнює $\vert x[i] - \text{n}[i] \vert$. Разом із p(flip) ми також обчислюємо ваговий коефіцієнт ймовірності перевертання за допомогою модифікованої функції ReLU, описаної вище.

Ймовірність перевертання ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

	Q3	Q2	Q1	Q0
p(flip) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(flip))	0	0.03	0.007	0.31

Нарешті, використовуючи зважені ймовірності вище, ми можемо перевернути одну з незайнятих орбіталей Q2, Q1 або Q0. Виходячи з наведених значень, Q0 буде перевернуто найімовірніше, і можливою відновленою конфігурацією може бути $\vert \text{1001} \rangle$. Повний самоузгоджений процес відновлення конфігурацій можна підсумувати таким чином:

Перша ітерація: Припустимо, що бітрядки (конфігурації або детермінанти Слейтера), згенеровані квантовим комп'ютером, утворюють множину $\widetilde{\chi}$, яка включає конфігурації як з правильною ($\widetilde{\chi}_{correct}$), так і з неправильною ($\widetilde{\chi}_{incorrect}$) кількістю частинок у кожному спіновому секторі.

1. Конфігурації з ($\widetilde{\chi}_{correct}$) випадково вибірково відбираються для створення батчів $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ векторів для проєкції на підпростір. Кількість батчів і зразків у кожному батчі — параметри, що задаються користувачем. Чим більша кількість зразків у кожному батчі, тим більший вимір підпростору і тим обчислювально вимогливішою стає діагоналізація. З іншого боку, надто мала кількість зразків може пропустити опорні вектори основного стану і призвести до хибної оцінки.
2. Запустити розв'язувач власних станів (тобто проєкцію на підпростір і діагоналізацію) на батчах і отримати наближені власні стани $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. З наближених власних станів побудувати перше наближення для $n$.

Подальші ітерації:

1. Використовуючи $n$, виправити конфігурації з неправильною кількістю частинок у $\widetilde{\chi}_{incorrect}$. Назвемо їх $\widetilde{\chi}_{correct\_new}$. Тоді $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ утворює нову множину конфігурацій з правильною кількістю частинок.
2. $\widetilde{\chi}_{R}$ вибірково відбирається для створення батчів $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Розв'язувач власних станів запускається з новими батчами і генерує нові оцінки основних станів $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. З наближених власних станів побудувати уточнене наближення для $n$.
5. Якщо критерій зупинки не виконано, повернутись до кроку 2.1.

4.2 Оцінка основного стану

Спочатку перетворимо кількості на матрицю бітрядків та масив ймовірностей для постобробки.

Кожен рядок матриці представляє один унікальний бітрядок. Оскільки кубіти в Qiskit індексуються справа наліво в бітрядку, стовпець 0 відповідає кубіту N-1, а стовпець N-1 — кубіту 0, де N — кількість кубітів.

Альфа-орбіталі представлені в діапазоні індексів стовпців (N, N/2] (права половина), а бета-орбіталі — в діапазоні (N/2, 0] (ліва половина).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Є кілька параметрів, які контролює користувач і які важливі для цієї техніки:

• iterations: кількість самоузгоджених ітерацій відновлення конфігурацій
• n_batches: кількість батчів конфігурацій, що використовуються різними викликами розв'язувача власних станів
• samples_per_batch: кількість унікальних конфігурацій у кожному батчі
• max_davidson_cycles: максимальна кількість циклів Девідсона, що виконуються кожним розв'язувачем

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We post-select here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_davidson=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Обговорення результатів

Перший графік показує, що вже після кількох ітерацій ми отримуємо оцінку енергії основного стану з похибкою ~24 мГа (хімічна точність зазвичай прийнята на рівні 1 ккал/моль $\approx$ 1.6 мГа). Другий графік відображає середню заповненість кожного просторового орбіталя після фінальної ітерації. Можна побачити, що як електрони зі спіном угору, так і зі спіном донизу з великою імовірністю займають перші п'ять орбіталів у наших розв'язках.

Хоча отримана оцінка енергії основного стану є непоганою, вона не вкладається в межі хімічної точності ($\pm \approx 1.6$ мГа). Цю різницю можна пояснити малою розмірністю підпростору, використаного вище для проєкції та діагоналізації. Оскільки ми використали samples_per_batch=500, підпростір охоплює щонайбільше $500$ векторів, яким бракує векторів із опори основного стану. Збільшення параметра samples_per_batch має покращити точність ціною більших класичних обчислювальних ресурсів і часу виконання.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-6)
axs[0].axhline(
    y=chem_accuracy, color="#BF5700", linestyle="--", label="chemical accuracy"
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

Вправа для читача

Поступово збільшуй параметр samples_per_batch (наприклад, від $1000$ до $10000$ з кроком $1000$; в межах, дозволених пам'яттю твого комп'ютера) і порівнюй отримані оцінки енергії основного стану.

Посилання

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.