化学ハミルトニアンのエネルギー推定のためのSQD

このレッスンでは、SQDを応用して分子の基底状態エネルギーを推定する。

特に、 $4$ -step Qiskitパターンのアプローチを使用して、以下のトピックについて説明します：

1. ステップ1：問題を量子回路と演算子にマッピングする
   • $N_2$ の分子ハミルトニアンを設定する。
   • 化学にヒントを得た、ハードウェアに優しいローカル・ユニタリー・クラスター・ジャストロー（LUCJ） [1] について説明する
2. ステップ2：ターゲット・ハードウェアに最適化する
   • ハードウェア実行のために、ゲート数とアンサッツのレイアウトを最適化する
3. ステップ3：ターゲット・ハードウェア上での実行
   • 最適化された回路を実際のQPUで実行し、部分空間のサンプルを生成する。
4. ステップ4：結果の後処理
   • 自己無撞着コンフィギュレーション・リカバリー・ループの導入 [2］
     • 粒子数の事前知識と最新の反復で計算された平均軌道占有率を使用して、ビット列サンプルのフルセットを後処理します。
     • 回収されたビット列から確率的にサブサンプルのバッチを作成する。
     • 各サンプル部分空間上の分子ハミルトニアンを投影し、対角化する。
     • すべてのバッチで見つかった基底状態の最小エネルギーを保存し、平均軌道占有率を更新する。

レッスンではいくつかのソフトウェアを使用します。

• PySCF で分子を定義し、ハミルトニアンを設定する。
• ffsim パッケージを使ってLUCJアサッツを構築する。
• Qiskit をハードウェア実行用にトランスパイルする。
• Qiskit IBM Runtime QPU上で回路を実行し、サンプルを収集する。
• Qiskit addon SQD 部分空間射影と行列対角化を用いた構成回復と基底状態エネルギー推定。

1. 問題を量子回路と演算子に写像する

分子ハミルトニアン

分子ハミルトニアンは一般的な形をとる：

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ は、 $p$ -番目の基底セット要素とスピン $\sigma$ に関連するフェルミオン生成/消滅演算子である。 $h_{pr}$ と $(pr|qs)$ は、1体と2体の電子積分である。 pySCF, を使って分子を定義し、基底セット 6-31g のハミルトニアンの1体積分と2体積分を計算する。

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

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

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

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

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

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

Output:

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

このレッスンでは、ヨルダン・ウィグナー（JW）変換を使ってフェルミオン波動関数を量子ビット波動関数にマッピングし、量子回路を使って準備できるようにする。 JW変換は、M個の空間軌道のフェルミオンのフォック空間を 2M 量子ビットのヒルベルト空間にマッピングする。つまり、空間軌道は2つのスピン軌道に分割され、1つはスピンアップ( $\alpha$ )電子に関連付けられ、もう1つはスピンダウン( $\beta$ )電子に関連付けられる。スピン軌道は占有軌道と非占有軌道がある。 通常、軌道の数に言及する場合、 空間軌道の数を使用する。 スピン軌道の数は2倍になる。 量子回路では、各スピン軌道を1量子ビットで表現する。 従って、1組の量子ビットはスピンアップ軌道（ $\alpha$ ）を表し、もう1組の量子ビットはスピンダウン軌道（ $\beta$ ）を表す。 例えば、 6-31g 基底セットの $N_2$ 分子は、 $16$ 空間軌道を持つ（つまり、 $16$ $\alpha$ + $16$ $\beta$ = $32$ スピン軌道）。 したがって、 $32$ -qubit量子回路が必要になる（後述するように、アンシラ量子ビットが追加で必要になるかもしれない）。 量子ビットは、電子配置または（スレーター）行列式を表すビット列を生成するために、計算ベースで測定される。 このレッスンでは、ビット列、コンフィギュレーション、行列式という用語を同じ意味で使う。 ビット列は、スピン軌道における電子の占有率を示している。ビット位置の $1$ は、対応するスピン軌道が占有されていることを意味し、 $0$ は、スピン軌道が空であることを意味する。 電子構造問題は粒子保存的であるため、決まった数のスピン軌道だけが占有されていなければならない。 $N_2$ 分子は $5$ スピンアップ電子 ( $\alpha$ ) と $5$ スピンダウン電子 ( $\beta$ ) を持つ。 したがって、 $\alpha$ と $\beta$ 軌道を表すビット列は、 $N_2$ 分子に対してそれぞれ5つの $1\text{s}$ を持たなければならない。

1.1 量子回路によるサンプル生成：LUCJアンザッツ

このレッスンでは、量子状態の準備とそれに続くサンプリングにLUCJ(Local Unitary Coupled Cluster Jastrow) ◆[1] ansatzを使用します。 最初に、完全なUCJアサッツの様々な構成要素と、その局所バージョンで行われる近似について説明する。 次に、ffsimパッケージを使用して、LUCJのansatzを構築し、Qiskitトランスパイラを使用してハードウェア実行用に最適化する。

UCJのansatzは次のような形をしている（ $L$ レイヤーまたはUCJ演算子の繰り返しの積の場合）

ここで、 $\vert \Phi_{0} \rangle$ は参照状態であり、一般的にはハートリーフォック（HF）状態とされる。 ハートリーフォック準位は最低軌道が占有されていると定義されるため、HF準位の準備にはXゲートを適用し、占有軌道に対応する量子ビットを1に設定する必要がある。 例えば、4つの空間軌道と2つのアップ・スピンと2つのダウン・スピンのHF状態準備ブロックは以下のようになる：

8つの量子ビット（4つは「アルファ軌道」、4つは「ベータ軌道」と呼ばれる）を示す

UCJ演算子 ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ の1回の繰り返しは、軌道回転( $e^{K^{(\mu)}}$ と $e^{-K^{(\mu)}}$ )に挟まれた対角クーロン進化( $e^{iJ^{(\mu)}}$ )で構成される。

軌道回転ブロックは、単一のスピン種（ $\alpha$ （アップスピン）/ $\beta$ （ダウンスピン））に対して機能する。 各電子種について、軌道回転は、1量子ビットの $R_{z}$ ゲートのレイヤーと、それに続く2量子ビットのジブンの回転ゲート（ $XX + YY$ ゲート）のシーケンスで構成される。

2量子ビットゲートは隣接するスピン軌道（最近接量子ビット）に作用するため、SWAPゲートを必要とせず、 IBM® QPUで実装可能である。

4つのアルファ軌道量子ビットと4つのベータ軌道量子ビットを示す回路

$e^{iJ^{(\mu)}}$ は対角クーロン演算子としても知られ、3つのブロックから構成される。 そのうち2つは同じスピンセクター（ $e^{iJ_{\alpha \alpha}^{(\mu)}}$ と $e^{iJ_{\beta \beta}^{(\mu)}}$ ）で、1つは2つのスピンセクターの間で働く（ $e^{iJ_{\alpha \beta}^{(\mu)}}$ ）。

$e^{iJ^{(\mu)}}$ のすべてのブロックは、数-数ゲート $U_{nn}(\phi)$ [1] で構成されている。 $U_{nn}(\phi)$ ゲートはさらに、 $R_{ZZ}(\frac{\phi}{2})$ ゲートに続いて、2つの別々の量子ビットに作用する2つの単一量子ビット $Rz(-\frac{\phi}{2})$ ゲートに分解することができる。

同スピン・コンポーネント（ $J_{\alpha \alpha}$ と $J_{\beta \beta}$ ）は、すべての可能な量子ビットのペア間に $U_{nn}$ ゲートを持つ。 しかし、超伝導QPUは接続性に制約があるため、隣接しない量子ビット間のゲートを実現するには量子ビットをスワップしなければならない。

例えば、 $N = 4$ 空間軌道の $e^{iJ_{\alpha \alpha}^{(\mu)}}$ （または $e^{iJ_{\beta \beta}^{(\mu)}}$ ）ブロックを考えてみよう。 直線的な量子ビットの接続性の場合、最後の3つのゲートは隣接しない量子ビット間で動作するため、直接実装することはできない（例えば、 Q0 と Q2 は直接接続されていない）。 従って、隣接させるためにはSWAPゲートが必要である（下図は $3$ SWAPゲートの例）。

次に、 $J_{\alpha \beta}$ は、異なるスピンセクタの同じインデックスを持つ軌道間（例えば、 $0\alpha$ と $0\beta$ 間）のゲートを実装します。同様に、量子ビットがQPU上で物理的に隣接していない場合、これらのゲートもSWAPを必要とします。

以上の議論から、UCJアサッツは、非隣接量子ビット相互作用のためにSWAPゲートを必要とするため、HW実行にはいくつかのハードルがある。 UCJアサッツの局所的変形であるLUCJは、対角クーロン作用素から $U_{nn}$。

同じ電子種ブロック、 $J_{\alpha \alpha}$ と $J_{\beta \beta}$ ）において、 $U_{nn}$ ゲートのみを最近接結合と互換性を保ち、LUCJバージョンでは非隣接量子ビット間のゲートを削除する。 下図は、隣接しないゲートを取り除いた後のLUCJブロックである。

次に、異なる電子種の間で働く $J_{\alpha \beta}$ ブロックのLUCJバージョンは、デバイスのトポロジーに基づいて異なる形状を取ることができる。

ここでも、LUCJバージョンは互換性のないゲートを取り除く。 下図は、グリッド、ヘキサゴン、ヘビーヘックス、リニアなど、異なる量子ビットのトポロジーに対する $J_{\alpha \beta}$ ブロックのバリエーションを示している。

• 正方形 ：すべての $\alpha$ と $\beta$ 軌道間に、SWAPなしで $U_{nn}$ ゲートを持つことができる。したがって、 $U_{nn}$ ゲートを削除する必要はない。
• ヘビーヘクス ： $\alpha$ - $\beta$ 相互作用は、 $4$ -番目のインデックスを持つ（0番目、4番目、8番目など）スピン軌道ごとに保持され、 アンシラ媒介されます。つまり、 $\alpha$ と $\beta$ 軌道を表す線形鎖の間にアンシラ量子ビットが必要です。 この配置は、限られた数のSWAPを必要とする。
• 六角形 ：0番目、2番目、4番目のインデックスを付けた軌道など、他のすべての軌道は、 $\alpha$ と $\beta$ が隣接する2つの直線鎖に並べられたとき、最近接になる。
• リニア ：1つの $\alpha$ と1つの $\beta$ 軌道だけが接続されている。つまり、 $J_{\alpha \beta}$ ブロックにはゲートが1つしかない。

さまざまな量子ビット配置の

LUCJバージョンを構築するためにUCJのansatzからゲートを削除すると、よりHW互換性が高くなるが、ansatzは表現力を失う。 従って、LUCJアナザッツを使用する場合、修正UCJ演算子の繰り返し( $L$ )が必要になる可能性がある。

1.2 LUCJ アンスァッツ初期化

LUCJはパラメータ化されたアサッツであり、ハードウェアの実行前にパラメータを初期化する必要がある。 アナザッツを初期化する一つの方法は、古典的なクラスターシングル・ダブルス結合（CCSD）法の t1 と t2 の振幅を使うことである。 t1 の振幅は単一励起作用素の係数であり、 t2 の振幅は二重励起作用素の係数である。

LUCJアサッツを t1 と t2 の振幅で初期化すると適切な結果が得られるが、アサッツのパラメータはさらなる最適化が必要かもしれないことに注意。

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

t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 LUCJアプローチの構築 ffsim

ffsim パッケージを使って、上記で計算された t1 と t2 の振幅を用いたansatzを作成し、初期化する。 私たちの分子は閉殻ハートリーフォック状態を持っているので、UCJアサッツのスピンバランス変形を使用します、 UCJOpSpinBalanced.

IBM ハードウェアはヘビーヘキストポロジーを持つので、 [1] で使用され、上で説明したジグザグパターンを量子ビット相互作用に採用する。 このパターンでは、同じスピンを持つ軌道（量子ビット）がライン・トポロジーで結ばれている（赤丸と青丸）。 ヘビーヘキソトポロジーのため、異なるスピンの軌道は4番目の軌道ごとに、つまり0番目、4番目、8番目......というようにつながっている（紫色の円）。

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

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

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

nelec = (num_elec_a, num_elec_b)

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

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

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

層が繰り返されるLUCJ解析は、隣接するいくつかのブロックを統合することで最適化できる。 n_reps=2 の場合を考えてみよう。 真ん中の2つの軌道回転ブロックは、1つの軌道回転ブロックに統合することができる。 ffsim パッケージには、このような隣接ブロックをマージして回路を最適化するパス・マネージャー（ ffsim.qiskit.PRE_INIT ）がある。

2. 対象ハードウェア向けに最適化

まず、好きなバックエンドを取得する。 バックエンド用に回路を最適化し、最適化された回路を同じバックエンドで実行し、部分空間のサンプルを生成する。

from qiskit_ibm_runtime import QiskitRuntimeService

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

次に、ansatzを最適化し、ハードウェアと互換性を持たせるために、以下のステップを推奨する。

• 上記のジグザグパターン（間にアンシラ量子ビットを挟んだ2つの直線チェーン）に従った物理量子ビット（initial_layout ）をターゲットハードウェアから選択する。 このパターンで量子ビットを配置することで、より少ないゲートで効率的なハードウェア互換回路を実現できる。
• Qiskitの generate_preset_pass_managerbackend initial_layout関数を使用してステージドパスマネージャーを生成します。
• ステージド・パス・マネージャーの pre_init ステージを ffsim.qiskit.PRE_INIT に設定する。 ffsim.qiskit.PRE_INIT には、ゲートを軌道回転に分解し、軌道回転をマージするQiskitトランスパイラ・パスが含まれており、最終的な回路に含まれるゲートの数が少なくなります。
• サーキットでパスマネージャーを実行する。

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

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

initial_layout = spin_a_layout + spin_b_layout

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

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

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

Output:

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

3. 対象ハードウェア上で実行する

ハードウェア実行向けに回路を最適化した後、ターゲットハードウェア上で実行し、基底状態のエネルギー推定のためのサンプルを収集する準備が整いました。 回路は1つしかないため、「 qiskit-ibm-runtime ジョブ実行モード 」を使用して、この回路を実行します。

from qiskit_ibm_runtime import SamplerV2 as Sampler

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

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

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

4. 結果の後処理

SQDワークフローの後処理部分は、以下の図を使って要約することができる。

LUCJアサッツを計算基底でサンプリングすると、ノイズの多いコンフィギュレーションのプール $\tilde{\mathcal{\chi}}$ が生成され、後処理ルーチンで使用される。 これには、（詳細は後述するが） コンフィギュレーション・リカバリーと呼ばれる、電子数が正しくないコンフィギュレーションを確率的に修正する方法が含まれる。 次に、正しい電子番号 $\tilde{\mathcal{\chi}}_{R}$ を持つコンフィギュレーションのみがサブサンプリングされ、各ユニークなコンフィギュレーションの出現頻度に基づいて複数のバッチに分配される。 サンプルの各バッチは部分空間( $\mathcal{S^{(k)}}$ )を定義する。次に、分子ハミルトニアン（ $H$ ）を部分空間に投影する：

投影された各ハミルトニアン $H_{\mathcal{S}^{(k)}}$、固有値と固有ベクトルを計算するために対角化され、固有状態が再構成される。 このレッスンでは、 PySCF のDavidsonの方法を対角化に使う qiskit-addon-sqd パッケージを使って、ハミルトニアンを射影し、対角化する。

次に、バッチから最も低い固有値（エネルギー）を収集し、平均軌道占有率 $\text{n}$ も計算する。平均軌道占有率情報は、ノイズ軌道を確率的に修正するための軌道復元ステップで使用される。

次に、自己無撞着構成回復ループを詳細に説明し、 $N_2$ ハミルトニアンの基底状態エネルギーを推定するために、上記のステップを実装する具体的なコード例を示す。

4.1 構成復旧：概要

ビット列（スレーター行列式）の各ビットはスピン軌道を表す。 ビット列の右半分はスピンアップ軌道を表し、左半分はスピンダウン軌道を表す。 1 はその軌道が電子によって占有されていることを意味し、 0 はその軌道が空であることを意味する。 私たちは粒子（アップスピン電子とダウンスピン電子）の正しい数を先験的に知っている。 $N_x$ 個の電子（つまり、ビット列には $N_x$ 個の $1$ s がある）を含む行列式 $x$ があるとする。 正しい粒子数は $N$ である。もし $N_x \neq N$ ならば、ビット列はノイズによって壊れていることがわかる。 自己無撞着コンフィギュレーション・ルーチンは、平均軌道占有率情報を活用して、 $|N_x - N|$ ビットを確率的に反転させることにより、ビット列の修正を試みる。 平均軌道占有率( $n$ )は、ある軌道が電子によって占有される確率を示す。 $N_x < N$ の場合、電子の数が少なくなるため、 $0$ s を $1$ s に反転させる必要がある。

反転の確率は、 i-番目のスピン軌道について $|x[i] - avg\_occupancy[i]|$。 [2] では、修正された ReLU 関数を用いた重み付き反転確率を使用している。

ここで、 $h$ は、 ReLU 関数の「コーナー」の位置を定義し、パラメータ $\delta$ は、コーナーにおける ReLU 関数の値を定義する。 $\delta = 0$ の場合、 $w$ は真の ReLU 関数となり、 $\delta >0$ の場合は修正された ReLU となる。 論文では、著者らは $\delta = 0.01$、 $h =$ アルファ（またはベータ）粒子の数/アルファ（またはベータ）スピン軌道の数 $= N/M$ （充填係数）を使用している。

平均軌道占有率( $n$ )は事前にはわからない。 基底状態推定の最初の反復は、両方のスピン種において正しい粒子数のみを持つ構成から始まります。 最初の反復の後、基底状態の推定値が得られ、その推定値を使用して、 $n$ の最初の推測を構築することができる。この推測 $n$ を使用して、コンフィギュレーションを回復し、基底状態の推定の次の反復を実行し、 $n$ の推測を自己無撞着に改良する。このプロセスは、停止基準が満たされるまで繰り返される。

$N = 2$ と $x = |1000\rangle$ ( $N_x = 1$ ) について以下の例を考えてみよう。粒子数を補正するために 0s のいずれかを 1 に反転する必要があり、その選択肢は 1100、 1010、 1001 である。 反転する確率に基づき、いずれかの選択肢が回収されたコンフィギュレーション （または正しいパーティクル数を持つビット列）として選択される。

最初の反復で2つのバッチを実行し、それらから推定された基底状態が次のとおりだとする：

計算基底状態とその振幅を用いて、スピン軌道(qubit)ごとの電子占有率( occupancy )を計算することができます(確率=｜振幅｜ $^2$ )。以下では、推定された基底状態に現れる各ビット列の量子ビットごとの占有率を表にして、一括して全軌道占有率を計算します。 なお、Qiskitの順序規則に従い、右端のビットは qubit-0 ( Q0 ) を表し、左端のビットは Q3 を表す。

占有率( Batch0 )：

稼働率 ( Batch1 )

稼働率（バッチ平均）

上記で計算した平均軌道占有率を用いて、コンフィギュレーション $x = \vert 1000 \rangle$ における異なる軌道のフリップ確率を求めることができる。 Q3。 で表される軌道はすでに占有されており、フリップする必要がないため、そのp(flip)を $0$ とする。残りの未占有の軌道については、フリップの確率はそれぞれ $\vert x[i] - \text{n}[i] \vert$。 p(flip)とともに、前述の修正 ReLU 関数を用いて、反転に関連する確率の重みも計算する。

フリップの確率 ( $x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )

最後に、上記の重み付けされた確率を使って、占有されていない Q2、 Q1、 Q0 軌道の一つを反転させることができる。 上記の値から、 Q0 が反転する可能性が最も高く、回収可能なコンフィギュレーションは $\vert \text{1001} \rangle$ となる。

完全な自己無撞着コンフィギュレーション・リカバリー・プロセスは、以下のように要約できる：

最初の反復： 量子コンピュータによって生成されたビット列（コンフィギュレーションまたはスレーター行列式）が、各スピンセクターの粒子数が正しい（ $\widetilde{\chi}_{correct}$ ）コンフィギュレーションと正しくない（ $\widetilde{\chi}_{incorrect}$ ）コンフィギュレーションの両方を含む集合 $\widetilde{\chi}$ を形成しているとします。

1. ( $\widetilde{\chi}_{correct}$ ) の構成がランダムにサンプリングされ、部分空間投影のためのベクトルのバッチ $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ が作成される。 バッチ数と各バッチのサンプル数は、ユーザー定義のパラメーターである。 各バッチのサンプル数が多ければ多いほど、部分空間の次元が大きくなり、対角化の計算量が多くなる。 一方、サンプル数が少なすぎると、基底状態のサポートベクトルを見落とし、誤った推定につながる可能性がある。
2. バッチに対して固有状態ソルバー（つまり、部分空間への射影と対角化）を実行し、近似固有状態を得る。 $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. 近似固有状態から、 $n$ の最初の推測を構築する。

その後の反復：

1. $n$ を使って、 $\widetilde{\chi}_{incorrect}$ の間違った粒子番号のコンフィギュレーションを修正する。 $\widetilde{\chi}_{correct\_new}$ とする。すると、 $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ は正しい粒子番号を持つ新しいコンフィギュレーションの集合を形成する。
2. $\widetilde{\chi}_{R}$ をサンプリングしてバッチを作成する $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$。
3. 固有状態ソルバーは新しいバッチで実行され、基底状態の新しい推定値を生成する $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$。
4. 近似的な固有状態から、 $n$ の精密な推測を構築する。
5. 停止基準が満たされない場合は、ステップ 2.1 に戻る。

4.2 基底状態推定

まず、カウントをビット列行列と確率配列に変換し、後処理を行う。

行列の各行は、一意なビット列を表す。 Qiskitでは量子ビットはビット列の右からインデックスされるので、列 0 は量子ビット N-1 を表し、列 N-1 は量子ビット 0 を表し、 N は量子ビットの数である。

アルファ軌道は列インデックス範囲 (N, N/2] （右半分）で表され、ベータ軌道は列範囲 (N/2, 0] （左半分）で表される。

from qiskit_addon_sqd.counts import counts_to_arrays

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

このテクニックにとって重要なユーザー制御オプションがいくつかある：

• iterations:自己無撞着構成回復の反復回数
• n_batches:固有状態ソルバーへのさまざまな呼び出しによって使用されるコンフィギュレーションのバッチ数
• samples_per_batch:各バッチに含まれるユニークなコンフィギュレーションの数
• max_davidson_cycles:各Eigensolverが実行するDavidsonサイクルの最大数

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

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

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

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

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

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

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

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

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

Output:

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

4.3 結果の検討

最初のプロットは、数回の反復の後、基底状態エネルギーを ~24 mH 以内で推定していることを示している（化学的精度は一般的に1 kcal/mol $\approx$ 1.6 mH ）。 2番目のプロットは、最終反復後の各空間軌道の平均占有率を示している。 スピン・アップ電子とスピン・ダウン電子の両方が、解の最初の5つの軌道を高い確率で占めていることがわかる。

推定された基底状態のエネルギーはまずまずだが、化学的精度の限界（ $\pm \approx 1.6$ mH ）には達していない。 このギャップは、上で射影と対角化に使った部分空間の次元が小さいことに起因している。 samples_per_batch=500 を使用したので、部分空間は最大 $500$ ベクトルでスパンされる。これは基底状態サポートからの欠損ベクトルである。 samples_per_batch パラメータを増やすと、より古典的な計算リソースと実行時間を犠牲にして精度が向上するはずである。

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

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

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

import matplotlib.pyplot as plt

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

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

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

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

Output:

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

読者のための練習問題

samples_per_batch （例えば、 $1000$ から $10000$ まで、 $1000$ のステップで）徐々にパラメータを増やし、推定された基底状態のエネルギーを比較する。

参照

[1] M. モッタら "相関電子状態に対する物理的直観とハードウェア効率の橋渡し：電子構造に対する局所ユニタリークラスターJastrowアサッツ" (2023) 化学だ。 科学..、 2023, 14, 11213.

[2] J. ロブレド＝モレノら、 「量子中心スパコンで厳密解を超える化学」（2024年）。 arXiv:quant-ph/2405.05068.