基底状態推定のためのSqDRIFTアルゴリズム

使用量の目安: Heron r3プロセッサーで180秒(注: これはあくまで目安です。実際の実行時間は異なる場合があります。)

学習内容

• トロッター分解と比べて、より浅い回路を作成する方法を学ぶ

• qDRIFTとSQDを使用した基底状態推定のエンドツーエンドのワークフローを順に確認する

• qiskit-fermions を他のQiskitアドオンと組み合わせて、このようなワークフローを実装する方法を学ぶ

このチュートリアルは、教育目的でPythonノートブックとして提供されています。

前提条件

• サンプルベース量子対角化(SQD)の概要を読む

• サンプルベース・クリロフ量子対角化(SKQD)のレッスンを読む

背景

SqDRIFTはSKQDの変種で、ビット列をサンプリングする元となるアンザッツを選ぶ必要をなくし、ターゲットのハミルトニアンから直接構築した時間発展回路のアンサンブルを使用します。これは、ハミルトニアンの係数に基づいてより小さな時間発展演算子をサブサンプリングすることで実現され、qDRIFTトロッター分解法として知られています。

このチュートリアルでは、Qiskit Fermionsを使用して、qDRIFTアルゴリズムのためのより自然なフェルミオン回路を作成し、続いてフェルミオンのレイアウトパスと合成パスを使用してから、ハードウェアで実行するために従来のQiskitパイプラインに回路を渡します。

ハミルトニアンを次の形とします。

$$H = \sum_{i=1}^{N} c_i h_i$$

ここで、一般性を失うことなく、$c_i > 0$ であること、および $h_i$ の最大固有値の絶対値が $1$ に等しいことを要求します。符号付きまたは複素数の係数はすべて $h_i$ に吸収されるため、係数 $c_i$ は厳密に正の重みであり、$h_i$ が各項の方向を担います。ここで $N$ はハミルトニアンの項の数(グループ化後はグループの数)であり、ハミルトニアンの性質であって、1つの回路にサンプリングされる演算子の数(以下では $n$ と表記)とは別のものです。

qDRIFTアルゴリズムは、ターゲット時間 $t$ に対して、ある演算子 $V_k$ を実現します。ここで $k$ は $1 \cdots K$ を動き、$k_{th}$ 番目のSqDRIFT回路を表し、次のように定義されます。

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

ここで $n$ は回路ごとにサンプリングされる演算子の数、$K$ はアンサンブル内の回路の数です。積は $N$ 個のハミルトニアン項すべてではなく $n$ 回の抽出にわたって取られ、項は復元抽出で引かれるため、同じ $h_i$ が1つの $V_k$ の中に複数回現れることがあります。

その量:

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

は係数の $L_1$ ノルムであり、そのため $n$ 個のステップはそれぞれ、どの項が引かれたかに関係なく同じ継続時間 $\lambda t / n$ だけ発展します。ステップ角が一様であることがqDRIFTの特徴です。係数は、その項がどれだけ回転されるかではなく、どれだけ頻繁に引かれるかによって結果に影響します。添字は次の分布からサンプリングされます。

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

したがって、系列 $(k_1, \ldots, k_n)$ は、この分布から引かれた項の添字のランダムな列です。$c_i$ は正で和が $\lambda$ になるため、これは正規化された確率分布であり、ランダムな抽出にわたる結果のチャネルの期待値は、$H$ のもとでの時間発展を近似し、その誤差は $n$ が増えるにつれて減少します。近似誤差は項の数 $N$ ではなく $\lambda$ に依存することに注意してください。

(SqDRIFTの論文では項の数を $\mathcal{N}$、系列の長さを $N$ と書いていますが、ここでは2つを明確に区別するために $N$ と $n$ を使用します。)

このチュートリアルでは、このようなランダム化された回路のアンサンブルを生成する方法を示します。これらの回路を作成した後、異なる演算子に対してクリロフ部分空間を作成する方法と同様に、異なる時間パラメータを持つ複数のそのような演算子からビット列をサンプリングします。これにより、基底状態ベクトルとサンプリングされたビット列との重なりが高くなります。

要件

このチュートリアルを始める前に、次のものがインストールされていることを確認してください。

• Python(>=3.10)の仮想環境
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0(名前が複数形であることに注意)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

必要なパッケージはすべて、次のコマンドでインストールできます。

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

セットアップ

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

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

シミュレーターの例

ステップ1: 古典的な入力を量子問題にマッピングする

FCIDumpの読み込みと準備

このチュートリアルでは、窒素(N2)の電子構造ハミルトニアンを読み込みます。フェルミオン演算子を作成する方法は他にもあります。ドキュメントqiskit_fermions.operators.libraryを参照してください。

このFCIDumpについて。 ファイル N2_sto_3g は、実験的な平衡結合長である1.09 $\AA$ の原子間距離における、最小のSTO-3G基底での窒素分子($N_2$)を記述しています。ヘッダーでは NORB=10、NELEC=14、MS2=0 が宣言されており、10個の空間軌道(したがって20個のスピン軌道で、Jordan-Wigner変換では20量子ビット)と、スピン一重項の14個の電子、すなわち7個の $\alpha$ 電子と7個の $\beta$ 電子を表します。すべての軌道には対称性ラベル1が与えられており、点群対称性は利用されていません。全空間のSTO-3Gダンプであるため、凍結される軌道はなく、相関空間は十分に小さいので、次のセルで示すように、比較用の厳密なFCI参照エネルギーを古典的に計算できます。

同等のファイルは、PySCFで再生成できます。

from pyscf import gto, scf, tools

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

積分は収束したSCF軌道に依存するため、再生成されたファイルは、同梱のものと軌道の位相や順序が異なる場合がありますが、全エネルギーには影響しません。

ファイルの入手。 FCIDumpは、このGitHubリポジトリにあります。下のセルを実行すると、チュートリアルの残りの部分が想定する場所に取得できます。

まず、pyscfが提供する cisolver を使用して参照エネルギーを取得します。これは、扱っている分子の真の基底状態エネルギーです。そのために、最初に norb と nelec を宣言します。これらはそれぞれ軌道数と電子数です。次に h1e と h2e を宣言します。これらはそれぞれ1電子積分と2電子積分です。これらはすべて、後でSQDにも使用されます。

import os
from urllib.request import urlopen

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

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

ハミルトニアンの読み込み

必要なデータが揃ったので、qiskit-fermions と互換性のある形式で、FCIファイルからハミルトニアンを読み込みます。

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

qiskit-fermions を用いたフェルミオンワークフロー

まず、フェルミオン回路に特化したトランスパイラーパスとゲートを提供する qiskit-fermions を使用して、ハミルトニアンをフェルミオン回路モデルにマッピングします。これらは後で、このワークフローにおいてQiskitの従来のトランスパイラーパスの前に使用されます。

項のグループ化

結果の再現性を確保するため、まず canonical_order を使用して、構造のみに基づいて項を並べ替えます。したがって、canon リスト内の演算子の順序は固定されます。これにより、今後使用する QDriftTrotterization パスがqDRIFT演算子を作成するためにランダムな添字をサンプリングするため、作成される演算子の再現性が確保されます。

このステップでは、同一の係数を持つ関連する項をグループ化することで、電子構造ハミルトニアンに存在する多くの対称性を利用します。これによりqDRIFTプロトコルがサンプリングする演算子の係数分布は変わりますが、収束保証には影響しません。重要なのは、対称性によって関連する項をグループ化すると、パウリ項の好ましい打ち消しが生じ、それらの作用のもとで状態を時間発展させるときの回路全体の深さが短くなることです。

qiskit-fermions は、このグループ化を行う group_terms_by_electronic_structure 関数を提供しています。

group_terms_by_electronic_structure は、項が正規順序になっていることを前提としている点に注意してください。

対角項のフィルタリング

回路の生成に使用するハミルトニアンから対角項を取り除き、$n$ 個のqDRIFTサンプリングスロットが、配置間で占有を移動させる項に使われるようにします。このような項は、次のステップで Evolution ゲートが構築される前の、この時点でハミルトニアンからフィルタリングしておくのが最適です。

対象となる項は、占有数基底で対角的な項、つまり数演算子の積 $a^\dagger_i a_i$ です。この記述に当てはまる項は3種類あります。

• 定数のエネルギーオフセット。数演算子を1つも含まない積で、その時間発展は大域位相にしか寄与しません。

• 個別の数演算子 $n_i$。その時間発展は1量子ビットの $Z$ 回転に帰着します。

• $n_i n_j$ のような高次の積。

これらはそれ自体では、占有数配置の間で占有を移動させることはなく、すでに存在する配置の位相にのみ作用します。ただし、これらは無関係ではありません。それらの相対位相は、回路の後半で励起項が生み出す干渉に影響するため、フィルタリングすると実際に生成される時間発展が変わり、サンプリング分布が変わる可能性があります。これは回路生成ステップにおける意図的な近似であり、サンプリングを励起項に集中させるために行われるもので、サンプリング分布に影響を与えないステップではありません。qDRIFTの収束保証を損なわない上記の対称性グループ化とは異なり、このフィルターは時間発展させる演算子そのものを変えます。そのため、回路はもはやハミルトニアン全体のもとでの時間発展を近似せず、qDRIFTの誤差上界は元の演算子ではなくフィルタリング後の演算子に適用されます。ここでこれが許容されるのは、回路が配置を提案するためのサンプリングのヒューリスティックにすぎないからです。フィルターは回路の構築に使うハミルトニアンにのみ適用され、後の古典的対角化では対角項を含む完全なハミルトニアンを使用するため、エネルギー推定自体からは項が失われません。SQDの精度はその古典ステップに依存し、配置がどのように提案されたかにかかわらず、その古典ステップはサンプリングされた部分空間において変分的であり続けます。

filter_diagonal_terms() 関数は、このような項を演算子からその場で取り除きます。正規順序の構造、すなわち生成モードの多重集合が消滅モードの多重集合と一致するという条件からそれらを識別するため、すでに正規順序になっている演算子に対してのみ有効です。この前提は実行時にチェックされません。

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

print(len(canon.groups))

5060

ハミルトニアンの項をグループ化できたので、回路のアンサンブルを生成するために、次のパラメータを決めます。

• 生成する回路の数: num_circuits
• 励起グループで表した各回路の長さ: num_exc
• 異なる時間発展時間のための係数: times

フェルミオン回路の作成

次に、各時間ステップのフェルミオン回路を作成します。各回路は、先ほど宣言した発展時間を持つ単一の時間発展ゲートで構成されます。時間発展演算子はハミルトニアンです。後でこれらの回路にトランスパイラーパスを実行して、qDRIFT回路を作成します。

アンザッツの準備

InitializeModes クラスを使用して、Hartree-Fock状態を準備します。窒素の場合、この処理は、最初の num_elec_a 個の量子ビットに、続いて num_elec_b 個の量子ビットにXゲートを適用するだけで、窒素ではどちらも7に等しくなります。この状態は、窒素の7個の $\alpha$ 電子と7個の $\beta$ 電子を表します。

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

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

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

ステップ2: 量子ハードウェアでの実行に向けて問題を最適化する

回路ができたので、まず qiskit-fermions で利用できるパスを使ってフェルミオンレベルの最適化を行い、続いて選択したバックエンド向けに回路をトランスパイルします。これはシミュレーターの実験なので、まずAerSimulator向けに行います。 各グループの重みの計算

このステップでは、ハミルトニアンの係数に比例する確率で項を確率的にサンプリングする qDRIFT サンプリングを行います。これは qDRIFT トランスパイラー・パスが実行してくれます。ハミルトニアンに長距離結合や2次を超える項が含まれていても、限られたQubitの接続性のもとでハードウェア上でより効率的に実行できる、より浅いCircuitを作成できるようになります。 項のグループ化の後、重みに基づいて演算子をサンプリングします。各演算子 $h_i$ について、重み $W_{h_i}$ は次のように定義されます。

$$W_{h_i} = |c_i| / \lambda$$

ステップ1で項をグループ化したため、ここでの各 $h_i$ はグループ全体を表します。$c_i$ はグループ $i$ に含まれる項の係数の絶対値の平均であり、グループ内の各項は、係数をその符号のみに縮小して時間発展させます。

フェルミオン最適化とハードウェアネイティブ最適化

関数 generate_preset_jw_pass_manager() は MultiStagePassManager を返します。これは FermionicCircuit を受け取り、ハードウェア上で動作させるためにトランスパイルできる、最適化された最終的なCircuitを生成します。デフォルトの最適化ステージを、QDriftTrotterization パスを含む FermionicPassManager に置き換えます。

• QDriftTrotterization パスは、内部で重みの計算とサンプリングを行い、サンプリングに使用するCircuitを生成します

• RelabelModes パスは、フェルミオンモードを並べ替えてQubit間の接続性を最適化し、Gateの深さを減らすために使用できる、もう1つの最適化パスです。詳細は API リファレンスを参照してください

MultiStagePassManager の残りのステージは自動的に実行され、フェルミオンからQubitへの完全なマッピングを処理します。

• F2QLayout: プリセットのパスマネージャーは TrivialF2QLayout パスを適用し、$n$ 個のフェルミオンビットを $n$ 個のQubitに自明にマッピングします。

• F2QSynth: フェルミオンベースのCircuit命令をQubitベースの命令にマッピングするトランスパイル・パスです。

qdrift = QDriftTrotterization(num_exc, rng=19)

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

フェルミオンレベルの最適化が完了したので、シミュレーター上で実行するためにCircuitをトランスパイルできます。

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

ステップ3: Qiskit プリミティブを使用して実行する

Circuitが用意できたので、AerSimulator 上で Qiskit プリミティブを使って実行できます。異なるCircuitのカウントをすべて結合します。最後に SQD で後処理する前に、それらをブール値ベクトルに変換します。

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

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

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

ステップ4: 後処理を行い、目的の古典形式で結果を返す

SQD にビット列を使用する

選択したビット列に対して対角化スキームを実行し、分子の基底状態エネルギーに対応する最低固有値を求めます。コールバック関数を作成し、初期占有数を宣言し、パラメーターを設定してから、対角化スキームを実行します。コールバック関数は、各イテレーションで現在のイテレーションと現在の固有値の推定値を出力するために使用されます。

最後に、基底状態の推定値を得るために、得られたエネルギーに nuclear_repulsion_energy を加えます。

注: ノイズのないシミュレーターであっても、部分空間の次元はイテレーション間で固定されません。各サブサンプルは異なる構成の集合を抽出し、回復ステップがイテレーションの間にプールを再形成するため、報告される次元はサブサンプルごとに変化します。ノイズなしのサンプリングだけでは、選択された部分空間の次元は固定されません。一方、ハードウェアでの実行では、ノイズを含むショットが粒子数対称性を破り、構成回復によってそれらが追加の基底ベクトルになるため、系統的により大きな部分空間が得られる傾向があります。このため、ハードウェアのセクションでは、ビット列を除外する別のステップも導入します。

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

ハードウェアの例

この例では20 Qubit（10個の空間軌道）を使用します。この選択は、すばやく実行できるチュートリアルのための便宜上のものであり、この手法の厳密な上限ではありません。

古典ステップのコストは、Qubit数によって直接決まるわけではありません。SQD は、サンプリングされた構成が張る部分空間に射影したハミルトニアンを対角化します。したがって、古典コストを左右するのは、その選択された部分空間の次元（ここでは samples_per_batch、num_batches、およびCircuitが実際に生成する異なる構成の数によって決まります）と、射影されたハミルトニアンを適用するために必要なスパース線形代数です。完全CI空間は軌道数と電子数に応じて組合せ的に増大しますが、選択された部分空間はその小さく調整可能な一部であり、そのサイズは直接制御できます。したがって、Qubit数と古典的な難しさはある程度独立に変えることができます。小さな系を非常に大きな部分空間で対角化するよりも、広い軌道空間を控えめな部分空間にサンプリングするほうが低コストになることもあります。

したがって実際には、実現可能な系のサイズは、目的の精度に必要な部分空間の次元と、固有値ソルバーが利用できるメモリおよびコア数によって決まります。より大きな軌道空間では、化学的精度に達するために通常より大きな部分空間が必要となり、それがやがて分散リソースの利用を動機づけます。このステップをスケールアウトする方法については、qiskit-addon-sqd-hpc を参照してください。固定のカットオフを仮定するのではなく、報告される部分空間の次元とイテレーションごとのエネルギー収束を観察し、エネルギーの改善が止まるか、利用可能なメモリを使い切るまで部分空間のサイズを増やすのが実用的なアプローチです。

注: ハードウェアのノイズに起因するサンプリング誤差のため、ハードウェアでの実行で対角化のために作成される部分空間は、シミュレーターを使用した場合よりも大きくなります。対角化したい部分空間の次元は増えますが、SQD のノイズに対する頑健性により、このワークフローは依然として正確な答えを与えてくれます。

疑似ビット列の除外

ここでは追加のステップを実行することもできます。Circuit実行から得られたすべてのビット列がそろったら、SQD を実行する前に無効なビット列を除外するか、除外せずに先へ進むことができます。ハードウェアでの実行では、一般に除外を行わないほうが望ましいです。対称性が破れたショットを構成回復に利用でき、それらを有効な構成に修復して、ショットを単に捨てる代わりに部分空間を広げられるためです。

窒素は $\alpha$ 電子と $\beta$ 電子をそれぞれ7個しか持てないため、出力の前半と後半で1の数が7より多いまたは少ないビット列はすべて破棄できます。ビット列が有効かどうかを確認し、有効でなければ破棄する関数を定義します。疑似ビット列を除外したら、残りを対角化スキームに渡します。下の PRUNE フラグを使用して、2つの動作を切り替えます。

除外は、最終的な部分空間を形作るいくつかの選択肢の1つにすぎないことに注意してください。ほかにも、Circuit数、時間発展の時間の集合、対角項のフィルタリングがあります。除外した実行と除外しなかった実行の比較は、ほかのすべてを固定した場合にのみ意味を持ちます。このチュートリアルの C++ 版では、回復ではなく事後選択を行い、それ以外のパラメーターも異なるため、これについてより詳しく説明しています。

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

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

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

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

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

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

次のステップ

おすすめ

この内容に興味を持たれた方は、次の資料にも関心があるかもしれません。

• フェルミオン格子モデルのサンプルベース・クリロフ量子対角化 - 変分アンザッツの代わりに時間発展Circuitを使用する関連チュートリアル。
• 化学ハミルトニアンのサンプルベース量子対角化 - 量子化学シミュレーション用の局所ユニタリー・クラスター・ジャストロウ (LUCJ) Circuitの構築方法に関するチュートリアル。
• SqDRIFT 論文 - このチュートリアルの元となった文献。（この論文で議論されている最適化の一部は現在開発中であり、使用しているライブラリーの進化に応じて、このチュートリアルは将来変更される可能性があります。）