Langkau ke kandungan utama

Pepenjuruan kuantum berasaskan sampel terkumpul bagi Hamiltonian nuklear

Anggaran penggunaan: 2.5 minit pada pemproses Heron (NOTA: Ini hanyalah anggaran. Masa larian anda mungkin berbeza.)

Mencari versi Fortran?

Buku nota ini membentangkan pelaksanaan Python. Pelaksanaan Fortran terdapat dalam direktori pendamping Fortran repositori dokumentasi ini. Versi Python menambah langkah pemulihan konfigurasi-swakonsisten, yang tidak dilakukan oleh pemacu Fortran.

Hasil pembelajaran​

  • Belajar bagaimana Hamiltonian model-cengkerang nuklear, yang dijadualkan dalam asas orbital terganding-JJ, menjadi Hamiltonian qubit dalam skema-mm, di mana satu qubit ialah satu keadaan zarah-tunggal.

  • Bina ansatz pengujaan tetap, bukan-variasional yang sudutnya datang daripada teori gangguan peringkat-kedua, jadi tiada gelung pengoptimuman klasik.

  • Bandingkan pengujaan qubit dan fermionik dan ukur bagaimana pilihan tersebut menjejaskan kedalaman dua-qubit ensembel.

  • Jalankan pemulihan konfigurasi swakonsisten dengan qiskit-addon-sqd apabila kuantiti yang dikekalkan ialah bilangan nukleon, MJM_J, dan pariti dan bukannya bilangan elektron dan spin.

  • Aplikasikan satu aliran kerja daripada masalah 24-qubit yang boleh anda semak dengan tepat kepada masalah 40-qubit dengan hampir dua juta keadaan asas, melampaui keupayaan pepenjuruan-tepat tutorial ini.

Prasyarat​

Sebelum bermula, semak topik berikut:

Latar belakang​

Model cengkerang nuklear menganggap nukleus sebagai beberapa nukleon valens yang bergerak dalam satu set kecil orbital zarah-tunggal di atas teras lengai, berinteraksi melalui daya dua-jasad empirik yang dipadankan dengan spektrum yang diukur. Ia digunakan secara meluas dalam struktur nuklear tenaga-rendah. Kos pengkomputerannya ialah kombinatorik: asasnya ialah setiap cara mengagihkan proton dan neutron valens ke atas keadaan yang tersedia, dan pertumbuhan ini mengehadkan ruang model yang boleh dicapai oleh pepenjuruan tepat.

Pepenjuruan kuantum berasaskan sampel terkumpul (SQD terkumpul) [1] membahagikan masalah itu kepada dua. Satu litar kuantum digunakan hanya untuk mencadangkan keadaan asas yang penting. Ia diukur dalam asas pengkomputeran, dan setiap rentetan-bit yang diukur menamakan satu penentu Slater. Hamiltonian kemudiannya dibina dan dipenjurukan secara klasik dalam rentang penentu tersebut. Oleh kerana langkah klasik ialah pepenjuruan tepat dalam subruang, ia mengembalikan sempadan atas variasional bagi tenaga keadaan-asas sebenar, dan sempadan itu hanya boleh menurun apabila penentu ditambah.

Pembahagian kerja ini menjadikan kaedah tersebut bertoleransi-hingar, dengan satu batasan penting. Hingar mengubah yang mana penentu yang litar cadangkan. Ia tidak memasuki Hamiltonian klasik, jadi ia tidak boleh mengalihkan eigennilai satu-satu subruang: satu tembakan yang melanggar kuantiti yang dikekalkan dibuang atau dibaiki, dan satu tembakan yang bertahan ialah vektor asas yang sah tidak kira bagaimana ia dihasilkan. Oleh itu hingar merugikan kualiti subruang anda, bukan ketepatan, dan nombor yang anda laporkan tetap sempadan atas dalam kedua-dua keadaan.

Struktur nuklear menyediakan beberapa nombor kuantum tepat untuk penapisan sampel. Penentu fizikal mesti membawa bilangan proton valens yang betul dan bilangan neutron valens yang betul, unjuran momentum sudut jumlah MJM_J yang betul, dan pariti yang betul. Setiap satu boleh disemak dengan ujian integer ke atas rentetan-bit. Pecahan sampel yang ditolak bergantung kepada kekangan dan ruang model.

The 24-qubit sd-shell register: three proton orbitals and three neutron orbitals, each split into 2j+1 magnetic substates, one qubit per substate, with the reference determinant of neon-20 filled on the maximal magnetic substates of the 0d5/2 orbital.

Setiap qubit ialah satu keadaan zarah-tunggal skema-mm (n,ℓ,j,mj,tz)(n, \ell, j, m_j, t_z), dan ∣1⟩|1\rangle bermaksud terduduki. Daftar tersebut menggunakan susunan tetap: proton dahulu, kemudian neutron; dalam satu spesies, orbital dalam susunan fail; dalam satu orbital, mjm_j menurun. Oleh itu, dua separuh rentetan-bit ialah konfigurasi proton dan konfigurasi neutron. Ini ialah dwi-bahagian yang dijangkakan oleh alat pasca-pemprosesan SQD terkumpul.

Aliran kerja​

Workflow diagram: a reference determinant feeds an ensemble of shallow excitation circuits, which are sampled on a QPU to produce bitstrings; the bitstrings are repaired and post-selected on proton and neutron number, recombined into a product subspace where the magnetic projection and parity are imposed, and diagonalized to give a variational upper bound; average occupancies from the resulting eigenvector feed back into the next repair.

Dua peringkat dalam rajah mengendalikan simetri nuklear.

Pembaikan dan pasca-pemilihan mengendalikan sampel yang terjejas oleh hingar perkakasan. Dua bilangan nukleon separuh-daftar ialah berat Hamming, jadi qiskit-addon-sqd mengendalikannya secara langsung: recover_configurations membaiki rentetan-bit yang rosak dengan menyongsangkan bit yang paling tidak konsisten dengan anggaran semasa bagi kependudukan orbital purata, dan bukannya membuang tembakan tersebut.

Subruang hasil darab memperkenalkan MJM_J. Oleh kerana MJ=Mp+MnM_J = M_p + M_n menggandingkan kedua-dua separuh, ia bukan sifat mana-mana satu, jadi ia tidak boleh digunakan untuk menapis keseluruhan tembakan: satu rentetan-bit yang separuh protonnya dan separuh neutronnya masing-masing sah masih menyumbang dua separuh-konfigurasi baik walaupun MJM_J jumlahnya salah. Oleh itu subruang itu dirangkumi oleh setiap hasil darab konfigurasi proton yang disampel dengan konfigurasi neutron yang disampel, mengekalkan hasil darab yang jatuh dalam sektor MJM_J dan pariti sasaran. Ini ialah pembinaan subruang SQD terkumpul, dan ia bermaksud beberapa ribu rentetan-bit boleh merangkumi subruang yang jauh lebih besar daripada bilangan sampel.

Dua persamaan penentu​

Hamiltonian model-cengkerang ialah terma satu-jasad ditambah interaksi dua-jasad,

H=∑pεp ap†ap+14∑pqrs⟨pq∥rs⟩ ap†aq†asar,\begin{equation} \tag{1} H = \sum_{p} \varepsilon_p\, a_p^\dagger a_p + \tfrac{1}{4}\sum_{pqrs} \langle pq \| rs \rangle\, a_p^\dagger a_q^\dagger a_s a_r , \end{equation}

dengan p,q,r,sp,q,r,s melabelkan keadaan skema-mm dan tz=−1t_z = -1 untuk proton, +1+1 untuk neutron. Interaksi empirik seperti USDA [2] dan GXPF1 [3] dijadualkan bukan dalam skema-mm tetapi dalam asas terganding-JJ, sebagai unsur matriks ⟨ab;J∣V∣cd;J⟩\langle ab; J | V | cd; J \rangle antara keadaan dua-jasad antisimetri ternormal bagi orbital a,b,c,da,b,c,d. Memulihkan unsur skema-mm ialah penggandingan semula Clebsch-Gordan,

⟨pq∥rs⟩=1+δab1+δcd∑J⟨jpmp jqmq∣JM⟩⟨jrmr jsms∣JM⟩⟨ab;J∣V∣cd;J⟩,\begin{equation} \tag{2} \langle pq \| rs \rangle = \sqrt{1 + \delta_{ab}}\sqrt{1 + \delta_{cd}} \sum_{J} \langle j_p m_p\, j_q m_q | J M \rangle \langle j_r m_r\, j_s m_s | J M \rangle \langle ab; J | V | cd; J \rangle , \end{equation}

dengan faktor 1+δ\sqrt{1+\delta} membatalkan konvensyen normalisasi keadaan yang dijadualkan. Segala-galanya yang lain dalam tutorial ini dibina berdasarkan dua persamaan ini.

Tiga larian​

NukleusPetalaQubitAsas yang dibenarkan simetriBoleh disemak dengan tepat?
Berskala kecil20Ne^{20}\mathrm{Ne} (2p + 2n)sdsd24640Ya
Berskala besar44Ti^{44}\mathrm{Ti} (2p + 2n)pfpf404,000Ya
Berskala besar48Cr^{48}\mathrm{Cr} (4p + 4n)pfpf401,963,461Tidak

Larian skala-kecil ialah panduan langkah-demi-langkah. Kedua-dua larian skala-besar menggunakan daftar 40-qubit: yang pertama masih cukup kecil untuk dipenjurukan dengan tepat pada komputer riba, jadi anda boleh membandingkan keputusan perkakasan dengan rujukan tepat. Yang kedua melampaui keupayaan pepenjuruan-tepat tutorial ini.

Setiap larian di sini dilaksanakan pada QPU. Itu ialah pilihan yang dibuat untuk tutorial ini dan bukannya keperluan kaedah tersebut: kesemua tiga larian berkongsi satu backend dan satu bajet get supaya anda boleh membandingkan prestasinya pada saiz masalah yang berbeza.

Keperluan​

Pasang pakej berikut sebelum bermula:

  • Qiskit SDK v2.0 atau lebih baharu (pip install qiskit)

  • qiskit-ibm-runtime v0.40 or later (pip install qiskit-ibm-runtime)

  • Tambahan SQD v0.12 atau lebih baharu (pip install qiskit-addon-sqd)

  • NumPy, SciPy, dan Matplotlib (pip install numpy scipy matplotlib)

Anda juga memerlukan akaun IBM Quantum® dengan bukti kelayakan disimpan secara tempatan, dan akses kepada QPU dengan sekurang-kurangnya 40 qubit.

Tiada pakej simulator diperlukan, dan tiada fail data perlu dimuat turun. Dua fail interaksi yang digunakan tutorial ini dibenamkan dalam sel persediaan berikut dan ditulis ke direktori sementara apabila anda menjalankannya.

Persediaan​

Bahagian ini mengimport alat dan mentakrifkan pembantu model-cengkerang yang diperlukan oleh aliran kerja, dalam susunan aliran kerja menggunakannya. Fizik di sebalik setiap satu diterbitkan dalam Lampiran; komen menerangkan peranan setiap fungsi dalam aliran kerja.

Dua fail interaksi dibongkar dahulu. Kedua-duanya ialah set parameter yang diterbitkan, dibenamkan di sini supaya buku nota ini berdikari: usda.snt ialah Hamiltonian cengkerang-sdsd USDA [2] dan gxpf1.snt ialah Hamiltonian cengkerang-pfpf GXPF1 [3].

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-sqd qiskit-ibm-runtime scipy
from __future__ import annotations

import base64
import gzip
import itertools
import tempfile
from dataclasses import dataclass
from functools import lru_cache
from math import factorial, sqrt
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

from qiskit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.quantum_info import Operator, SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.counts import bit_array_to_arrays
from qiskit_addon_sqd.subsampling import postselect_by_hamming_right_and_left
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2

from scipy.linalg import eigh

_USDA_SNT_GZ = (
"H4sIAI5wqWoC/4WY224bRwyG7/UU9F0CONsZco5AUyBuC/QmQNG06KWhWGql1LYMy+kJefiOTrucGXJrQAgsf5nl8ecOr+CX"
"D9+9g/0K9pv1/T389rx7gPLrsH+Cz/vVctg+vgAsrgB+HeDdAD9t7zYv6+dr+DDA+z8223/X17B8XMFN+ev9+m+4ed799XjA"
"f9z8sy/4+s8BvoWYrsEERwbhFRqTXi+uCvOwW63vYf+0vFvDAgDo/AFIh8/hK7DH30354Omvb8o3V8c/vIUnMKtboK/wiGKF"
"+gnFEfVn9PQUe8bthNIRtftbsGfUtQbAGXUFfawM8K0BF9SP6MWA0BpwQcMRvRhwBSX86+fl3ct2d4jq4+eHa3hYv2x2q7fX"
"AJuPy+fb3cP69+Uh4luAT8djv95++eGV/fj6y6dvFudnmcX4lGOoBmttNOVncTL2FLs3NGT0l++ndJTv0cR8/t6NUanP8WMI"
"6nPC6C8/5wp+vnn//QLOP9ank3k2wRszkDn/Zwv15xiv41l2SDmjhuEFM4PJ0QkYcgzM4A1Jp1GDkcvtaWMAR9tosAa9huHk"
"AtoUBYwaF8hGarAxOywgxnoeN+Se2smF4H3QMPZQ653XMJpOi6GyrcLcZJtzLgoYcdsOyXKibV1Aom0fikJO0VLUMDeW3kAl"
"qQLWpd66qGE02eZ9lXqSPD3UWyl5DeOeGnINJnmaUggCVrtgh2DJNm3fVy8WF3LSMFa9hYga5iYX0IQsYK2nmTpJkoq81Acy"
"jLTTnPEaxqo3JGcbvZMeGnJVbw6mf2cUqcIsO80nFdOFq8JoSlY+Zb7FfNNZlnw2PqOE8SL3PtlSJhJGXGqsDRjbh4bmocWy"
"oiIxSxjv09JYgaxlGLZxA9EFFOImuIBt3EB04TK3Z5S8wizrrFCJQ4Xpgl9hzDYbu0LCNm6HzkJHweUkYazITalyF7NjGAnh"
"FZJFQniFZFHnqWgbdZ6qtoX50VZhvEJ81TJe7AUcQozH8yQMJ8ymUS07jLin6LXT3FRIhmwSTmtbJpEjr2AsvM7E1ra+ZaiM"
"NoxewZinBlPQTmPJyikn7TQ2630SbWuz4CLlwzuGhCHPqbUmRZQwpkiBDn1vsoQ55oKjbEvNNVgbXls0WvA0NKnPpku91Fkx"
"pahhXPDJtxVCQkBcMtJD27I8VJGAhWYYlVuZzc5e+jSIZWnLa57nOQ1iWZY+jSE3GP5/QIKi5E1Agqi9BbMZiUjEWNwChqKq"
"KTZY11ku5ca2Pqc4YKk1o2DstBjJC1htW2nAgEZ4aDdPvQmmmLeY4tWFt9Y3VAdlpW+oDcpa3zpM1jdUB2XV9agOyqrrcWZQ"
"sq7vMMffBseux5l5yjoLZ+Yp6yzU5mmt5DgzT5mSozZPoVIknJmnTJFQHZRV12OrlvpFQFbycoGKOWiYfl9QlXySGtSUvL9W"
"VJi/YG4grC4CmuCX6CZrKUsY64UyJolCQoZJqWeKhDOCzxQJNcGXbZMEX7EtzF+gvCI1npxjmKaW5Q06ewVjOY0u2+Y0kt5D"
"JhlETS2P+5CqyIMYkCI1NK50SO3TarSR2qfVaKOZPmXjg2beVJmS01zqp/CSmvoqvNSOD/kaS+qUoSGaTM2yz83fdjtMvsZy"
"xOv7N44Fff/GT5q5tXWYfB3jWND3b9yumUvK5SQnvK6w/VuHyfu3DpP3bx0m7984FvT9W4fJ+zceWq/v3zpM3r9xLOj7tw6T"
"928X053QgGz/1mHy/u3yZ8lTtn/jWND3b74NiDx2vRbeep56Lbz1oPRaeOv9G0dmxofXPK33bx0m798C1B9RuP4DPJ0WcLAa"
"AAA="
)

_GXPF1_SNT_GZ = (
"H4sIAI5wqWoC/4WcQY8dtw3H7/kU41sCeF8lUqKkQw9NgbaXoEHQQ26BkdiIEWdt2E6L9tNXb9+ORIn8PwdZILB/S3I4IsW/"
"ZiYvjr//+P3f4nG8ef/x+PDm4dOvr9+9++rF8d3l+Mf7x99fvTz+dTn++fnTH7/1//z28pfL8e3H9/95fHm8evzl+lffvf3f"
"H59e/fb26L9zHN//+t9Pl+OH1/++HH89JL88gkQO8esfvjm+phDomyvW//39/S+v3x2fPrz6+fVX/dfS889B4fpz/aN43P7p"
"f3Bw/ynH8fD0Zy+uf/fn48MR3vxU/kRXlp7Z+PzDiqUnNn74iW8sb3azYvm0m29s2uxGxabTbryx2cZ7nGzu7KOKV2y8g5Un"
"dsZbbLyDLafd53irjXew9bR7jffF8fbx8+uPr37+/Pb941dPv3P93ZH4W/If6kUohXCm+Jbmh3yR0jicybwl9CFeuFILy908"
"HtIlcrlZyCNJ2q6MdGi7ZVy4tlvHJS52Y71FnHp8D+HCIYTlQuYFhdsv0yVxzZB6vsxwaZwJUumkIjWB1POdDhcqlR2KVo85"
"SoDU8NgdNodiZevqMUqB1LCVWsO2ZMQlyYs+bXEFk4m5bMJJlVgSpGjaouhQvFKxUIbUjIvI85hWW4FNvubyDnMxMkFKRc8V"
"Uml65OhQW1x9fe3XOAtsZPXJ36TIX18SaoQUnysnZsHUXPfcKqTyaSu36FHbWg0tN0jxWIXEAVLzbqcqkMpj3T8lwlBp9RhD"
"qpAaHpnzni9nRXPh6lDbysmJBVLzDtU71KztnDxqW1+JhDaKnTtUa4DUzH00Hp0VXQMHSM3c16irlpf7GEcN1dogpfpEDpAa"
"WRUOBKl0UqkIQ2qs+5IYU6OvthbFodIaF6/9i51VeF0TLZSN2lfh9T729utQOqtXKi/dl/21mihhW+M+BpIKqblypDWHSntc"
"mSA1O3nMe+6dHp0DEaRmXCU0h9p7dOYCqVm1MehOnvy7nWphSKnuu+Q++X2Ca4uQGrYkxp3yMkFLL0ygavPovmR64dhFY5QM"
"KVK1XRxqX4UydlHCXU44hY1ydndKuUBK11BzqH1NNOKNMvt27PMXVUWBemSqCVIjrr4nYFtzTSQJkJp3O7se09p9YygCKVL7"
"UN4opx5DNnE5q5AlerZ2jyExpEYmggTtMYFrzC1Bak58jetGeSvHenSip9ROik0m5lqdsxzjrCYhTM1uImOHYZzVPvnWjXIm"
"0RoaKyqhjjlqiEE3uWoYHqopbT/KY0pDD+Xtx9d8lvI0n6U8zWcpT/NpQqDms5Sn+TRVoOazlKf5LOVpPk1VqPl05AI1n6U8"
"zaepAjWfpTzNp6kKNZ/OQoGaz1Ke5rOUp/k0VaHm03mvUPOd1+esL6X5LOVpPkt5ms9SnubTVIGaz1Ke5rOUp/ks5Wk+TVWo"
"+SzlaT69mgVqPk0VqPks5Wk+S3maT1MVaj5dPwVqPkt5mk9TFWo+S3ma77x/8z56ms9SnuazlKf5LOVpPkt5ms9SnubTVIWa"
"z1Ke5jupfRWumk9TBWo+S3maz1Ke5rOUp/k0VaHms5Sn+XTnLVDzWcrTfJqqUPNZytN85/1z7rbSfJbyNN/5t06fUJrPUp7m"
"OykvE7T0wgqqdmo+Mb3Q03yCO6bSfII7ptJ8gruc0nxi+pen+eROL5yaT3CXU5pPx1Sh5hNcj0rzCajHVfMJqMdV88mdesyu"
"xwo1n+B6VJpPcD0qzSe4HpXmE1yPSvMJrkelwERVR4Waz1Ke5hNTQ3THoxO90nzFZMLTfAVnVWm+grOqNF/BWVWaTxMVar6C"
"r1FpvgK6yar56vbja76onormu8/5ViqOrNblZH6l0NPAlRq9sKaIbaFnhis15tUuH/E1oieLK1XOTIgsz3ROap9XpQ/SJWqK"
"bCYcKm6K4tp9iUqqzbd1j4qO+g2tpTV6sll1qLipk6fnw32/isDWPSruPborwz7x9eakKbbry6Gi7RN9nrhuHr4tPXXsVLT9"
"vvTdKlUQ18i9Q0Wn+2ZOlYpva0yiDhXthNz69p4S+7ZG7h0q2jma+/BIIfm2RnU4VLQ9R6hG4iWryVurhorOlNbzEJr4tubk"
"bim6e25iqbFbJY4MKfREfaXGNXKK2aF2BSapD6Oy5162yd1S5Og0iq3k5tua1WEpsus+ZO4ZS76tqU4sRc7+mAtJyL6teR8t"
"RVbDUO9yXM2a2O6jQ5GjdCi3HKpva04dluK7J1uWGtFLjhFS6J2HlZpKhzOm0JsRK5XnmjD7Izv5av0eBd7zVb5IsVU6ffTN"
"oQXf1peoeveNjZWaNURFZ0L8fTsX1s9hwO4eeHkOA3Z3Zct6ZOccgPAM0O91hR4ZPEcWMAPkdBWQvsdZj72xZuTRs2U95qFY"
"a8FZHVQpy/Mhd55YbGmPew2VGB2PRvOVsNlib0orYbO1U8fi0cbF6uQhw7jUycMybYs/dSiKwGxyLB5tXEmdM5WM4lKnGLFm"
"FFeaPTruK8ehlEcbV55n2yIwrnlOXorAuPKsxyIwrjz7l3greptzSgve+tqo2lLbcu9MQ4oiPDMpjzaueR8rUUZxqaql/T46"
"k5WiCM9fyuNpy5uZpM9WkvfaxpS1Ndd96JKV2LflUdrW1idaKOT0HNMnTCdnT/OZvcOZv5RHG5fqAJFhXLO2WykwrtlzJkV4"
"llMebVxp7kOUMopr7h1dSmcUl9LuxesT2z6kPOq4ttkkk7sj1/1JWWWnHvfdfVKE50LlUV+jObPqs2hO+zViytpSJ+CFSyu+"
"LY+yttRM3gf3HH1bHmVtzR5NXZQXcI0epW1tdyh2d86a2O52F9vR6V8FUgTm1dWjjWuepZUgMK45k2cKMC6PIjD7rh6j+gpi"
"0Y904b675/n2h3u6QteTedkoc7qy2LIeaUxpkXNAHgeV53khozOYxZb1qDpTrvAa+TxBClECvEbPlvWY5o4cA7zGNDSfFIbX"
"6NmyHmd1SEnwGgfVJ6YKr9GzZT3KebdDZIYe5bzGQLFAj54t7XHXfDkmJ6v789r17e7iTEMrxXhmUh5tXHNmainCuObMVBOO"
"y6MYzEyrx9MW2VMf4cC9NWlbzpmVoqytucPkIO12SsZ4/lKUtTX7V6LGsfi2phadlLa1nX89naUHc437uYnQni/nlEzWN/S8"
"s7TVo41LvUsW9+pwqFBl79HsPYkVr/sK9GjjYnVWSzBfqjoowXx5FN+ZC6dHG5d6WybjfM15ouYE8+VRDObC1aOOa3v+WEJq"
"Tlzb2VCRELOtbYEU47NH5dHGpd5nagzj0s+HAozLoxjMq6tHnfutHmNjblz23G/1qChra3rMOT6f3jGefRVlbamZKfR1COKa"
"PWdS1tZcXy3mnMm3pd4IGpS1Nd9Uai0Wir6t+fxxUtrWPvu27M1MG8XJ1KMz+yqK8Vmt8mjjmn1CKmcUF88zhVIzimtOtZNi"
"fO6rPGpbW9XGLrg50G5rfwNhUueu7swmEtI55yQ0m/RpqI3JPd15ujVtWY9TSbcxpSU8dSShAj16tk4b3r6dmUlI2/L27Ulp"
"W7vKLFK3uByK4nI2VP19SFEJ7EOrRxvXyESr413YhPchLsssV/0dRlEJ70PKo7a1n/vWa1rLbsu8mTooa0up8tpyEd+Wep42"
"KG1rf6su5ebkvpqztLzl3qltRSVc28ojmTM+9F2TwNm3Qgp9/eSeKnaNXFtwqH3nk1KaJE05M6ZDOdNjv9NJIvm2Zr+3lFeP"
"EvqGnH1bM/eWMtNQvFCRHnzUlJlzXMqZYFKqHHLxbdFUFIbiu+8XWmo+50vL6R06l1u/PBO/asN8B5zAmdXTDtP7r8h+jeWL"
"lFOPocYumti39SWq3v0iTvxdtDw9hhmUq9N6PbaqvzJyqyNeMi0nIm51LLbojgIrOUKP+l2yBj16tghrqxib1gqg0pLknfIq"
"bdoirJqkMMNrnDNATPgaPVuE9VBv0J7HbU107V42W44eUhSjPrF4JKyHUokZxqXeTV+exAI9pCjGp/zKIwE9dO2+XSTX1vbc"
"b6pJUYT1UOCWWw2+rdmZJkVYD1Ef0Uoi39acySdFWA9RztLVjm9LPQ0cFAENc3unJjrVsWuY0mK297EYxbrbcs7vlUfCGib3"
"bQ/GNWffMijCGkZRDHr0w+KRsIbpnem6dnZbe9VOitTsW7e3I+e7Pgnv7tfHQxvl7O7KlvU4zu95fouUUNV21RSGFk2oHhdb"
"ZGbf+b4cXf8nRlnb8maASVlbc+UkSnSb0hKox4eFIjNHz7c/aHydkvC6T9Lydo3Je+u8ZSdf+xsu06O2ta2vPr70W1l2W9v6"
"UhSbPoG+kwYnSDVnTKGvqYvfv3JuDVLom2vvpOZpvh9fizHuX71jVp79vqDcG8rrJolLnecTxe8TgMLfBTDoJteOGUVntYJr"
"5PHOVsKnUaEuJw8oE9MWA5V5e2tA9wmQrw5xQB6RrernS33VkEB1PCmKW5f7P6zudeW7TwAA"
)

DATA = Path(tempfile.mkdtemp(prefix="nuclear_sqd_"))
for name, blob in (("usda.snt", _USDA_SNT_GZ), ("gxpf1.snt", _GXPF1_SNT_GZ)):
(DATA / name).write_bytes(gzip.decompress(base64.b64decode(blob)))

if not (DATA / name).is_file():
raise RuntimeError(f"{name} did not unpack to {DATA}")

Ruang model dan daftar qubit​

Fail .snt menyimpan ruang model, tenaga zarah-tunggal, dan unsur matriks dua-jasad terganding-JJ. Bagi interaksi bergantung-jisim yang digunakan di sini, medan ketiga dan keempat pengepala dua-jasad menentukan jisim rujukan ArefA_{\mathrm{ref}} di mana interaksi tersebut dipadankan dan eksponen kebergantungan jisimnya. Kedua-dua fail membawa eksponen −0.3-0.3, dengan Aref=18A_{\mathrm{ref}} = 18 untuk USDA dan 4242 untuk GXPF1, jadi unsur matriks yang dijadualkan mesti diskala semula dengan (A/Aref)−0.3(A/A_{\mathrm{ref}})^{-0.3} untuk nukleus yang dikira [2], [3]. Tenaga zarah-tunggal tidak diskala semula. Melangkau langkah ini mengubah tenaga korelasi sebanyak beberapa peratus.

Tenaga yang berikut ialah tenaga valens, diukur daripada teras lengai; ia bukan tenaga pemisahan eksperimen.

@dataclass(frozen=True)
class Orbital:
idx: int
n: int
ell: int
j2: int
tz: int # j2 = 2j; tz = -1 proton, +1 neutron

@dataclass(frozen=True)
class SPState:
"""One m-scheme single-particle state, i.e. one qubit."""

orb: int
j2: int
mj2: int
tz: int
ell: int
spe: float # mj2 = 2 * m_j

@dataclass
class ModelSpace:
orbitals: list
spes: dict
tbmes: dict
core_z: int
core_n: int
mass_number: int
a_ref: int
mass_exponent: float
mass_factor: float

def read_snt(path, n_protons, n_neutrons):
"""Parse a .snt interaction file, applying its mass dependence for this nucleus.

The two-body header line is ``n_tbme method A_ref exponent``. When ``method`` is 1 the
tabulated matrix elements are rescaled by ``(A / A_ref) ** exponent``, where A is the mass
number of the whole nucleus -- the core plus the valence nucleons. A is derived from the
file's own core numbers rather than passed in, so it cannot silently disagree with the
valence counts the rest of the workflow uses. Single-particle energies are not rescaled.
"""
rows = [
ln.split("!")[0].split() for ln in Path(path).read_text().splitlines()
]
rows = iter([r for r in rows if r])

n_p_orb, n_n_orb, core_z, core_n = (int(x) for x in next(rows)[:4])
orbitals = [
Orbital(*(int(x) for x in next(rows)[:5]))
for _ in range(n_p_orb + n_n_orb)
]

spes = {}
for _ in range(int(next(rows)[0])): # "i i <i|H(1b)|i>"
field = next(rows)
spes[int(field[0])] = float(field[2])

n_tbme, method, a_ref, exponent = next(rows)[:4]
n_tbme, method, a_ref, exponent = (
int(n_tbme),
int(method),
int(a_ref),
float(exponent),
)
mass_number = core_z + core_n + n_protons + n_neutrons
factor = (mass_number / a_ref) ** exponent if method == 1 else 1.0

tbmes = {}
for _ in range(n_tbme): # "a b c d J value"
field = next(rows)
tbmes[tuple(int(x) for x in field[:5])] = float(field[5]) * factor

return ModelSpace(
orbitals,
spes,
tbmes,
core_z,
core_n,
mass_number,
a_ref,
exponent,
factor,
)

def m_scheme_states(ms):
"""The qubit register: protons then neutrons, orbitals in file order, m_j descending."""
return [
SPState(o.idx, o.j2, m2, o.tz, o.ell, ms.spes[o.idx])
for tz in (-1, +1)
for o in ms.orbitals
if o.tz == tz
for m2 in range(o.j2, -o.j2 - 1, -2)
]

Penggandingan semula Clebsch-Gordan​

Persamaan (2) memerlukan pekali Clebsch-Gordan bagi momentum sudut separuh-integer. Setiap hujah dilepasi sebagai dua kali nilai fizikalnya, jadi j=5/2j = 5/2 memasuki sebagai 5 dan aritmetik kekal tepat.

Interaction.v_ms mengendalikan carian unsur matriks interaksi. Fail .snt menyimpan setiap unsur matriks sekali, jadi carian mungkin memerlukan fasa pertukaran-pasangan antisimetri −(−1)ja+jb−J-(-1)^{j_a + j_b - J} pada mana-mana belah, dan bra serta ket mungkin disimpan dalam mana-mana susunan.

@lru_cache(maxsize=None)
def clebsch_gordan(j1_2, j2_2, J_2, m1_2, m2_2, M_2):
"""<j1 m1 j2 m2 | J M>. Every argument is twice its physical value."""
if m1_2 + m2_2 != M_2 or not abs(j1_2 - j2_2) <= J_2 <= j1_2 + j2_2:
return 0.0
if abs(m1_2) > j1_2 or abs(m2_2) > j2_2 or abs(M_2) > J_2:
return 0.0
if (j1_2 + j2_2 - J_2) % 2 or (j1_2 - m1_2) % 2 or (j2_2 - m2_2) % 2:
return 0.0

f, half = factorial, lambda x: x // 2
prefactor = sqrt(
(J_2 + 1)
* f(half(j1_2 + j2_2 - J_2))
* f(half(j1_2 - j2_2 + J_2))
* f(half(-j1_2 + j2_2 + J_2))
/ f(half(j1_2 + j2_2 + J_2) + 1)
* f(half(J_2 + M_2))
* f(half(J_2 - M_2))
* f(half(j1_2 - m1_2))
* f(half(j1_2 + m1_2))
* f(half(j2_2 - m2_2))
* f(half(j2_2 + m2_2))
)
total = 0.0
for k in range(half(j1_2 + j2_2 - J_2) + 1):
d = [
half(j1_2 + j2_2 - J_2) - k,
half(j1_2 - m1_2) - k,
half(j2_2 + m2_2) - k,
half(J_2 - j2_2 + m1_2) + k,
half(J_2 - j1_2 - m2_2) + k,
]
if all(x >= 0 for x in d):
total += (-1) ** k / (
f(k) * f(d[0]) * f(d[1]) * f(d[2]) * f(d[3]) * f(d[4])
)
return prefactor * total

class Interaction:
"""Antisymmetrized m-scheme two-body matrix elements <pq||rs>, per Eq. (2)."""

def __init__(self, model_space, sp):
self.ms, self.sp, self._cache = model_space, sp, {}

def _tbme(self, oa, ob, oc, od, J, j_ab_2, j_cd_2):
"""<oa ob; J|V|oc od; J>, allowing for how the file happens to order each pair."""
table = self.ms.tbmes
# |ba; J> = -(-1)^(j_a + j_b - J) |ab; J> for a normalized antisymmetrized pair;
# dropping the leading minus makes v_ms symmetric instead of antisymmetric, and
# the Hamiltonian then fails the rotational-invariance check in Step 1.
phase_ab = -1.0 if (j_ab_2 // 2 - J) % 2 == 0 else 1.0
phase_cd = -1.0 if (j_cd_2 // 2 - J) % 2 == 0 else 1.0
for keys, phase in (
(((oa, ob, oc, od), (oc, od, oa, ob)), 1.0),
(((ob, oa, oc, od), (oc, od, ob, oa)), phase_ab),
(((oa, ob, od, oc), (od, oc, oa, ob)), phase_cd),
(((ob, oa, od, oc), (od, oc, ob, oa)), phase_ab * phase_cd),
):
for key in keys:
value = table.get(key + (J,))
if value is not None:
return value * phase
return 0.0

def v_ms(self, p, q, r, s):
"""<pq||rs>, zero unless M_J and charge are conserved."""
cached = self._cache.get((p, q, r, s))
if cached is not None:
return cached

P, Q, R, S = (self.sp[i] for i in (p, q, r, s))
value = 0.0
if P.mj2 + Q.mj2 == R.mj2 + S.mj2 and P.tz + Q.tz == R.tz + S.tz:
M = P.mj2 + Q.mj2
# sqrt(1 + delta): undo the normalization of the tabulated pair states
c12 = sqrt(2.0) if (P.tz == Q.tz and P.orb == Q.orb) else 1.0
c34 = sqrt(2.0) if (R.tz == S.tz and R.orb == S.orb) else 1.0
for J2 in range(
max(abs(P.j2 - Q.j2), abs(R.j2 - S.j2)),
min(P.j2 + Q.j2, R.j2 + S.j2) + 1,
2,
):
cg_bra = clebsch_gordan(P.j2, Q.j2, J2, P.mj2, Q.mj2, M)
cg_ket = clebsch_gordan(R.j2, S.j2, J2, R.mj2, S.mj2, M)
if abs(cg_bra) < 1e-12 or abs(cg_ket) < 1e-12:
continue
value += (
c12
* c34
* cg_bra
* cg_ket
* self._tbme(
P.orb,
Q.orb,
R.orb,
S.orb,
J2 // 2,
P.j2 + Q.j2,
R.j2 + S.j2,
)
)

self._cache[(p, q, r, s)] = value
return value

Elemen matriks dan ujian simetri​

Penentu ialah tuple tersusun bagi indeks qubit yang diduduki. Dua penentu yang berbeza pada lebih daripada dua keadaan terduduki mempunyai unsur matriks lenyap; jika tidak, peraturan Slater-Condon memberikan hasil tambah pendek ke atas interaksi, didarab dengan tanda fermionik yang mengira berapa banyak keadaan terduduki terletak antara operator dalam susunan daftar yang tetap.

symmetry_allowed ialah ujian integer di mana keempat-empat nombor kuantum tepat direduksikan kepadanya. Ia digunakan kedua-dua untuk menapis sampel dan untuk menyenaraikan asas tepat bagi larian yang cukup kecil untuk disemak.

def matrix_element(inter, det_a, det_b):
"""<A|H|B> for two determinants, each a sorted tuple of occupied qubit indices."""
set_a, set_b = set(det_a), set(det_b)
out_a, out_b = sorted(set_a - set_b), sorted(set_b - set_a)
if len(out_a) != len(out_b) or len(out_a) > 2:
return 0.0

if not out_a: # diagonal: one-body plus two-body
return sum(inter.sp[i].spe for i in det_a) + sum(
inter.v_ms(i, j, i, j)
for i, j in itertools.combinations(det_a, 2)
)

if len(out_a) == 1: # one state moves, p -> q
p, q = out_a[0], out_b[0]
crossings = sum(1 for k in set_a if min(p, q) < k < max(p, q))
return (-1.0) ** crossings * sum(
inter.v_ms(p, j, q, j) for j in det_a if j not in (p, q)
)

(p, r), (q, s) = out_a, out_b # two states move
crossings = sum(1 for k in set_a if p < k < r) + sum(
1 for k in set_b if q < k < s
)
return (-1.0) ** crossings * inter.v_ms(p, r, q, s)

def subspace_hamiltonian(inter, dets):
"""Dense real-symmetric H projected onto the span of `dets`."""
H = np.zeros((len(dets), len(dets)))
for a, det_a in enumerate(dets):
H[a, a] = matrix_element(inter, det_a, det_a)
for b in range(a + 1, len(dets)):
H[a, b] = H[b, a] = matrix_element(inter, det_a, dets[b])
return H

def ground_state(inter, dets):
"""Lowest eigenvalue and eigenvector of H over `dets`."""
values, vectors = np.linalg.eigh(subspace_hamiltonian(inter, dets))
return values[0], vectors[:, 0]

def symmetry_allowed(
sp, det, n_protons, n_neutrons, mj2_target=0, parity_target=0
):
"""The four exact shell-model quantum numbers, as integer tests on one determinant."""
n_p = sum(1 for i in det if sp[i].tz == -1)
return (
n_p == n_protons
and len(det) - n_p == n_neutrons
and sum(sp[i].mj2 for i in det) == mj2_target
and sum(sp[i].ell for i in det) % 2 == parity_target
)

def full_basis(sp, n_protons, n_neutrons, **targets):
"""Every symmetry-allowed determinant. Only tractable for small model spaces."""
protons = [i for i, s in enumerate(sp) if s.tz == -1]
neutrons = [i for i, s in enumerate(sp) if s.tz == +1]
return [
p + n
for p in itertools.combinations(protons, n_protons)
for n in itertools.combinations(neutrons, n_neutrons)
if symmetry_allowed(sp, p + n, n_protons, n_neutrons, **targets)
]

def count_basis(sp, n_protons, n_neutrons, mj2_target=0, parity_target=0):
"""How many determinants `full_basis` would return, without enumerating them.

A dynamic program over (occupied count, sum of 2*m_j, parity) per species. This stays
cheap when the basis itself is far too large to build, which is how the largest run below
can report the size of the space it is sampling from.
"""

def species(states, k):
table = {(0, 0, 0): 1}
for s in states:
for key, value in list(table.items()):
count, m_sum, parity = key
if count < k:
nxt = (count + 1, m_sum + s.mj2, (parity + s.ell) % 2)
table[nxt] = table.get(nxt, 0) + value
totals = {}
for (count, m_sum, parity), value in table.items():
if count == k:
totals[(m_sum, parity)] = (
totals.get((m_sum, parity), 0) + value
)
return totals

left = species([s for s in sp if s.tz == -1], n_protons)
right = species([s for s in sp if s.tz == +1], n_neutrons)
return sum(
a * b
for (mp, pp), a in left.items()
for (mn, pn), b in right.items()
if mp + mn == mj2_target and (pp + pn) % 2 == parity_target
)

Penentu rujukan​

Ansatz dibina di atas satu penentu, jadi penentu itu sepatutnya yang terbaik yang tersedia. Mengisi tenaga zarah-tunggal terendah mengabaikan interaksi dua-jasad. Dalam ruang model ini, pilihan itu memberikan tenaga 1–2 MeV di atas penentu bertenaga-terendah.

Menghadkan kepada pengisian yang terdiri daripada pasangan (+mj,−mj)(+m_j, -m_j) terbalik-masa memaksa MJ=0M_J = 0 dengan tepat dan hanya meninggalkan (npairsk)\binom{n_{\mathrm{pairs}}}{k} calon setiap spesies (paling banyak beberapa ribu), jadi yang terbaik boleh dijumpai dengan mencari semuanya pada diagonal penuh ⟨Φ∣H∣Φ⟩\langle \Phi | H | \Phi \rangle. Untuk keadaan seri, keutamaan diberikan kepada pasangan yang paling kuat sejajar, di mana daya pemasangan J=0J = 0 paling kuat. Dalam setiap kes dalam tutorial ini yang boleh disemak terhadap penyenaraian penuh, carian mengembalikan penentu diagonal-terendah global, yang juga merupakan komponen tunggal terbesar keadaan-asas tepat.

def reference_determinant(sp, inter, n_protons, n_neutrons):
"""Lowest-diagonal determinant built from time-reversed (+m_j, -m_j) orbital pairs."""
if n_protons % 2 or n_neutrons % 2:
raise ValueError(
"an odd valence count has no time-reversed paired reference at M_J = 0"
)

def species_pairs(tz):
return [
(
q,
next(
p
for p, t in enumerate(sp)
if t.tz == tz and t.orb == s.orb and t.mj2 == -s.mj2
),
)
for q, s in enumerate(sp)
if s.tz == tz and s.mj2 > 0
]

best = None
for chosen_p in itertools.combinations(species_pairs(-1), n_protons // 2):
protons = tuple(q for pair in chosen_p for q in pair)
for chosen_n in itertools.combinations(
species_pairs(+1), n_neutrons // 2
):
det = tuple(
sorted(protons + tuple(q for pair in chosen_n for q in pair))
)
# break ties toward the most aligned pairs, where J = 0 pairing is strongest
score = (
matrix_element(inter, det, det),
-sum(abs(sp[q].mj2) for q in det),
)
if best is None or score < best[0]:
best = (score, det)
return best[1]

Kumpulan pengujaan dan kedudukan perturbatifnya​

Korelasi dibawa oleh pengujaan dua-zarah–dua-lubang (2p2h2p2h) daripada rujukan. Dua peraturan pemilihan mengurangkan kumpulan sebelum mana-mana litar dibina: pengujaan mesti mengekalkan MJM_J, dan pasangan lubang serta pasangan zarah mesti boleh bergandingan kepada JJ jumlah yang sama, yang merupakan ketaksamaan segitiga.

Pengujaan yang tinggal disusun mengikut skor peringkat-kedua Epstein-Nesbet bagi interaksi konfigurasi terpilih [4],

sα=∣⟨Φref∣H∣α⟩∣2∣Δα∣,Δα=Href,ref−Hαα,\begin{equation} \tag{3} s_\alpha = \frac{|\langle \Phi_{\mathrm{ref}} | H | \alpha \rangle|^2}{|\Delta_\alpha|}, \qquad \Delta_\alpha = H_{\mathrm{ref},\mathrm{ref}} - H_{\alpha\alpha}, \end{equation}

yang menganggarkan berapa banyak tenaga korelasi yang dibawa oleh setiap pengujaan. Dua nombor yang sama menetapkan sudut litar: dengan V=⟨Φref∣H∣α⟩V = \langle \Phi_{\mathrm{ref}} | H | \alpha \rangle, amplitud peringkat-pertama ialah tα=V/Δαt_\alpha = V / \Delta_\alpha. Lampiran menerangkan mengapa amplitud peringkat-pertama ialah pilihan yang digunakan dalam tutorial ini dan bukannya sudut dua-aras tepat.

def excitation_pool(sp, occ):
"""2p2h quadruples (h1, h2, v1, v2): same-species pairs, then proton-neutron pairs."""
holes = {tz: [i for i in occ if sp[i].tz == tz] for tz in (-1, +1)}
virtuals = {
tz: [i for i, s in enumerate(sp) if s.tz == tz and i not in occ]
for tz in (-1, +1)
}
pool = [
(h1, h2, v1, v2)
for tz in (-1, +1)
for h1, h2 in itertools.combinations(holes[tz], 2)
for v1, v2 in itertools.combinations(virtuals[tz], 2)
]
pool += [
(h1, h2, v1, v2)
for h1 in holes[-1]
for h2 in holes[+1]
for v1 in virtuals[-1]
for v2 in virtuals[+1]
]
return pool

def conserves_symmetry(sp, op):
"""Keeps M_J, and the hole and particle pairs share a reachable total J."""
h1, h2, v1, v2 = op
if sp[v1].mj2 + sp[v2].mj2 != sp[h1].mj2 + sp[h2].mj2:
return False
return max(abs(sp[v1].j2 - sp[v2].j2), abs(sp[h1].j2 - sp[h2].j2)) <= min(
sp[v1].j2 + sp[v2].j2, sp[h1].j2 + sp[h2].j2
)

def en_denominator(inter, occ, holes, virtuals, floor=0.1):
"""Epstein-Nesbet gap: bare gap, spectator rearrangement, and the pair's own term."""
gap = sum(inter.sp[h].spe for h in holes) - sum(
inter.sp[v].spe for v in virtuals
)
for k in occ:
if k in holes:
continue
gap += sum(inter.v_ms(h, k, h, k) for h in holes)
gap -= sum(inter.v_ms(v, k, v, k) for v in virtuals)
gap += inter.v_ms(holes[0], holes[1], holes[0], holes[1])
gap -= inter.v_ms(virtuals[0], virtuals[1], virtuals[0], virtuals[1])
return gap if abs(gap) >= floor else (floor if gap >= 0 else -floor)

def rank_pool(inter, occ, pool):
"""Sort by descending PT2 score; return (operator, coupling, first-order amplitude)."""
ranked = []
for op in pool:
h1, h2, v1, v2 = op
coupling = inter.v_ms(v1, v2, h1, h2)
gap = en_denominator(inter, occ, (h1, h2), (v1, v2))
ranked.append((coupling**2 / abs(gap), op, coupling, coupling / gap))
ranked.sort(key=lambda row: (-row[0], row[1])) # deterministic on ties
return [
(op, coupling, amplitude) for _, op, coupling, amplitude in ranked
]

Blok pengujaan qubit​

Di bawah pemetaan Jordan-Wigner, operator pengujaan 2p2h2p2h yang mengekalkan zarah menjadi hasil tambah lapan rentetan Pauli, setiap satu membawa rentetan operator ZZ antara indeks paling luar. Rentetan ZZ menguatkuasakan antisimetri fermionik, dan ia mahal: pengujaan proton-neutron merentangi sempadan antara dua separuh daftar dan merangkumi rentetan pariti merentangi sempadan tersebut.

Menjatuhkan rentetan ZZ memberikan operator pengujaan-qubit Yordanov et al. [5]. Keadaan yang disediakan oleh operator ini mempunyai amplitud yang berbeza, tetapi ia menghubungkan pasangan penentu yang sama tepat, jadi set penentu yang boleh dicapai litar tidak berubah. SQD terkumpul menggunakan penentu ini untuk pepenjuruan klasik. Langkah 2 membandingkan sokongan kedua-dua pembinaan tersebut dan mengukur kos perkakasannya.

Membina bentuk Pauli daripada aj†=12(Xj−iYj)⊗Z<ja_j^\dagger = \tfrac{1}{2}(X_j - i Y_j) \otimes Z_{<j}, dengan rentetan ZZ sebagai pilihan, mengekalkan kedua-dua pembinaan hanya berbeza satu bendera. Kesemua lapan terma satu penjana berkomut, jadi satu langkah PauliEvolutionGate ialah eksponen tepat dan bukan penghampiran Trotter kepadanya.

def _ladder(num_qubits, q, dagger, parity):
"""Pauli form of a_q or a_q^dagger. `parity` toggles the Jordan-Wigner Z string."""
prefix = (
["Z"] * q + ["I"] * (num_qubits - q) if parity else ["I"] * num_qubits
)
x_part, y_part = list(prefix), list(prefix)
x_part[q], y_part[q] = "X", "Y"
return SparsePauliOp(
["".join(reversed(x_part)), "".join(reversed(y_part))],
coeffs=[0.5, 0.5 * (-1j if dagger else 1j)],
)

def excitation_generator(num_qubits, op, parity=False):
"""Hermitian H with exp(-i theta H) = exp(theta (T - T^dagger)) for T = a+ a+ a a."""
h1, h2, v1, v2 = op
T = SparsePauliOp("I" * num_qubits)
for q, dagger in ((v1, True), (v2, True), (h2, False), (h1, False)):
T = (T @ _ladder(num_qubits, q, dagger, parity)).simplify()
return (1j * (T - T.adjoint())).simplify()

def excitation_block(op, theta, parity=False):
"""(window, circuit) for one excitation.

A qubit excitation touches only its four qubits. A fermionic excitation also carries Z
operators on every qubit between the outermost indices, so its window is the whole span --
which is exactly where its extra cost comes from.
"""
window = list(range(min(op), max(op) + 1)) if parity else sorted(op)
local = tuple(window.index(i) for i in op)
generator = excitation_generator(len(window), local, parity=parity)
return window, PauliEvolutionGate(generator, time=theta).definition

def excitation_ansatz(
num_qubits, occ, operators, amplitudes, measure=True, parity=False
):
"""X gates for the reference determinant, then one evolution block per excitation."""
qc = QuantumCircuit(num_qubits)
for q in occ:
qc.x(q)
for op, theta in zip(operators, amplitudes):
if abs(theta) < 1e-12:
continue
window, block = excitation_block(op, theta, parity=parity)
qc.compose(block, qubits=window, inplace=True)
if measure:
qc.measure_all()
return qc

Bajet kedalaman dan ensembel litar​

Satu litar dalam yang mengandungi setiap pengujaan tersusun boleh melebihi masa koheren perkakasan. Menyebarkan kumpulan ke atas ensembel litar cetek dan mengumpulkan tembakan mereka ke dalam satu set penentu mengubah Langkah 2 menjadi masalah pembungkusan: setiap pengujaan mempunyai kos yang diukur, setiap litar mempunyai bajet, dan persoalannya ialah berapa banyak kumpulan tersusun yang muat.

Bajet diukur dalam kedalaman dua-qubit (lapisan get dua-qubit pada laluan genting) dan bukan dalam kiraan get mentah, kerana kedalaman menentukan tempoh litar dan oleh itu berapa banyak koheren peranti yang dibelanjakan. Jumlah kiraan dilaporkan bersama, kerana itu ialah proksi yang lebih baik bagi ralat get terkumpul; kedua-duanya menjawab soalan yang berbeza dan tiada satu pun menggantikan yang lain.

Kedua-dua kuantiti diekstrak mengikut ariti: satu arahan yang bertindak ke atas tepat dua qubit, tidak kira apa backend itu memanggil get penggandingannya. Padanan mengikut nama get sebaliknya boleh mengembalikan sifar untuk set asas yang tidak dikenali, secara salah meletakkan keseluruhan kumpulan dalam satu litar tanpa melebihi bajet yang dikira.

Mengisi litar mana sahaja yang paling kosong pada masa itu, mengikut susunan kedudukan, mengekalkan setiap litar hampir dengan bajet. Kos diukur pada sasaran backend sebenar, satu pengujaan pada satu masa, kerana kos yang dibaca daripada litar abstrak bukan kos yang dihasilkan oleh transpiler.

DIRECTIVES = ("barrier", "delay")

def is_two_qubit(instruction):
"""True for an operation on exactly two qubits, excluding directives.

Selecting by arity rather than by gate name keeps this correct on any backend, whatever its
two-qubit basis gate happens to be called -- cz on today's Heron devices, ecr on Eagle, or
something newer tomorrow. A gate-name allow-list silently returns zero on anything it has
not heard of, which would collapse the whole pool into one circuit and pass every budget
check. Barriers are excluded because a barrier spanning two qubits is not a gate.
"""
return (
len(instruction.qubits) == 2
and instruction.operation.name not in DIRECTIVES
)

def two_qubit_count(qc):
"""How many two-qubit gates the circuit contains: the accumulated-gate-error proxy."""
return sum(1 for instruction in qc.data if is_two_qubit(instruction))

def two_qubit_depth(qc):
"""Layers of two-qubit gates on the critical path: the duration and decoherence proxy.

This is what the budget is measured in. Two gates on disjoint qubit pairs run in the same
layer, so depth tracks how long the circuit takes -- and therefore how much coherence it
spends -- while the count above tracks how much gate error it accumulates. Both are
reported; only depth is budgeted.
"""
return qc.depth(filter_function=is_two_qubit)

def excitation_costs(num_qubits, ranked, pm, parity=False):
"""Transpiled two-qubit depth of each excitation on its own."""
return [
two_qubit_depth(
pm.run(
excitation_ansatz(
num_qubits, (), [op], [amp], measure=False, parity=parity
)
)
)
for op, _, amp in ranked
]

def pack_ensemble(
num_qubits, occ, ranked, costs, budget, n_circuits, parity=False
):
"""Fill n_circuits in rank order, always adding to whichever is currently emptiest."""
bins, loads = [[] for _ in range(n_circuits)], [0] * n_circuits
for (op, _, amplitude), cost in zip(ranked, costs):
emptiest = min(range(n_circuits), key=lambda b: loads[b])
if loads[emptiest] + cost > budget:
break # every circuit is full
bins[emptiest].append((op, amplitude))
loads[emptiest] += cost
circuits = [
excitation_ansatz(
num_qubits,
occ,
[o for o, _ in b],
[a for _, a in b],
parity=parity,
)
for b in bins
]
return circuits, bins

def pack_to_budget(
num_qubits, occ, ranked, costs, budget, n_circuits, pm, attempts=6
):
"""Pack, transpile, and shrink the target until the assembled circuits really fit.

Costs are measured one excitation at a time, but excitations that share qubits neither add
nor parallelize cleanly once the transpiler routes them together, so the assembled depth is
not the sum of its measured parts. This loop closes that gap against the real transpiler,
and it runs entirely before any job is submitted -- a budget failure must never cost shots.
"""
target = budget
for attempt in range(attempts):
circuits, bins = pack_ensemble(
num_qubits, occ, ranked, costs, target, n_circuits, parity=False
)
isa = pm.run(circuits)
worst = max(two_qubit_depth(c) for c in isa)
if worst <= budget:
return circuits, bins, isa
target = max(min(costs), int(target * budget / worst * 0.95))
raise RuntimeError(
f"could not fit {n_circuits} circuits inside a two-qubit depth of {budget} in "
f"{attempts} attempts; raise N_CIRCUITS or DEPTH_BUDGET and re-run this cell. "
"No QPU time was spent."
)

Pemprosesan pasca: baiki, gabungkan semula, pepenjurukan​

Tiga pembantu melakukan kerja Langkah 4.

half_configurations membahagikan setiap baris yang disampel kepada separuh proton dan separuh neutron, dan mengekalkan setiap separuh yang mempunyai bilangan nukleon yang betul. Satu baris dengan separuh proton yang sah menyumbang separuh itu walaupun separuh neutronnya mempunyai bilangan nukleon yang salah. Setiap separuh membawa jumlah berat sampel bagi baris yang ia muncul di dalamnya, yang menentukan kedudukannya jika subruang perlu dipangkas.

grow_subspace menggabungkan semula separuh-separuh tersebut kepada setiap hasil darab yang jatuh dalam sektor MJM_J dan pariti sasaran, menambah kepada subruang yang diberikan dan bukannya membina semula. Itu mengekalkan subruang berturutan bersarang, yang menjadikan jujukan tenaga menurun-monoton dan bukan sekadar berubah-ubah di sekitar sempadan.

recovery_loop ialah pemulihan konfigurasi swakonsisten bagi kertas SQD terkumpul [1]: baiki dua bilangan nukleon separuh-daftar terhadap anggaran kependudukan semasa, gabungkan semula, penjurukan, dan ambil anggaran kependudukan seterusnya daripada eigenvektor.

Semak konvensyen susunan-bit dengan teliti untuk mengelakkan keputusan yang salah. qiskit-addon-sqd menulis lajur 0 matriks rentetan-bitnya sebagai indeks qubit tertinggi, jadi menyongsangkan satu baris memberikan kependudukan mengikut indeks qubit; separuh "kanan"nya ialah indeks qubit rendah, iaitu blok proton. Selaras dengannya, recover_configurations mengambil num_elec_a sebagai bilangan proton dan kependudukan purata tersusun (protons, neutrons) mengikut indeks qubit. Tambahan tersebut menganggap bit ii berpasangan dengan bit i+Ni + N; dalam daftar ini, qubit proton ii dan qubit neutron i+Ni + N ialah keadaan (n,ℓ,j,mj)(n, \ell, j, m_j) yang sama, jadi andaian itu bermakna secara fizikal di sini dan bukan secara kebetulan.

def half_configurations(
bitstring_matrix, probabilities, sp, n_protons, n_neutrons
):
"""Split each row into proton and neutron halves, keeping each half on its own weight.

Column 0 of the addon's matrix is the highest qubit index, so reversing a row gives
occupation indexed by qubit.
"""
protons, neutrons = {}, {}
for row, weight in zip(
bitstring_matrix, np.asarray(probabilities, dtype=float)
):
occupied = np.flatnonzero(row[::-1])
p = tuple(int(i) for i in occupied if sp[i].tz == -1)
n = tuple(int(i) for i in occupied if sp[i].tz == +1)
if len(p) == n_protons:
protons[p] = protons.get(p, 0.0) + weight
if len(n) == n_neutrons:
neutrons[n] = neutrons.get(n, 0.0) + weight
return protons, neutrons

def product_subspace(sp, protons, neutrons, n_protons, n_neutrons, **targets):
"""Every (proton half) x (neutron half) product that lands in the target sector."""
return sorted(
d
for d in (
tuple(sorted(tuple(p) + tuple(n)))
for p in protons
for n in neutrons
)
if symmetry_allowed(sp, d, n_protons, n_neutrons, **targets)
)

def grow_subspace(
sp,
kept_protons,
kept_neutrons,
offered_protons,
offered_neutrons,
n_protons,
n_neutrons,
max_dimension=None,
**targets,
):
"""Add as many offered halves as the dimension cap allows, never dropping a kept one."""
kept_p, kept_n = list(kept_protons), list(kept_neutrons)
new_p = [c for c in offered_protons if c not in set(kept_p)]
new_n = [c for c in offered_neutrons if c not in set(kept_n)]

if max_dimension is None:
kept_p, kept_n = kept_p + new_p, kept_n + new_n
return (
product_subspace(
sp, kept_p, kept_n, n_protons, n_neutrons, **targets
),
kept_p,
kept_n,
)

basis = product_subspace(
sp, kept_p, kept_n, n_protons, n_neutrons, **targets
)
step = max(1, (len(new_p) + len(new_n)) // 24)
taken_p = taken_n = 0
while taken_p < len(new_p) or taken_n < len(new_n):
try_p, try_n = (
min(taken_p + step, len(new_p)),
min(taken_n + step, len(new_n)),
)
candidate = product_subspace(
sp,
kept_p + new_p[:try_p],
kept_n + new_n[:try_n],
n_protons,
n_neutrons,
**targets,
)
if len(candidate) > max_dimension:
if step == 1:
break
step = max(1, step // 2)
continue
basis, taken_p, taken_n = candidate, try_p, try_n
return basis, kept_p + new_p[:taken_p], kept_n + new_n[:taken_n]

def occupancies(sp, dets, vector):
"""Average occupancy of each qubit in a subspace eigenvector, as (protons, neutrons)."""
half = len(sp) // 2
occ = np.zeros(len(sp))
for weight, det in zip(np.abs(vector) ** 2, dets):
for q in det:
occ[q] += weight
return occ[:half], occ[half:]

def sample_occupancies(sp, bitstring_matrix, probabilities):
"""The same quantity estimated directly from sampled bitstrings."""
half = len(sp) // 2
weights = np.asarray(probabilities, dtype=float)
occ = (weights[:, None] * bitstring_matrix[:, ::-1]).sum(
axis=0
) / weights.sum()
return occ[:half], occ[half:]
def recovery_loop(
inter,
sp,
bitstring_matrix,
probabilities,
reference,
n_protons,
n_neutrons,
max_iterations=4,
energy_tol=1e-4,
max_dimension=None,
seed=None,
**targets,
):
"""Self-consistent configuration recovery, diagonalizing in the product subspace.

`num_elec_a` is the proton number and `num_elec_b` the neutron number, matching the
addon's right/left bipartition of the bitstring matrix.
"""
half = len(sp) // 2
p_ref = tuple(i for i in reference if sp[i].tz == -1)
n_ref = tuple(i for i in reference if sp[i].tz == +1)

survivors, survivor_probs = postselect_by_hamming_right_and_left(
bitstring_matrix,
np.asarray(probabilities, dtype=float).copy(),
hamming_right=n_protons,
hamming_left=n_neutrons,
)

if len(survivors):
guess = sample_occupancies(sp, survivors, survivor_probs)
else: # nothing survived: start from the reference itself
guess = (
np.array([1.0 if q in p_ref else 0.0 for q in range(half)]),
np.array(
[1.0 if q + half in n_ref else 0.0 for q in range(half)]
),
)

weights_p, weights_n = {p_ref: np.inf}, {n_ref: np.inf}
kept_p, kept_n = [p_ref], [n_ref]
history, best = [], None

for iteration in range(max_iterations):
# keep the occupancy estimate strictly inside (0, 1): the heuristic divides by it
clipped = tuple(np.clip(a, 1e-4, 1.0 - 1e-4) for a in guess)
recovered, recovered_probs = recover_configurations(
bitstring_matrix,
probabilities,
clipped,
n_protons,
n_neutrons,
rand_seed=None if seed is None else seed + iteration,
)

new_p, new_n = half_configurations(
recovered, recovered_probs, sp, n_protons, n_neutrons
)
for config, weight in new_p.items():
weights_p[config] = weights_p.get(config, 0.0) + weight
for config, weight in new_n.items():
weights_n[config] = weights_n.get(config, 0.0) + weight

def order(w):
return sorted(w, key=lambda c: (-w[c], c))

basis, kept_p, kept_n = grow_subspace(
sp,
kept_p,
kept_n,
order(weights_p),
order(weights_n),
n_protons,
n_neutrons,
max_dimension=max_dimension,
**targets,
)
energy, vector = ground_state(inter, basis)

history.append(
dict(
iteration=iteration + 1,
energy=energy,
dimension=len(basis),
protons=len(kept_p),
neutrons=len(kept_n),
recovered=len(recovered),
survivors=len(survivors),
)
)
print(
f" iteration {iteration + 1}: {len(kept_p)} proton x {len(kept_n)} neutron "
f"halves -> dimension {len(basis)}, E = {energy:.6f} MeV"
)

if best is None or energy < best[0]:
best = (energy, basis, vector)
guess = occupancies(
sp, basis, vector
) # the self-consistent update
if (
len(history) > 1
and abs(history[-2]["energy"] - energy) < energy_tol
):
break

return dict(
energy=best[0], basis=best[1], vector=best[2], history=history
)

Backend, bajet, dan parameter larian​

Setiap larian yang berikut menggunakan backend yang sama, pengurus laluan yang sama, dan bajet kedalaman yang sama, jadi ketiga-tiganya boleh dibandingkan secara langsung. Bajet menghubungkan mereka bersama: setiap litar dalam setiap ensembel mesti muat di dalamnya, dan ia menentukan berapa banyak kumpulan yang boleh disampel sama sekali.

Nilai di sini dipilih dengan mengukur kos yang ditranspilasi terhadap sasaran Heron. Pada kedalaman dua-qubit 300 dan 16 litar, kedua-dua ensembel 24-qubit dan 40-qubit muncul jauh di bawah 100 mikrosaat setiap litar, berbanding masa koheren beberapa ratus mikrosaat. Menambah bajet merangkumi lebih banyak kumpulan tetapi meningkatkan tempoh litar. Ukur tradeoff ini untuk backend anda.

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=40
)

pass_manager = generate_preset_pass_manager(
optimization_level=3, backend=backend, seed_transpiler=42
)
costing_manager = generate_preset_pass_manager(
optimization_level=1, backend=backend, seed_transpiler=42
)

DEPTH_BUDGET = 300 # two-qubit depth per circuit
N_CIRCUITS = 16 # circuits per ensemble
SHOTS = 10_000 # shots per circuit
MAX_DIMENSION = 4_000 # largest subspace the dense solver here will build
JOB_TAGS = ["TUT_SBQDNH"] # initials of the title's content words

# derive the two-qubit basis gate from the target by arity, not from a hard-coded name
two_qubit_basis = sorted(
name
for name in backend.target.operation_names
if backend.target.operation_from_name(name).num_qubits == 2
)
if not two_qubit_basis:
raise RuntimeError(
f"{backend.name} exposes no two-qubit gate; pick another backend"
)

print(
f"{backend.name}: {backend.num_qubits} qubits, two-qubit basis gate {two_qubit_basis[0]}"
)
print(
f"two-qubit depth budget {DEPTH_BUDGET}, {N_CIRCUITS} circuits x {SHOTS:,} shots per run"
)
print(
f"three runs: {3 * N_CIRCUITS} circuits, {3 * N_CIRCUITS * SHOTS:,} shots in total"
)
ibm_phoenix: 120 qubits, two-qubit basis gate cz
two-qubit depth budget 300, 16 circuits x 10,000 shots per run
three runs: 48 circuits, 480,000 shots in total

Contoh perkakasan berskala kecil​

Bahagian ini mengikuti aliran kerja empat-langkah pada QPU, menggunakan backend yang sama dan bajet get yang sama seperti larian skala-besar. Masalah yang lebih kecil menyediakan rujukan tepat untuk menyemak keputusan.

Masalah skala-kecil ialah 20Ne^{20}\mathrm{Ne}: dua proton valens dan dua neutron valens dalam cengkerang sdsd di atas teras 16O^{16}\mathrm{O}, dengan interaksi USDA [2]. Tiga orbital setiap spesies memberikan 24 qubit, dan asas yang dibenarkan-simetri yang lengkap ialah 640 penentu, cukup kecil untuk membandingkan anggaran tenaga dengan jawapan tepat.

Langkah 1: Memetakan input klasik kepada masalah kuantum​

Baca interaksi, bina daftar, dan bina penentu rujukan. Jadual berikut menunjukkan maklumat daftar daripada Background, dibaca terus daripada fail interaksi.

N_PROTONS, N_NEUTRONS = 2, 2

ms_sd = read_snt(DATA / "usda.snt", N_PROTONS, N_NEUTRONS)
sp_sd = m_scheme_states(ms_sd)
inter_sd = Interaction(ms_sd, sp_sd)
occ_sd = reference_determinant(sp_sd, inter_sd, N_PROTONS, N_NEUTRONS)

# post-selection splits the register in half, so the two species must contribute equally
n_proton_states = sum(1 for s in sp_sd if s.tz == -1)
if n_proton_states != len(sp_sd) - n_proton_states:
raise ValueError(
"this workflow needs equal proton and neutron state counts"
)

SHELL_LABEL = {0: "s", 1: "p", 2: "d", 3: "f", 4: "g"}
print(
f"core Z={ms_sd.core_z} N={ms_sd.core_n} plus {N_PROTONS}p + {N_NEUTRONS}n valence "
f"-> A={ms_sd.mass_number} on {len(sp_sd)} qubits"
)
print(
f"interaction: {len(ms_sd.tbmes)} J-coupled matrix elements, fitted at "
f"A_ref={ms_sd.a_ref}, rescaled by (A/A_ref)^{ms_sd.mass_exponent:g} = "
f"{ms_sd.mass_factor:.6f}\n"
)

print(
f"{'orbital':>9} {'SPE (MeV)':>10} {'proton qubits':>14} {'neutron qubits':>15}"
)
for o in (o for o in ms_sd.orbitals if o.tz == -1):
twin = next(
t
for t in ms_sd.orbitals
if t.tz == +1 and (t.n, t.ell, t.j2) == (o.n, o.ell, o.j2)
)
qp = [q for q, s in enumerate(sp_sd) if s.orb == o.idx]
qn = [q for q, s in enumerate(sp_sd) if s.orb == twin.idx]
print(
f"{f'{o.n}{SHELL_LABEL[o.ell]}{o.j2}/2':>9} {ms_sd.spes[o.idx]:>10.4f} "
f"{f'{qp[0]}-{qp[-1]}':>14} {f'{qn[0]}-{qn[-1]}':>15}"
)

print(f"\nreference determinant occupies qubits {occ_sd}")
print(
f" M_J = {sum(sp_sd[i].mj2 for i in occ_sd) / 2:g}, "
f"parity = {(-1) ** (sum(sp_sd[i].ell for i in occ_sd) % 2):+d}, "
f"energy = {matrix_element(inter_sd, occ_sd, occ_sd):.6f} MeV"
)
core Z=8 N=8 plus 2p + 2n valence -> A=20 on 24 qubits
interaction: 158 J-coupled matrix elements, fitted at A_ref=18, rescaled by (A/A_ref)^-0.3 = 0.968886

orbital SPE (MeV) proton qubits neutron qubits
0d3/2 2.1117 0-3 12-15
0d5/2 -3.9257 4-9 16-21
1s1/2 -3.2079 10-11 22-23

reference determinant occupies qubits (4, 9, 16, 21)
M_J = 0, parity = +1, energy = -29.765549 MeV

Jalankan dua semakan ke atas Hamiltonian sebelum meneruskan. Kedua-duanya murah dan boleh mendedahkan ralat penggandingan-semula yang mungkin tidak dapat dikesan oleh satu pengiraan tenaga sahaja.

Hamiltonian yang tak varian secara putaran menyusun eigenkeadaannya ke dalam multiplet JJ, jadi setiap eigennilai sektor MJ=2M_J = 2 mesti juga muncul dalam spektrum MJ=0M_J = 0 pada tenaga yang sama. Jurang antara keadaan-asas dan keadaan terendah yang membawa MJ=2M_J = 2 ialah tenaga pengujaan 2+2^+, yang diukur: 1.6341.634 MeV untuk 20Ne^{20}\mathrm{Ne} [6]. Interaksi cengkerang-sdsd empirik dijangka bersetuju dalam lingkungan beberapa ratus keV.

basis_exact_sd = full_basis(sp_sd, N_PROTONS, N_NEUTRONS)
if len(basis_exact_sd) != count_basis(sp_sd, N_PROTONS, N_NEUTRONS):
raise AssertionError("the basis counter disagrees with the enumeration")
E_REF_SD = matrix_element(inter_sd, occ_sd, occ_sd)
E_EXACT_SD, _ = ground_state(inter_sd, basis_exact_sd)

# the M_J = 2 sector: its spectrum must be contained in the M_J = 0 spectrum
basis_mj2 = full_basis(sp_sd, N_PROTONS, N_NEUTRONS, mj2_target=4)
spectrum_0 = np.linalg.eigvalsh(
subspace_hamiltonian(inter_sd, basis_exact_sd)
)
spectrum_2 = np.linalg.eigvalsh(subspace_hamiltonian(inter_sd, basis_mj2))
contained = sum(
1 for e in spectrum_2 if np.min(np.abs(spectrum_0 - e)) < 1e-7
)
if contained != len(spectrum_2):
raise AssertionError(
f"rotational invariance broken: only {contained}/{len(spectrum_2)} "
"M_J=2 eigenvalues appear in the M_J=0 spectrum"
)

print(
f"rotational invariance: all {contained} M_J=2 eigenvalues found in the M_J=0 spectrum"
)
print(
f"E(2+) - E(0+) = {spectrum_2[0] - E_EXACT_SD:.3f} MeV (experiment: 1.634 MeV)\n"
)
print(f"reference determinant {E_REF_SD:11.6f} MeV")
print(
f"exact diagonalization {E_EXACT_SD:11.6f} MeV (dimension {len(basis_exact_sd)})"
)
print(f"correlation energy to find {E_EXACT_SD - E_REF_SD:11.6f} MeV")
rotational invariance: all 497 M_J=2 eigenvalues found in the M_J=0 spectrum
E(2+) - E(0+) = 1.747 MeV (experiment: 1.634 MeV)

reference determinant -29.765549 MeV
exact diagonalization -40.472331 MeV (dimension 640)
correlation energy to find -10.706782 MeV

Seterusnya, bina kumpulan operator. Mengaplikasikan dua peraturan pemilihan memberikan hasil yang penting: bagi rujukan ini, dalam ruang model ini, tiada langsung pengujaan tunggal yang dibenarkan.

Sebabnya khusus dan boleh disemak. Pengujaan 1p1h1p1h mengekalkan MJM_J hanya jika keadaan zarah mempunyai mjm_j yang sama dengan lubang tersebut. Rujukan menduduki dua keadaan ∣mj∣|m_j| terbesar dalam orbital terendah (mj=±5/2m_j = \pm 5/2 bagi 0d5/20d_{5/2}), dan tiada orbital lain dalam cengkerang sdsd mencapai ∣mj∣=5/2|m_j| = 5/2, kerana 0d3/20d_{3/2} berhenti pada 3/23/2 dan 1s1/21s_{1/2} pada 1/21/2. Oleh itu, tiada pengujaan tunggal bertahan, dan korelasi dibawa sepenuhnya oleh pengujaan 2p2h2p2h. Ini ialah sifat rujukan dan cengkerang tersebut, bukan hukum umum; sel berikut mengiranya dan bukan menganggapnya.

raw_pool_sd = excitation_pool(sp_sd, occ_sd)
pool_sd = [op for op in raw_pool_sd if conserves_symmetry(sp_sd, op)]
ranked_sd = rank_pool(inter_sd, occ_sd, pool_sd)

singles_sd = [
(h, v)
for h in occ_sd
for v in range(len(sp_sd))
if v not in occ_sd and sp_sd[h].tz == sp_sd[v].tz
]
singles_mj_sd = [
(h, v) for h, v in singles_sd if sp_sd[h].mj2 == sp_sd[v].mj2
]

print(
f"1p1h: {len(singles_sd):4d} raw -> {len(singles_mj_sd):3d} conserve M_J"
)
print(
f"2p2h: {len(raw_pool_sd):4d} raw -> {len(pool_sd):3d} conserve M_J and couple to a common J\n"
)

print(
f"{'rank':>4} {'holes':>9} {'particles':>11} {'<ref|H|a> (MeV)':>16} {'amplitude':>10}"
)
for r, (op, coupling, amplitude) in enumerate(ranked_sd[:8], start=1):
print(
f"{r:>4} {f'{op[0]},{op[1]}':>9} {f'{op[2]},{op[3]}':>11} "
f"{coupling:>16.4f} {amplitude:>10.4f}"
)

# what is the best this ansatz could possibly do? Apply every excitation once and recombine.
reachable = {occ_sd}
for op, _, _ in ranked_sd:
h1, h2, v1, v2 = op
reachable |= {
tuple(sorted(set(d) - {h1, h2} | {v1, v2}))
for d in reachable
if {h1, h2} <= set(d) and not {v1, v2} & set(d)
}
ceiling = product_subspace(
sp_sd,
{tuple(i for i in d if sp_sd[i].tz == -1) for d in reachable},
{tuple(i for i in d if sp_sd[i].tz == +1) for d in reachable},
N_PROTONS,
N_NEUTRONS,
)
print(
f"\nthe pool reaches {len(reachable)} determinants, whose product subspace spans "
f"{len(ceiling)} of {len(basis_exact_sd)}"
)
1p1h: 40 raw -> 0 conserve M_J
2p2h: 490 raw -> 78 conserve M_J and couple to a common J

rank holes particles <ref|H|a> (MeV) amplitude
1 4,9 0,3 -1.8714 0.1375
2 16,21 12,15 -1.8714 0.1375
3 4,21 3,12 1.6775 -0.1029
4 9,16 0,15 1.6775 -0.1029
5 4,9 10,11 -0.8728 0.1168
6 16,21 22,23 -0.8728 0.1168
7 4,21 3,17 1.0622 -0.0882
8 4,21 8,12 -1.0622 0.0882

the pool reaches 412 determinants, whose product subspace spans 640 of 640

Langkah 2: Optimumkan masalah untuk pelaksanaan perkakasan kuantum​

Transpilasi mendedahkan kos perkakasan bagi rentetan ZZ Jordan-Wigner dan penjimatan daripada menggunakan pengujaan qubit. Sel pertama mengukur kedua-dua pembinaan terhadap sasaran backend sebenar dan menyemak dakwaan, yang diperkenalkan dalam Setup, bahawa menjatuhkan rentetan ZZ mengubah amplitud tetapi bukan set penentu yang boleh dicapai litar.

Bandingkan dua akibat penggantian ini. Pengujaan qubit mengenakan kos yang sama tanpa mengira jarak antara indeksnya, jadi pengujaan proton-neutron, yang merentangi sempadan antara dua separuh daftar dan membentuk kebanyakan kumpulan, tidak lagi mempunyai kos tambahan ini. Keseluruhan kumpulan kemudian muat di dalam bajet, yang bermaksud had ke atas keputusan ialah pensampelan dan bukan kedalaman litar.

# 1. do the two constructions reach the same determinants?
# Apply one block to the reference on the window it spans and read off which basis states
# acquire amplitude. Column 0 of the unitary is the image of |0...0>, and the X gates that
# place the reference are part of the circuit, so that column is exactly what is wanted.
# A fermionic block's window is its whole span, and building a unitary on it costs 4^n, so
# probe the narrowest excitations in the pool rather than the highest-ranked ones.
PROBE_SPAN = 12
narrow = sorted(ranked_sd, key=lambda row: max(row[0]) - min(row[0]))
probes = [op for op, _, _ in narrow if max(op) - min(op) + 1 <= PROBE_SPAN][
:3
]
if len(probes) < 2:
raise RuntimeError(
f"no excitation spans {PROBE_SPAN} qubits or fewer; raise PROBE_SPAN"
)

print(
f"{'excitation':>16} {'span':>5} {'reachable determinants':>22} {'same as fermionic?':>19}"
)
for probe_op in probes:
probe_window = list(range(min(probe_op), max(probe_op) + 1))
probe_local = tuple(probe_window.index(i) for i in probe_op)
probe_occ = tuple(
probe_window.index(i) for i in occ_sd if i in probe_window
)

supports = {}
for parity in (True, False):
unitary = Operator(
excitation_ansatz(
len(probe_window),
probe_occ,
[probe_local],
[0.7],
measure=False,
parity=parity,
)
).data
supports[parity] = frozenset(
np.flatnonzero(np.abs(unitary[:, 0]) > 1e-10).tolist()
)

if len(supports[True]) < 2:
raise AssertionError(
f"{probe_op}: the block did not move any amplitude, so this "
"comparison would be vacuous"
)
if supports[True] != supports[False]:
raise AssertionError(
f"{probe_op}: the two constructions reach different determinants"
)
print(
f"{str(probe_op):>16} {len(probe_window):>5} {len(supports[True]):>22} {'yes':>19}"
)

print(
"\n-> identical support; the amplitudes differ, and pooled SQD only consumes the support\n"
)

# 2. what does each one cost on this backend?
cost_qeb = excitation_costs(
len(sp_sd), ranked_sd, costing_manager, parity=False
)
cost_jw = excitation_costs(
len(sp_sd), ranked_sd, costing_manager, parity=True
)

def species(op):
return "same" if len({sp_sd[i].tz for i in op}) == 1 else "pn"

print(
f"{'excitation':>10} {'count':>5} {'QEB 2q depth':>14} {'fermionic 2q depth':>19}"
)
for group in ("same", "pn"):
q = [
c
for (op, _, _), c in zip(ranked_sd, cost_qeb)
if species(op) == group
]
j = [
c for (op, _, _), c in zip(ranked_sd, cost_jw) if species(op) == group
]
print(
f"{group:>10} {len(q):>5} {f'{min(q)}-{max(q)}':>14} {f'{min(j)}-{max(j)}':>19}"
)
print(
f"{'pool total':>10} {len(ranked_sd):>5} {sum(cost_qeb):>14} {sum(cost_jw):>19}"
)
print(
f"\nfermionic / qubit-excitation cost ratio: {sum(cost_jw) / sum(cost_qeb):.2f}x"
)
print(
f"\nensemble capacity: {N_CIRCUITS} circuits at two-qubit depth {DEPTH_BUDGET}"
)
excitation span reachable determinants same as fermionic?
(4, 9, 5, 8) 6 2 yes
(16, 21, 17, 20) 6 2 yes
(4, 9, 6, 7) 6 2 yes

-> identical support; the amplitudes differ, and pooled SQD only consumes the support

excitation count QEB 2q depth fermionic 2q depth
same 26 40-48 48-144
pn 52 48-48 48-256
pool total 78 3728 9112

fermionic / qubit-excitation cost ratio: 2.44x

ensemble capacity: 16 circuits at two-qubit depth 300
circuits_sd, bins_sd, isa_sd = pack_to_budget(
len(sp_sd),
occ_sd,
ranked_sd,
cost_qeb,
DEPTH_BUDGET,
N_CIRCUITS,
pass_manager,
)

PACKED_SD = sum(len(b) for b in bins_sd)
worst_sd = max(two_qubit_depth(c) for c in isa_sd)
worst_count_sd = max(two_qubit_count(c) for c in isa_sd)

print(
f"packed {PACKED_SD} of {len(ranked_sd)} excitations into {N_CIRCUITS} circuits"
)
print(f" excitations per circuit {[len(b) for b in bins_sd]}")
print(f" two-qubit depth {[two_qubit_depth(c) for c in isa_sd]}")
print(f" two-qubit gates {[two_qubit_count(c) for c in isa_sd]}")
print(
f"\nworst circuit: two-qubit depth {worst_sd} of a {DEPTH_BUDGET} budget, "
f"{worst_count_sd} two-qubit gates"
)
packed 78 of 78 excitations into 16 circuits
excitations per circuit [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4]
two-qubit depth [228, 182, 177, 220, 214, 220, 226, 224, 222, 222, 181, 179, 136, 179, 171, 181]
two-qubit gates [234, 231, 227, 228, 225, 225, 233, 229, 233, 226, 230, 223, 232, 226, 176, 185]

worst circuit: two-qubit depth 228 of a 300 budget, 234 two-qubit gates

Langkah 3: Laksanakan menggunakan primitif Qiskit​

Hantar satu job setiap masalah, dengan keseluruhan ensembel sebagai satu senarai litar tunggal. Twirling get dan pengukuran serta penyahgandingan dinamik didayakan untuk mengurangkan kesan hingar perkakasan. Manfaatnya bergantung kepada litar dan backend.

ID setiap job dicetak. Gunakan service.job("JOB_ID") untuk mendapatkan semula job yang telah selesai dan keputusannya tanpa menggunakan masa QPU tambahan.

def sample(isa_circuits, shots, tags):
"""Submit one Sampler job; return the per-circuit bit arrays and the measured QPU seconds."""
sampler = SamplerV2(mode=backend)
sampler.options.environment.job_tags = tags
sampler.options.twirling.enable_gates = True
sampler.options.twirling.enable_measure = True
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"

job = sampler.run(isa_circuits, shots=shots)
print(
f"job {job.job_id()}: {len(isa_circuits)} circuits x {shots:,} shots "
f"on {backend.name}"
)
return [pub.data.meas for pub in job.result()]

def pool_samples(bit_arrays, sp):
"""Merge the ensemble's bit arrays into one bitstring matrix and probability vector."""
matrices, weights, total = [], [], 0
for bit_array in bit_arrays:
matrix, probabilities = bit_array_to_arrays(bit_array)
matrices.append(matrix)
weights.append(probabilities * bit_array.num_shots)
total += bit_array.num_shots
counts = np.concatenate(weights)
matrix = np.vstack(matrices)
# the same bitstring can appear in more than one circuit; merge duplicate rows
unique, inverse = np.unique(matrix, axis=0, return_inverse=True)
merged = np.zeros(len(unique))
np.add.at(merged, inverse.ravel(), counts)
return unique, merged / merged.sum(), total
bit_arrays_sd = sample(isa_sd, SHOTS, JOB_TAGS + ["20Ne"])
matrix_sd, probs_sd, shots_sd = pool_samples(bit_arrays_sd, sp_sd)

survivors_sd, _ = postselect_by_hamming_right_and_left(
matrix_sd,
probs_sd.copy(),
hamming_right=N_PROTONS,
hamming_left=N_NEUTRONS,
)
shot_survival_sd = float(
probs_sd[
(matrix_sd[:, len(sp_sd) // 2 :].sum(axis=1) == N_PROTONS)
& (matrix_sd[:, : len(sp_sd) // 2].sum(axis=1) == N_NEUTRONS)
].sum()
)

reference_bits = "".join(
"1" if q in occ_sd else "0" for q in range(len(sp_sd))
)[::-1]
print(f"\n{shots_sd:,} shots -> {len(matrix_sd):,} distinct bitstrings")
print(
f" {shot_survival_sd:6.1%} of shots carry the right proton and neutron numbers"
)
print(f" {len(survivors_sd):,} distinct bitstrings do")

order = np.argsort(-probs_sd)
half = len(sp_sd) // 2
print(f"\n{'neutrons':>{half}} | {'protons':<{half}} share")
for i in order[:4]:
bits = "".join("1" if b else "0" for b in matrix_sd[i])
tag = " <- reference determinant" if bits == reference_bits else ""
print(f"{bits[:half]} | {bits[half:]} {probs_sd[i]:6.2%}{tag}")

if len(survivors_sd) == 0:
raise RuntimeError(
"no shot carried the right nucleon numbers; check the backend and "
"the transpiled circuits before spending more QPU time"
)
job dap30qtr85ps73fg21p0: 16 circuits x 10,000 shots on ibm_phoenix

160,000 shots -> 17,221 distinct bitstrings
31.5% of shots carry the right proton and neutron numbers
973 distinct bitstrings do

neutrons | protons share
001000010000 | 001000010000 20.08% <- reference determinant
000000010000 | 001000010000 2.25%
001000010000 | 001000000000 2.19%
001000010000 | 000000010000 1.93%

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang dikehendaki​

Tukar sampel kuantum kepada anggaran tenaga menggunakan kekangan simetri-nuklear yang diterangkan dalam bahagian Latar Belakang.

Pemulihan konfigurasi membaiki dua bilangan nukleon. recover_configurations mengambil setiap shot yang mempunyai bilangan proton atau neutron yang salah dan menukar bit yang paling tidak konsisten dengan anggaran kependudukan orbital purata semasa, dan bukannya membuangnya. Pada laluan pertama anggaran kependudukan datang daripada shot yang telah terselamat; selepas itu ia datang daripada eigenvektor subruang sebelumnya, yang menjadikan prosedur ini konsisten-kendiri.

MJM_J dan pariti dikenakan pada produk yang digabung semula, bukan pada keseluruhan shot. Setiap shot yang dibaiki menyumbang separuh proton dan separuh neutron, dan subruang direntangi oleh setiap produk konfigurasi proton yang disampel dengan konfigurasi neutron yang disampel yang jatuh pada MJ=0M_J = 0 dengan pariti yang betul. Menapis keseluruhan shot berdasarkan MJM_J jumlah pula akan membuang dua separuh yang baik demi satu nombor kuantum yang sebenarnya kepunyaan gabungan mereka.

Empat semakan nombor kuantum menolak pecahan sampel yang berbeza. Dua bilangan nukleon mengambil kira sebahagian besar penapisan. Pariti secara automatik dipenuhi di dalam satu kelompok major: setiap orbital sdsd mempunyai ℓ\ell genap dan setiap orbital pfpf mempunyai ℓ\ell ganjil, jadi sebaik sahaja bilangan nukleon betul, pariti tidak boleh salah. Semakan pariti dikekalkan kerana ruang model silang-kelompok akan menjadikannya kekangan bebas. Semakan MJM_J mengekalkan produk dalam sektor momentum-sudut sasaran. Nilai mempunyai empat nombor kuantum tepat ialah kerana ia murah dan tepat, bukan kerana setiap satu ialah penapis besar.

Diagonalisasi memberikan sempadan atas variasional. Kerana subruang setiap iterasi mengandungi subruang sebelumnya, jujukan tenaga menurun secara monotonik, dan setiap entri di dalamnya ialah sempadan atas yang tegar bagi tenaga keadaan-dasar sebenar, tanpa mengira hingar dalam sampel yang menghasilkannya.

result_sd = recovery_loop(
inter_sd,
sp_sd,
matrix_sd,
probs_sd,
occ_sd,
N_PROTONS,
N_NEUTRONS,
max_iterations=4,
max_dimension=MAX_DIMENSION,
seed=42,
)

E_SQD_SD = result_sd["energy"]
recovered_sd = 100 * (E_SQD_SD - E_REF_SD) / (E_EXACT_SD - E_REF_SD)

print(f"\nreference determinant {E_REF_SD:11.6f} MeV")
print(
f"pooled SQD upper bound {E_SQD_SD:11.6f} MeV "
f"(subspace dimension {len(result_sd['basis'])} of {len(basis_exact_sd)})"
)
print(f"exact diagonalization {E_EXACT_SD:11.6f} MeV")
print(f"\ncorrelation energy recovered: {recovered_sd:.1f}%")

energies_sd = [h["energy"] for h in result_sd["history"]]
if any(b > a + 1e-9 for a, b in zip(energies_sd, energies_sd[1:])):
raise AssertionError(
"the subspaces are not nested; the bound should never rise"
)
if E_SQD_SD < E_EXACT_SD - 1e-7:
raise AssertionError(
f"pooled SQD returned {E_SQD_SD:.6f}, below the exact {E_EXACT_SD:.6f}; "
"a subspace bound cannot beat the full diagonalization"
)
iteration 1: 66 proton x 66 neutron halves -> dimension 640, E = -40.472331 MeV
iteration 2: 66 proton x 66 neutron halves -> dimension 640, E = -40.472331 MeV

reference determinant -29.765549 MeV
pooled SQD upper bound -40.472331 MeV (subspace dimension 640 of 640)
exact diagonalization -40.472331 MeV

correlation energy recovered: 100.0%

Nilaikan hasil​

Gunakan semakan berikut untuk menilai hasil anda pada backend kelas Heron dengan tetapan ini:

  • Kadar kelangsungan shot pada dua bilangan nukleon mengukur pecahan shot dengan bilangan proton dan neutron yang betul. Ia boleh menurun apabila register membesar. Kadar kelangsungan yang hampir sifar boleh menandakan masalah dengan pelaksanaan litar. Semak kedalaman ISA dalam Langkah 2 dan penentukuran backend, bukan pasca-pemprosesan.

  • Gelung pemulihan patut mencetak dimensi subruang yang kekal malar atau membesar dan tenaga yang kekal malar atau menurun dengan setiap iterasi. Jika iterasi 1 sudah mencapai MAX_DIMENSION, penyelesai klasik dan bukannya pensampelan adalah kekangan pengikat.

  • Pecahan yang dipulihkan untuk 20Ne^{20}\mathrm{Ne} patut tinggi, kerana siling ansatz yang dikira dalam Langkah 1 ialah ruang penentu 640 penuh; larian ini ialah tempat pensampelan, bukan keekspresifan, adalah satu-satunya halangan.

  • Dua pengesahan dalam sel sebelumnya menyemak sempadan variasional. Sempadan yang meningkat bermaksud subruang berhenti bersarang, dan sempadan di bawah tenaga tepat bermaksud ada sesuatu yang salah dengan Hamiltonian, bukan dengan perkakasan.

Secara berlawanan dengan intuisi, backend yang lebih berhingar boleh memberikan sempadan yang sedikit lebih baik daripada backend yang bersih, kerana ralat menghasilkan separuh-konfigurasi sah yang tidak akan pernah disampel oleh litar ideal, dan meluaskan subruang variasional tidak boleh menaikkan nilai eigen terendahnya. Simulasi berhingar boleh menunjukkan kesan yang sama; tutorial ini menunjukkannya dengan sampel perkakasan.

# IBM Carbon palette: Blue 60 and Blue 80 for data, Gray 100/70/30 for ink and rules
SURFACE, INK, MUTED, RULE = "#ffffff", "#161616", "#6f6f6f", "#c6c6c6"
SERIES, DEEP, PURPLE = "#0f62fe", "#002d9c", "#6929c4"

def convergence_plot(
history, e_ref, e_exact, title, colour=SERIES, full_dim=None
):
"""Energy against subspace dimension, scaled to the data rather than to the full window.

A good run lands within a fraction of a percent of the exact answer, so an axis spanning
reference-to-exact would squash every point onto one line. The axis is therefore scaled to
the data (plus the exact line, when there is one), and the right-hand axis carries the
fraction of the correlation energy so the absolute and relative readings sit side by side.
"""
dimensions = [h["dimension"] for h in history]
energies = [h["energy"] for h in history]

fig, ax = plt.subplots(figsize=(7.4, 4.3), facecolor=SURFACE)
ax.set_facecolor(SURFACE)
ax.plot(
dimensions,
energies,
"-o",
color=colour,
linewidth=2,
markersize=8,
markeredgecolor=SURFACE,
markeredgewidth=1.5,
zorder=3,
)
stacked = {}
for h in history:
# a converged loop repeats the same point; stack the labels so they do not overprint
key = (round(h["dimension"]), round(h["energy"], 9))
offset = 12 + 11 * stacked.get(key, 0)
stacked[key] = stacked.get(key, 0) + 1
ax.annotate(
str(h["iteration"]),
xy=(h["dimension"], h["energy"]),
xytext=(0, offset),
textcoords="offset points",
ha="center",
fontsize=8,
color=MUTED,
)

span = (max(dimensions) - min(dimensions)) or max(1, max(dimensions) // 4)
x_left, x_right = (
min(dimensions) - 0.14 * span,
max(dimensions) + 0.40 * span,
)
ax.set_xlim(x_left, x_right)

floor = min(energies) if e_exact is None else min(min(energies), e_exact)
height = max(max(energies) - floor, 1e-3)
ax.set_ylim(floor - 0.30 * height, max(energies) + 0.42 * height)

if e_exact is not None:
ax.axhline(
e_exact, color=MUTED, linestyle="--", linewidth=1, zorder=1
)
label = "exact" + (f", {full_dim:,} determinants" if full_dim else "")
ax.annotate(
f"{label} {e_exact:.3f} MeV".replace("-", "\u2212"),
xy=(x_left, e_exact),
xytext=(3, 5),
textcoords="offset points",
ha="left",
va="bottom",
color=MUTED,
fontsize=9,
)

# the reference determinant is far off this scale; state it rather than plotting it
ax.annotate(
f"reference determinant {e_ref:.3f} MeV".replace("-", "\u2212")
+ f" ({e_ref - max(energies):+.2f} MeV off the top of this axis)".replace(
"-", "\u2212"
),
xy=(x_right, max(energies) + 0.42 * height),
xytext=(-3, -12),
textcoords="offset points",
ha="right",
va="top",
color=MUTED,
fontsize=8.5,
)

if e_exact is not None and abs(e_exact - e_ref) > 1e-9:
right = ax.twinx()
low, high = ax.get_ylim()

def to_percent(e):
return 100 * (e - e_ref) / (e_exact - e_ref)

right.set_ylim(to_percent(low), to_percent(high))
right.set_ylabel("correlation energy recovered (%)", color=MUTED)
right.tick_params(colors=MUTED)
for side in ("top", "left"):
right.spines[side].set_visible(False)
right.spines["right"].set_color(MUTED)
right.spines["bottom"].set_color(MUTED)

ax.set_xlabel("subspace dimension", color=MUTED)
ax.set_ylabel("ground-state energy (MeV)", color=MUTED)
ax.set_title(title, color=INK, fontsize=11.5, loc="left", pad=12)
ax.grid(axis="y", color=RULE, alpha=0.55)
ax.tick_params(colors=MUTED)
for side in ("top", "right"):
ax.spines[side].set_visible(False)
for side in ("bottom", "left"):
ax.spines[side].set_color(MUTED)
fig.tight_layout()
return fig

convergence_plot(
result_sd["history"],
E_REF_SD,
E_EXACT_SD,
f"$^{{20}}$Ne: the bound falls as configuration recovery widens the subspace\n"
f"({backend.name}, {N_CIRCUITS} circuits x {SHOTS:,} shots)",
full_dim=len(basis_exact_sd),
)
plt.show()

Output of the previous code cell

Contoh perkakasan berskala besar​

Menskala naik hanya menukar input, jadi langkah seterusnya ialah menggabungkan empat peringkat menjadi satu fungsi dan menjalankannya dua kali, kedua-duanya pada register 40-qubit dalam kelompok pfpf di atas teras 40Ca^{40}\mathrm{Ca} dengan interaksi GXPF1 [3].

Dua larian ini menggambarkan aspek penskalaan yang berbeza:

  • 44Ti^{44}\mathrm{Ti}, dua proton valens dan dua neutron valens, mempunyai asas 4,000-penentu. Register ialah 40 qubit, tetapi masalah masih cukup kecil untuk didiagonalisasi secara tepat pada komputer riba, jadi anda boleh membandingkan hasil perkakasan dengan rujukan tepat selepas menaikkan saiz register.

  • 48Cr^{48}\mathrm{Cr}, empat proton valens dan empat neutron valens, mempunyai 1,963,461 penentu yang dibenarkan-simetri dalam 40 qubit yang sama. Penyelesai tumpat tutorial ini tidak boleh mendiagonalisasi ruang penuh itu, jadi larian ini mengembalikan sempadan atas yang tegar dan penentu rujukan yang diperbaikinya.

Perhatikan dua kuantiti merentas kedua-dua larian. Pecahan kolam yang muat di dalam belanjawan get tetap mengecil apabila kolam membesar, dan pack_ensemble melaporkan berapa banyak yang disertakan. Subruang berhenti dihadkan oleh pensampelan dan mula dihadkan oleh MAX_DIMENSION, matriks terbesar yang dibina oleh penyelesai klasik tumpat di sini. Pada skala ini, pengiraan pengeluaran akan menggunakan penyelesai interaksi konfigurasi terpilih (selected-CI).

Gabungkan langkah 1–4​

Fungsi berikut memanggil peringkat yang sama seperti dalam panduan langkah demi langkah, dalam susunan yang sama.

def sqd_run(snt_file, n_protons, n_neutrons, name, exact=True):
"""The whole workflow for one nucleus. Returns a record of every stage."""
# -------------------------Step 1-------------------------
ms = read_snt(DATA / snt_file, n_protons, n_neutrons)
sp = m_scheme_states(ms)
inter = Interaction(ms, sp)
if sum(1 for s in sp if s.tz == -1) * 2 != len(sp):
raise ValueError(
f"{name}: post-selection needs equal proton and neutron state counts"
)
reference = reference_determinant(sp, inter, n_protons, n_neutrons)
e_ref = matrix_element(inter, reference, reference)

raw = excitation_pool(sp, reference)
ranked = rank_pool(
inter, reference, [op for op in raw if conserves_symmetry(sp, op)]
)
print(
f"{name}: {len(sp)} qubits, {n_protons}p + {n_neutrons}n, A = {ms.mass_number}"
)
print(
f" 2p2h pool {len(raw)} raw -> {len(ranked)} symmetry-allowed; "
f"reference energy {e_ref:.6f} MeV"
)

# -------------------------Step 2-------------------------
costs = excitation_costs(len(sp), ranked, costing_manager)
circuits, bins, isa = pack_to_budget(
len(sp),
reference,
ranked,
costs,
DEPTH_BUDGET,
N_CIRCUITS,
pass_manager,
)
packed = sum(len(b) for b in bins)
worst = max(two_qubit_depth(c) for c in isa)
worst_count = max(two_qubit_count(c) for c in isa)
print(
f" packed {packed} of {len(ranked)} excitations; worst circuit two-qubit depth "
f"{worst}, {worst_count} two-qubit gates"
)

# -------------------------Step 3-------------------------
# a unique tag per run, so the jobs are findable later
bit_arrays = sample(isa, SHOTS, JOB_TAGS + [name])
matrix, probabilities, shots = pool_samples(bit_arrays, sp)
survival = float(
probabilities[
(matrix[:, len(sp) // 2 :].sum(axis=1) == n_protons)
& (matrix[:, : len(sp) // 2].sum(axis=1) == n_neutrons)
].sum()
)
print(
f" {shots:,} shots -> {len(matrix):,} distinct bitstrings, "
f"{survival:.1%} of shots with the right nucleon numbers"
)
if survival == 0.0:
raise RuntimeError(
f"{name}: no shot carried the right nucleon numbers"
)

# -------------------------Step 4-------------------------
result = recovery_loop(
inter,
sp,
matrix,
probabilities,
reference,
n_protons,
n_neutrons,
max_iterations=4,
max_dimension=MAX_DIMENSION,
seed=42,
)
energy = result["energy"]

full_dim = count_basis(sp, n_protons, n_neutrons) # cheap, even when huge
e_exact = None
if exact:
full = full_basis(sp, n_protons, n_neutrons)
if len(full) != full_dim:
raise AssertionError(
f"{name}: counted {full_dim} determinants but enumerated "
f"{len(full)}"
)
e_exact, _ = ground_state(inter, full)

print(f" reference {e_ref:11.6f} MeV pooled SQD {energy:11.6f} MeV")
if e_exact is not None:
print(
f" exact {e_exact:11.6f} MeV (dimension {full_dim}) -> "
f"{100 * (energy - e_ref) / (e_exact - e_ref):.1f}% of the correlation energy"
)
if energy < e_exact - 1e-7:
raise AssertionError(
f"{name}: pooled SQD bound is below the exact energy"
)
else:
print(
f" no exact reference: the symmetry-allowed basis is {full_dim:,} determinants"
)
print(
f" the bound captures {energy - e_ref:.6f} MeV of correlation energy"
)
print()

return dict(
name=name,
qubits=len(sp),
pool=len(ranked),
packed=packed,
two_qubit=worst,
two_qubit_gates=worst_count,
shots=shots,
distinct=len(matrix),
survival=survival,
dimension=len(result["basis"]),
full_dim=full_dim,
e_ref=e_ref,
e_sqd=energy,
e_exact=e_exact,
history=result["history"],
# the subspace and its eigenvector cannot be reconstructed from the summary --
# they depend on the sampled shots -- so keep them for the scaling analysis
interaction=inter,
states=sp,
reference=reference,
ranked=ranked,
basis=result["basis"],
vector=result["vector"],
)

pretty = {"20Ne": "$^{20}$Ne", "44Ti": "$^{44}$Ti", "48Cr": "$^{48}$Cr"}

small_scale = dict(
name="20Ne",
qubits=len(sp_sd),
pool=len(ranked_sd),
packed=PACKED_SD,
two_qubit=worst_sd,
two_qubit_gates=worst_count_sd,
shots=shots_sd,
distinct=len(matrix_sd),
survival=shot_survival_sd,
dimension=len(result_sd["basis"]),
full_dim=len(basis_exact_sd),
e_ref=E_REF_SD,
e_sqd=E_SQD_SD,
e_exact=E_EXACT_SD,
history=result_sd["history"],
interaction=inter_sd,
states=sp_sd,
reference=occ_sd,
ranked=ranked_sd,
basis=result_sd["basis"],
vector=result_sd["vector"],
)

44Ti^{44}\mathrm{Ti}: aliran kerja yang sama pada register 40-qubit​

Kelompok pfpf di atas 40Ca^{40}\mathrm{Ca} mempunyai empat orbital bagi setiap spesies dan 20 substat magnetik setiap satu, jadi register ialah 40 qubit. Dua proton valens dan dua neutron valens membentuk 44Ti^{44}\mathrm{Ti}, dengan 4,000 penentu yang dibenarkan-simetri — kira-kira enam kali ganda asas 20Ne^{20}\mathrm{Ne}, menggunakan 40 qubit dan bukannya 24.

Ini ialah yang lebih besar daripada dua contoh yang boleh diselesaikan secara tepat oleh notebook, jadi anda boleh membandingkan hasil perkakasan dengan rujukan tepat.

large_scale_verified = sqd_run("gxpf1.snt", 2, 2, "44Ti", exact=True)
44Ti: 40 qubits, 2p + 2n, A = 44
2p2h pool 1602 raw -> 174 symmetry-allowed; reference energy -44.309387 MeV
packed 96 of 174 excitations; worst circuit two-qubit depth 272, 285 two-qubit gates
job dap31a02fm4c73f67dp0: 16 circuits x 10,000 shots on ibm_phoenix
160,000 shots -> 48,170 distinct bitstrings, 18.6% of shots with the right nucleon numbers
iteration 1: 187 proton x 189 neutron halves -> dimension 3891, E = -47.849086 MeV
iteration 2: 190 proton x 190 neutron halves -> dimension 4000, E = -47.876666 MeV
iteration 3: 190 proton x 190 neutron halves -> dimension 4000, E = -47.876666 MeV
reference -44.309387 MeV pooled SQD -47.876666 MeV
exact -47.876666 MeV (dimension 4000) -> 100.0% of the correlation energy

48Cr^{48}\mathrm{Cr}: melangkaui kapasiti diagonalisasi-tepat tutorial ini​

Menambah dua proton dan dua neutron menggunakan register 40-qubit yang sama (4, 4 untuk 48Cr^{48}\mathrm{Cr}) dan meningkatkan saiz asas dengan faktor kira-kira 491, kepada 1,963,461 penentu yang dibenarkan-simetri. Matriks itu jauh melangkaui apa-apa yang akan dibina tutorial ini, jadi exact=False: tiada tenaga rujukan tepat, hanya sempadan variasional dan penentu rujukan yang diperbaikinya.

Dua perkara berubah pada skala ini, dan kedua-duanya kelihatan dalam cetakan. Kolam membesar kepada beberapa ratus pengujaan yang dibenarkan, jadi belanjawan get tetap kini merangkumi minoriti daripadanya dan bukan keseluruhannya. Selain itu, subruang produk yang direntangi oleh sampel adalah lebih besar daripada MAX_DIMENSION, jadi penyelesai tumpat memangkasnya mengikut berat sampel. Sempadan itu kekal tegar tetapi mungkin kurang tepat berbanding sempadan yang dikira daripada semua konfigurasi yang disampel. Pengiraan pengeluaran akan mengekalkan sampel dan menggunakan penyelesai yang menyokong subruang yang lebih besar.

large_scale_unverified = sqd_run("gxpf1.snt", 4, 4, "48Cr", exact=False)
48Cr: 40 qubits, 4p + 4n, A = 48
2p2h pool 5536 raw -> 582 symmetry-allowed; reference energy -93.041237 MeV
packed 96 of 582 excitations; worst circuit two-qubit depth 224, 279 two-qubit gates
job dap32a02fm4c73f67eog: 16 circuits x 10,000 shots on ibm_phoenix
160,000 shots -> 55,436 distinct bitstrings, 18.6% of shots with the right nucleon numbers
iteration 1: 223 proton x 223 neutron halves -> dimension 3977, E = -96.481598 MeV
iteration 2: 223 proton x 223 neutron halves -> dimension 3977, E = -96.481598 MeV
reference -93.041237 MeV pooled SQD -96.481598 MeV
no exact reference: the symmetry-allowed basis is 1,963,461 determinants
the bound captures -3.440361 MeV of correlation energy

Nilaikan hasil tanpa rujukan tepat​

Larian 48Cr^{48}\mathrm{Cr} tiada rujukan tepat dalam tutorial ini. Gunakan sampel sedia ada untuk menilai penumpuan dan bandingkan dengan garis dasar pemilihan klasik, tanpa masa QPU tambahan atau diagonalisasi ruang-penuh.

Adakah ia sudah menumpu? Susun semula penentu yang dikekalkan mengikut beratnya dalam eigenvektor yang menumpu dan subruang menjadi bersarang, jadi mendiagonalisasi blok d×dd \times d terkemuka bagi tangga dd menjejaki penurunan sempadan merentas dua dekad saiz subruang. Jika ia masih menurun secara mendadak pada dd terbesar, had dimensi penyelesai klasik ialah kekangan pengikat dan MAX_DIMENSION ialah parameter yang perlu ditingkatkan. Jika ia telah mendatar, menambah lebih banyak penentu yang dikekalkan memberikan sedikit penambahbaikan; kemajuan selanjutnya mungkin memerlukan pensampelan konfigurasi tambahan. Hamiltonian dibina sekali sahaja pada saiz penuh dan setiap anak tangga ialah blok utama daripadanya, jadi seluruh pengimbasan hanya memerlukan satu pembinaan matriks dan bukan satu bagi setiap anak tangga.

Bagaimana pensampelan kuantum berbanding dengan pemilihan klasik? Bandingkan dengan subruang yang bersaiz sama yang dipilih oleh prosedur pemilihan klasik: ambil kolam yang disusun mengikut teori gangguan dalam susunan skor, besarkan subruang produk kepada dimensi yang sama, dan diagonalisasikannya sebaliknya. Kedua-dua lengkung adalah sempadan atas tegar bagi Hamiltonian yang sama, jadi mana-mana yang lebih rendah pada dimensi yang sama memilih penentu yang lebih baik. Perbandingan ini menentukan sama ada pensampelan perkakasan meningkatkan anggaran tenaga berbanding garis dasar klasik ini.

Subruang ini tidak dipilih untuk keadaan teruja. Pemulihan konfigurasi mengarahkan subruang menggunakan kependudukan keadaan-dasar, jadi nilai eigen yang lebih tinggi jauh lebih tidak menumpu daripada yang terendah, dan tenaga pengujaan pertama datang jauh melebihi 2+2^+ yang diukur. Untuk mencapai keadaan teruja dengan betul memerlukan subruang yang dipilih khusus untuknya.

def subspace_scaling(
inter, basis, vector, points=18, smallest=32, largest=None
):
"""Nested Rayleigh-Ritz sweep: the lowest eigenvalue of the leading d x d block, for a ladder of d.

Reordering the basis by descending weight in the converged eigenvector makes every subspace in
the ladder a subset of the next, so the energies fall monotonically and each one is a valid
variational bound. H is built once at full size; each rung is a principal block.
"""
order = np.argsort(-(np.abs(vector) ** 2))
ordered = [basis[i] for i in order]
weights = (np.abs(vector) ** 2)[order]
if (
largest is not None
): # cap the ladder so two subspaces end at a common dimension
ordered, weights = ordered[:largest], weights[:largest]
H = subspace_hamiltonian(inter, ordered)
dimensions = np.unique(
np.geomspace(smallest, len(ordered), points).astype(int)
)
rows = [
(
int(d),
float(
eigh(H[:d, :d], eigvals_only=True, subset_by_index=[0, 0])[0]
),
)
for d in dimensions
]
return rows, np.cumsum(weights)

def classical_selection(
inter, sp, reference, ranked, target, n_protons, n_neutrons
):
"""The subspace classical perturbative ranking would pick, grown to `target` dimension.

Same product construction as the sampled subspace, and the same truncation discipline -- half
configurations are offered to `grow_subspace` in order of importance and it takes as many as
fit. The only difference from the sampled path is where the ordering comes from: PT2 score
here, measured sampling weight there. So the comparison isolates *which determinants got
chosen* and nothing else.

Truncating by any other rule would not be a fair baseline. Slicing an arbitrarily ordered
list, for instance, keeps determinants by accident rather than by importance and makes the
classical subspace look worse than classical selection really is.
"""
p_ref = tuple(i for i in reference if sp[i].tz == -1)
n_ref = tuple(i for i in reference if sp[i].tz == +1)
reached = {reference}
proton_order, neutron_order = [p_ref], [n_ref]
seen_p, seen_n = {p_ref}, {n_ref}
product_budget = 4 * target

for op, _, _ in ranked: # ranked is already in descending PT2 score
h1, h2, v1, v2 = op
fresh = {
tuple(sorted(set(d) - {h1, h2} | {v1, v2}))
for d in reached
if {h1, h2} <= set(d) and not {v1, v2} & set(d)
}
reached |= fresh
for det in fresh: # first appearance fixes a half's rank
half_p = tuple(i for i in det if sp[i].tz == -1)
half_n = tuple(i for i in det if sp[i].tz == +1)
if half_p not in seen_p:
seen_p.add(half_p)
proton_order.append(half_p)
if half_n not in seen_n:
seen_n.add(half_n)
neutron_order.append(half_n)
if len(proton_order) * len(neutron_order) > product_budget:
# Half-configuration products over-count the subspace, because only the
# symmetry-allowed ones survive `product_subspace`. Stopping on the product
# count alone can therefore leave the basis far short of `target`, so check
# the dimension actually realized and widen the budget if it falls short.
trial, _, _ = grow_subspace(
sp,
[p_ref],
[n_ref],
proton_order,
neutron_order,
n_protons,
n_neutrons,
max_dimension=target,
)
if len(trial) >= target:
break
product_budget *= 2

basis, _, _ = grow_subspace(
sp,
[p_ref],
[n_ref],
proton_order,
neutron_order,
n_protons,
n_neutrons,
max_dimension=target,
)
return basis

run = large_scale_unverified
if "basis" not in run:
raise RuntimeError(
"this cell needs the subspace and eigenvector that sqd_run now returns; "
"re-run the sqd_run definition and the 48Cr cell"
)

print(
f"{run['name']}: sweeping nested subspaces of the sampled basis "
f"(dimension {run['dimension']})"
)
sampled_rows, cumulative = subspace_scaling(
run["interaction"], run["basis"], run["vector"]
)

print(
f"{run['name']}: building the classically selected subspace at the same dimension"
)
classical_basis = classical_selection(
run["interaction"],
run["states"],
run["reference"],
run["ranked"],
run["dimension"],
4,
4,
)
# Both subspaces must be scored at the same dimension. Symmetry filtering can still leave
# the classical construction short of the target when the ranked pool runs out, so take the
# dimension both actually reach, cap both ladders there, and verify they agree.
common_dim = min(sampled_rows[-1][0], len(classical_basis))
if common_dim < sampled_rows[-1][0]:
sampled_rows, _ = subspace_scaling(
run["interaction"], run["basis"], run["vector"], largest=common_dim
)
classical_rows, _ = subspace_scaling(
run["interaction"],
classical_basis,
ground_state(run["interaction"], classical_basis)[1],
largest=common_dim,
)
if sampled_rows[-1][0] != classical_rows[-1][0]:
raise RuntimeError(
f"comparison dimensions differ: sampled {sampled_rows[-1][0]}, "
f"classical {classical_rows[-1][0]}"
)

advantage = sampled_rows[-1][1] - classical_rows[-1][1]
direction = "lower" if advantage < 0 else "higher"
verdict = "beats" if advantage < 0 else "does not beat"
descent = next(
e for d, e in reversed(sampled_rows) if d <= sampled_rows[-1][0] / 2
)
for fraction in (0.90, 0.99):
count = int(np.searchsorted(cumulative, fraction) + 1)
print(
f" {fraction:.0%} of the eigenvector norm sits on {count} determinants "
f"({count / run['full_dim']:.1e} of the {run['full_dim']:,}-determinant space)"
)
print(
f" bound still falling {1000 * (sampled_rows[-1][1] - descent):+.1f} keV "
f"over the last doubling of dimension"
)
print(
f" sampled {sampled_rows[-1][1]:.6f} MeV vs classically selected "
f"{classical_rows[-1][1]:.6f} MeV at a verified common dimension of "
f"{classical_rows[-1][0]:,}"
)
print(
f" -> the sampled subspace is {abs(advantage) * 1000:.0f} keV {direction}"
)
48Cr: sweeping nested subspaces of the sampled basis (dimension 3977)
48Cr: building the classically selected subspace at the same dimension
90% of the eigenvector norm sits on 107 determinants (5.4e-05 of the 1,963,461-determinant space)
99% of the eigenvector norm sits on 593 determinants (3.0e-04 of the 1,963,461-determinant space)
bound still falling -15.6 keV over the last doubling of dimension
sampled -96.481598 MeV vs classically selected -95.314510 MeV at a verified common dimension of 3,957
-> the sampled subspace is 1167 keV lower
fig, axes = plt.subplots(1, 2, figsize=(11.2, 4.0), facecolor=SURFACE)

# left: two nested convergence curves on the same axes
ax = axes[0]
ax.set_facecolor(SURFACE)
ax.plot(
[d for d, _ in sampled_rows],
[e for _, e in sampled_rows],
"-o",
color=SERIES,
linewidth=2,
markersize=5,
markeredgecolor=SURFACE,
markeredgewidth=1,
zorder=4,
label="sampled on the QPU",
)
ax.plot(
[d for d, _ in classical_rows],
[e for _, e in classical_rows],
"--s",
color=MUTED,
linewidth=1.6,
markersize=4,
markeredgecolor=SURFACE,
markeredgewidth=1,
zorder=3,
label="classically selected, same size",
)
ax.axhline(run["e_ref"], color=RULE, linestyle=":", linewidth=1.2, zorder=1)
ax.annotate(
f"reference determinant {run['e_ref']:.2f} MeV".replace("-", "\u2212"),
xy=(sampled_rows[-1][0], run["e_ref"]),
xytext=(-2, 4),
textcoords="offset points",
ha="right",
va="bottom",
color=MUTED,
fontsize=8,
)

# mark the gap between the two curves at the largest dimension, not either curve alone
edge = sampled_rows[-1][0]
ax.plot(
[edge, edge],
[classical_rows[-1][1], sampled_rows[-1][1]],
"-",
color=SERIES,
linewidth=1.0,
alpha=0.7,
zorder=2,
)
ax.annotate(
f"{abs(advantage) * 1000:.0f} keV {direction}\nat equal dimension",
xy=(edge, 0.5 * (classical_rows[-1][1] + sampled_rows[-1][1])),
xytext=(-8, 0),
textcoords="offset points",
ha="right",
va="center",
color=SERIES,
fontsize=8.5,
)

ax.set_xscale("log")
ax.set_xlim(sampled_rows[0][0] * 0.75, edge * 1.5)
ax.set_xlabel("subspace dimension", color=MUTED)
ax.set_ylabel("variational upper bound (MeV)", color=MUTED)
ax.set_title(
f"{pretty[run['name']]}: the bound, and the subspace it {verdict}",
color=INK,
fontsize=11,
loc="left",
pad=10,
)
legend = ax.legend(frameon=False, fontsize=8.5, loc="lower left")
for text in legend.get_texts():
text.set_color(MUTED)

# right: why a few thousand determinants can bound two million
ax = axes[1]
ax.set_facecolor(SURFACE)
ranks = np.arange(1, len(cumulative) + 1)
ax.plot(ranks, 100 * cumulative, "-", color=DEEP, linewidth=2, zorder=3)
for fraction, style, label_y in ((0.90, ":", 46), (0.99, "--", 24)):
count = int(np.searchsorted(cumulative, fraction) + 1)
ax.axvline(count, color=MUTED, linestyle=style, linewidth=1, zorder=1)
ax.annotate(
f"{fraction:.0%} of the norm\non {count} determinants",
xy=(count, label_y),
xytext=(7, 0),
textcoords="offset points",
ha="left",
va="center",
color=MUTED,
fontsize=8.5,
)
ax.set_xscale("log")
ax.set_xlim(0.8, len(cumulative) * 2.6)
ax.set_ylim(0, 104)
ax.set_xlabel("determinants, ordered by weight", color=MUTED)
ax.set_ylabel("cumulative share of the eigenvector (%)", color=MUTED)
ax.set_title(
f"Sparsity: {run['full_dim']:,} determinants in the sector",
color=INK,
fontsize=11,
loc="left",
pad=10,
)

for ax in axes:
ax.grid(axis="y", color=RULE, alpha=0.55)
ax.tick_params(colors=MUTED, labelsize=8.5)
for side in ("top", "right"):
ax.spines[side].set_visible(False)
for side in ("bottom", "left"):
ax.spines[side].set_color(MUTED)
fig.tight_layout()
plt.show()

Output of the previous code cell

Bandingkan tiga larian​

Tenaga mutlak tidak boleh dibandingkan merentas nukleus dan interaksi yang berbeza, jadi fokus pada pecahan tenaga korelasi yang dipulihkan merentas larian, di mana rujukan tepat tersedia. Bandingkan juga kedalaman litar dan pecahan shot yang dibuang.

runs = [small_scale, large_scale_verified, large_scale_unverified]

print(
f"{'run':>6} {'qubits':>6} {'pool':>9} {'2q depth':>8} {'2q gates':>8} "
f"{'shots kept':>10} {'dim':>6} {'of':>9} {'% corr':>7}"
)
for r in runs:
fraction = (
"--"
if r["e_exact"] is None
else f"{100 * (r['e_sqd'] - r['e_ref']) / (r['e_exact'] - r['e_ref']):.1f}%"
)
coverage = "{}/{}".format(r["packed"], r["pool"])
print(
f"{r['name']:>6} {r['qubits']:>6} {coverage:>9} "
f"{r['two_qubit']:>8} {r['two_qubit_gates']:>8} {r['survival']:>9.1%} "
f"{r['dimension']:>6} {(r['full_dim'] or 0):>9,} {fraction:>7}"
)

print()
for r in runs:
exact = (
f"exact {r['e_exact']:11.6f}"
if r["e_exact"] is not None
else "exact unavailable"
)
print(
f"{r['name']:>6} reference {r['e_ref']:11.6f} pooled SQD {r['e_sqd']:11.6f} {exact} MeV"
)
run qubits pool 2q depth 2q gates shots kept dim of % corr
20Ne 24 78/78 228 234 31.5% 640 640 100.0%
44Ti 40 96/174 272 285 18.6% 4000 4,000 100.0%
48Cr 40 96/582 224 279 18.6% 3977 1,963,461 --

20Ne reference -29.765549 pooled SQD -40.472331 exact -40.472331 MeV
44Ti reference -44.309387 pooled SQD -47.876666 exact -47.876666 MeV
48Cr reference -93.041237 pooled SQD -96.481598 exact unavailable MeV
# Left: how much of the correlation energy was recovered, where the exact answer is known.
# Right: the bound itself for the run that has nothing to score against.
scored = [r for r in runs if r["e_exact"] is not None]

fig, ax = plt.subplots(figsize=(6.4, 3.9), facecolor=SURFACE)
ax.set_facecolor(SURFACE)
labels = [
f"{pretty[r['name']]}\n{r['qubits']} qubits\n{r['full_dim']:,} determinants"
for r in scored
]
fractions = [
100 * (r["e_sqd"] - r["e_ref"]) / (r["e_exact"] - r["e_ref"])
for r in scored
]
shades = [SERIES, DEEP, PURPLE]
bars = ax.bar(
labels, fractions, width=0.46, color=shades[: len(scored)], zorder=3
)
for bar, fraction, r in zip(bars, fractions, scored):
ax.annotate(
f"{fraction:.1f}%",
xy=(bar.get_x() + bar.get_width() / 2, fraction),
xytext=(0, 5),
textcoords="offset points",
ha="center",
va="bottom",
color=INK,
fontsize=10,
)
ax.annotate(
f"dim {r['dimension']:,}",
xy=(bar.get_x() + bar.get_width() / 2, 3),
ha="center",
va="bottom",
color=SURFACE,
fontsize=8.5,
)
ax.axhline(100, color=MUTED, linestyle="--", linewidth=1, zorder=1)
ax.annotate(
"exact diagonalization",
xy=(-0.45, 100),
xytext=(0, 4),
textcoords="offset points",
ha="left",
va="bottom",
color=MUTED,
fontsize=8.5,
)
ax.set_ylim(0, 118)
ax.set_ylabel("correlation energy recovered (%)", color=MUTED)
ax.set_title(
f"Where the exact answer is known ({backend.name})",
color=INK,
fontsize=11,
loc="left",
pad=10,
)
ax.grid(axis="y", color=RULE, alpha=0.55)
ax.tick_params(colors=MUTED, labelsize=8.5)
for side in ("top", "right"):
ax.spines[side].set_visible(False)
for side in ("bottom", "left"):
ax.spines[side].set_color(MUTED)
fig.tight_layout()
plt.show()

# the same convergence view as the walkthrough, for the run with no exact reference
convergence_plot(
large_scale_unverified["history"],
large_scale_unverified["e_ref"],
None,
f"{pretty[large_scale_unverified['name']]}: "
f"{large_scale_unverified['full_dim']:,} determinants, no exact answer to score against\n"
f"({backend.name}, {N_CIRCUITS} circuits x {SHOTS:,} shots)",
colour=DEEP,
)
plt.show()

Output of the previous code cell

Output of the previous code cell

Ringkasan​

Satu aliran kerja, tidak berubah kecuali inputnya, dijalankan pada QPU pada tiga saiz masalah: masalah 24-qubit yang boleh anda semak secara tepat, masalah 40-qubit yang masih boleh anda semak secara tepat, dan masalah 40-qubit dengan hampir dua juta keadaan asas yang melangkaui kapasiti diagonalisasi-tepat tutorial ini.

Tiga larian ini menggambarkan perkara berikut:

  • Langkah kuantum hanya perlu mencadangkan penentu. Litar adalah tetap, disemai daripada teori gangguan peringkat kedua, dan tidak pernah dioptimumkan. Tiada apa-apa dalam aliran kerja yang memerlukan amplitudnya tepat, hanya sokongannya perlu berguna. Diagonalisasi klasik dalam subruang terpilih memberikan sempadan atas variasional, walaupun sempadan berubah mengikut konfigurasi yang disampel.

  • Pengujaan qubit mengurangkan kedalaman litar. Kerana hanya sokongan yang penting, blok pengujaan fermionik boleh digantikan dengan pengujaan qubit, yang kosnya tidak bertambah dengan jarak antara orbital yang dihubungkannya. Langkah 2 mengukur penjimatan itu pada backend sebenar, iaitu perbezaan antara litar yang muat selesa di dalam koheren dan yang tidak.

  • Pemulihan konfigurasi menggunakan semula sampel berhingar. Setiap shot dengan bilangan proton atau neutron yang salah dibaiki menggunakan anggaran kependudukan semasa dan bukannya dibuang, dan setiap separuh-konfigurasi yang dibaiki boleh menambah konfigurasi kepada subruang. Meluaskan subruang variasional tidak boleh menaikkan nilai eigen terendahnya. Tutorial ini menunjukkan pemulihan konfigurasi menggunakan sampel perkakasan.

  • Kekangan pengikat berubah apabila anda meningkatkan skala. Pada 24 qubit, ansatz boleh mencapai jawapan tepat, dan hanya pensampelan yang menghalang. Pada 40 qubit dengan empat nukleon valens bagi setiap spesies, belanjawan get hanya merangkumi minoriti kolam dan penyelesai klasik tumpat mengehadkan subruang. Mengetahui yang mana daripada tiga-tiga itu yang menghadkan anda ialah kemahiran praktikal yang diajar aliran kerja ini.

Langkah seterusnya​

Cadangan

Terokai sumber berkaitan ini:

Lanjutan untuk dipertimbangkan​

  • Gantikan penyelesai tumpat. MAX_DIMENSION ialah siling untuk segalanya pada skala 48Cr^{48}\mathrm{Cr}, dan np.linalg.eigh pada matriks tumpat adalah sebabnya. Membina Hamiltonian terunjur yang sama sebagai matriks jarang dan menggunakan penyelesai eigen berulang seperti scipy.sparse.linalg.eigsh, atau penyelesai Davidson atau selected-CI yang direka untuk interaksi dua-badan nuklear, boleh menyokong subruang yang lebih besar. Had praktikal bergantung pada kejarangan matriks, memori yang tersedia, dan penumpuan penyelesai, dan tutorial ini tidak menanda aras lanjutan itu. qiskit_addon_sqd.fermion.solve_sci addon SQD bukan pengganti terus: ia membungkus penyelesai struktur-elektronik dan mengharapkan kamiran satu- dan dua-badan dalam bentuk itu, jadi struktur produk proton ×\times neutron yang dikongsi tidak mencukupi dengan sendirinya. Menggunakannya bermakna memetakan interaksi model-kelompok Persamaan (1) kepada kamiran itu dan mengesahkan hasilnya terhadap tenaga tepat yang sudah dikira oleh notebook ini.

  • Tambah batching dan subpensampelan. Aliran kerja SQD berkumpul yang diterbitkan mendiagonalisasi beberapa subsampel bebas setiap iterasi dan mengekalkan yang terbaik. Tutorial ini menggunakan satu batch setiap iterasi, yang tidak memudaratkan sempadan variasional tetapi tidak memberikan maklumat varians yang menandakan sama ada lebih banyak shot akan membantu.

  • Keadaan teruja dan sektor lain. Nilai eigen yang lebih tinggi bagi setiap Hamiltonian subruang ialah sempadan atas untuk keadaan teruja dalam sektor simetri yang sama, dan menjalankan pada MJ≠0M_J \neq 0 mencapai sektor lain. Semakan 2+2^+ dalam Langkah 1 sudah menjadi separuh daripada pengiraan ini.

  • Ruang model silang-kelompok. Pariti secara automatik dipenuhi di dalam satu kelompok major, itulah sebabnya ia tidak melakukan apa-apa di sini. Ruang sdsd-pfpf mencampurkan pariti ℓ\ell, menjadikan pariti kekangan keempat yang sebenar, satu yang tidak akan dikesan sama ada oleh pembaikan berat-Hamming SQD mahupun pembinaan produk itu sendiri.

  • Nukleus jisim-ganjil. reference_determinant memerlukan bilangan valens genap bagi setiap spesies, kerana pengisian berpasangan terbalik-masa itulah yang memaksa MJ=0M_J = 0. Nukleus ganjil memerlukan sasaran MJM_J separuh-integer dan rujukan tidak berpasangan.

Lampiran​

Bahagian ini menerangkan sebab di sebalik pembantu yang diperkenalkan dalam bahagian Persediaan.

Mengapa penskalaan semula pergantungan-jisim tidak pilihan​

Interaksi model-kelompok empirikal dipadankan pada satu jisim dan diaplikasikan merentas rantaian isotop, dengan unsur matriks dua-badan diskalakan sebagai (A/Aref)p(A/A_{\mathrm{ref}})^{p}. Kedua-dua fail interaksi membawa p=−0.3p = -0.3, dengan Aref=18A_{\mathrm{ref}} = 18 bagi keluarga USD dan 4242 bagi GXPF1. Pada baris pengepala dua-badan sesebuah fail .snt, kedua-dua nombor itu terletak di tempat frekuensi osilator dan tenaga teras secara munasabah boleh berada, yang menjadikannya mudah disalah baca; membaca eksponen sebagai tenaga teras malar menambah ofset palsu kepada setiap unsur pepenjuru dan menggugurkan penskalaan semula, mengubah tenaga korelasi sebanyak beberapa peratus. Semakan simetri dalam Langkah 1 tidak dengan sendirinya mengesahkan skala tenaga. Membandingkan tenaga pengujaan 2+2^+, yang diukur dalam MeV, dengan eksperimen memberikan semakan tambahan pada penskalaan semula bergantung-jisim. Tenaga pengujaan ialah perbezaan antara paras, jadi ia tidak mengesan ofset malar yang dikenakan pada semua tenaga.

Mengapa rujukan dicari melalui carian dan bukan melalui pengisian​

Rujukan yang jelas ialah penentu yang mengisi tenaga zarah-tunggal terendah. Ia bukan penentu tenaga-terendah, kerana pepenjuru Persamaan (1) mengandungi terma dua-badan ∑i<j⟨ij∥ij⟩\sum_{i<j} \langle ij \| ij \rangle, dan interaksi pemasangan sangat memilih menduduki rakan terbalik-masa (+mj,−mj)(+m_j, -m_j) pada ∣mj∣|m_j| tersedia yang terbesar. Dalam kelompok sdsd, itu ialah perbezaan antara pasangan mj=±1/2m_j = \pm 1/2 dan pasangan mj=±5/2m_j = \pm 5/2 bagi 0d5/20d_{5/2}, dan bernilai kira-kira 1 MeV; dalam kelompok pfpf, ia bernilai lebih hampir kepada 2. Kerana tenaga rujukan menentukan sifar bagi metrik "tenaga korelasi yang dipulihkan", pilihan yang buruk menggelembungkan metrik itu dan memberikan titik mula yang kurang tepat.

Mengehadkan kepada pengisian berpasangan menjadikan carian menyeluruh murah, dengan (npairsk)\binom{n_{\mathrm{pairs}}}{k} calon bagi setiap spesies (paling banyak beberapa ribu), dan memastikan MJ=0M_J = 0. Dalam setiap kes dalam tutorial ini yang boleh disemak terhadap penyenaraian penuh, carian mengembalikan penentu pepenjuru-terendah global, yang juga merupakan komponen tunggal terbesar bagi keadaan dasar tepat.

Mengapa amplitud peringkat-pertama, bukan sudut dua-paras yang tepat​

Mendiagonalisasi Hamiltonian 2×22 \times 2 dalam ruang {∣Φref⟩,∣α⟩}\{|\Phi_{\mathrm{ref}}\rangle, |\alpha\rangle\} memberikan sudut percampuran θexact=12arctan⁡(2V/Δ)\theta_{\mathrm{exact}} = \tfrac{1}{2}\arctan(2V/\Delta); mungkin menggoda untuk menyebut itu pilihan yang betul bagi sepasang paras yang terpencil. Dalam ansatz ini, beberapa dozen blok pengujaan bertindak secara berjujukan pada rujukan yang sama, jadi mengoptimumkan setiap blok secara berasingan tidak semestinya mengoptimumkan litar komposit.

Peranan litar menentukan pilihan sudut. Kerana ∣12arctan⁡(2x)∣≤∣x∣|\tfrac{1}{2}\arctan(2x)| \le |x| bagi setiap xx nyata, sudut tepat sentiasa lebih kecil dalam magnitud berbanding amplitud peringkat-pertama t=V/Δt = V/\Delta, dan oleh itu sentiasa meninggalkan lebih banyak amplitud pada penentu rujukan. Litar yang mengekalkan lebih banyak amplitud pada rujukan mengembalikan rujukan lebih kerap dan penentu teruja yang berbeza lebih jarang. Bagi SQD berkumpul, output berguna bagi satu shot ialah penentu yang belum dilihat oleh langkah klasik, yang mendorong penggunaan sudut yang lebih besar dalam tutorial ini. Tiada satu sudut pun perlu tepat, kerana diagonalisasi klasik membuang amplitud litar sepenuhnya dan mendapatkan semula amplitudnya sendiri.

Mengapa SQD berkumpul boleh menggunakan pengujaan qubit​

Pengujaan fermionik T=av1†av2†ah2ah1T = a_{v_1}^\dagger a_{v_2}^\dagger a_{h_2} a_{h_1} dipetakan di bawah Jordan-Wigner kepada lapan rentetan Pauli, setiap satu membawa operator ZZ pada setiap qubit antara indeks terluar. Rentetan itu mengekod tanda fermionik, dan kosnya bertambah dengan jangkauan, yang bagi pengujaan proton-neutron ialah keseluruhan register.

Menghapuskannya memberikan operator pengujaan-qubit Yordanov et al. [5]. Ia ialah operator yang berbeza: keadaan yang disediakannya berbeza daripada yang fermionik dalam tanda amplitudnya, dan kedua-dua taburan pensampelan boleh berbeza dengan ketara. Apa yang tidak berubah ialah penentu mana yang mempunyai amplitud bukan sifar, kerana setiap blok masih berputar dalam ruang dua-dimensi yang sama {∣d⟩,∣d′⟩}\{|d\rangle, |d'\rangle\} bagi setiap penentu dd yang ditindaknya, dan ia masih mengekalkan kedua-dua bilangan nukleon, MJM_J, dan pariti dengan tepat. Set penentu yang boleh dicapai adalah oleh itu sama, dan set yang boleh dicapai itulah satu-satunya perkara yang digunakan SQD berkumpul; diagonalisasi klasik menetapkan amplitudnya sendiri tanpa mengira. Langkah 2 mengesahkan dakwaan sokongan-sama pada operator sebenar daripada kolam dan mengukur apa yang dijimatkan oleh penggantian itu.

Batasannya ialah berat pensampelan berbeza, jadi kedua-dua pembinaan tidak akan menemui penentu dalam susunan yang sama pada bilangan shot yang terhingga. Kerana kedudukan yang menentukan pengujaan mana yang memasuki litar adalah klasik dan tidak berubah, dan langkah klasik itu memberi berat semula segalanya, perbezaan dalam berat pensampelan ialah pertukaran demi kedalaman litar yang berkurangan.

Mengapa MJM_J kepunyaan peringkat produk​

Pasca-pemilihan dan pemulihan konfigurasi kedua-duanya bertindak pada berat Hamming: bilangan proton dalam satu separuh register, dan bilangan neutron dalam separuh yang lain. MJ=Mp+MnM_J = M_p + M_n bukan dalam bentuk itu. Ia ialah ciri konfigurasi proton dipasangkan dengan konfigurasi neutron. Satu shot yang separuh protonnya dan separuh neutronnya masing-masing membawa bilangan nukleon yang betul mengandungi dua separuh-konfigurasi yang boleh digunakan walaupun nilai MJM_J mereka tidak saling membatalkan, kerana separuh proton pada Mp=+1M_p = +1 adalah sempurna baik sebaik sahaja dipasangkan dengan separuh neutron pada Mn=−1M_n = -1. Menapis keseluruhan shot berdasarkan MJM_J jumlah membuang kedua-dua separuh itu, dan mengenakan MJM_J pada produk yang digabung semula mengekalkannya. Hujah yang sama menerangkan mengapa recover_configurations tidak memerlukan sebarang konsep MJM_J untuk berguna dalam kes ini.

Rujukan​

  1. J. Robledo-Moreno, M. Motta, H. Haas, et al., "Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer", Science Advances 11, eadu9991 (2025). arXiv:2405.05068

  2. B. A. Brown and W. A. Richter, "New USD Hamiltonians for the sd shell", Physical Review C 74, 034315 (2006). The embedded usda.snt file carries the USDA parameters as tabulated by W. A. Richter, S. Mkhize and B. A. Brown, "sd-shell observables for the USDA and USDB Hamiltonians", Physical Review C 78, 064302 (2008).

  3. M. Honma, T. Otsuka, B. A. Brown and T. Mizusaki, "Effective interaction for pf-shell nuclei", Physical Review C 65, 061301(R) (2002).

  4. B. Huron, J. P. Malrieu and P. Rancurel, "Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions", The Journal of Chemical Physics 58, 5745 (1973).

  5. Y. S. Yordanov, D. R. M. Arvidsson-Shukur and C. H. W. Barnes, "Efficient quantum circuits for quantum computational chemistry", Physical Review A 102, 062612 (2020).

  6. National Nuclear Data Center, Evaluated Nuclear Structure Data File, Brookhaven National Laboratory. Source of the measured 2+2^+ excitation energies quoted in Step 1.