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

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

C++版をお探しですか?

このチュートリアルではPythonを使用します。ソースコードとビルド手順を含むC++実装については、C++版SqDRIFTチュートリアルをご覧ください。

学習成果

トロッター化と比較して、より深さの小さい回路を作成する方法を学びます
qDRIFTとSQDを使用した基底状態推定のエンドツーエンドのワークフローを一通り確認します
qiskit-fermionsを他のQiskitアドオンと連携させて、このようなワークフローを実装する方法を学びます

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

前提条件

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

背景

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

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

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

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

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

qDRIFTアルゴリズムは、ターゲット時間 $t$ に対して、以下のように定義される演算子 $V_k$ を実現します。ここで $k$ は $1 \cdots K$ の範囲を取り、 $k$ 番目のSqDRIFT回路を表します。

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

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

次の量:

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

は係数の $L_1$ ノルムであり、したがって $n$ 個の各ステップは、どの項が抽出されたかにかかわらず、同じ時間 $\lambda t / n$ だけ発展します。ステップ角の均一性はqDRIFTの特徴的な性質です。係数は、その項がどれだけ回転されるかではなく、その項がどれだけの頻度で抽出されるかを通じて結果に影響を与えます。インデックスは次の分布からサンプリングされます。

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

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

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

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

要件

このチュートリアルを開始する前に、以下がインストールされていることを確認してください

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 — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

# 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個の空間軌道(したがってJordan-Wigner変換の下で20個のスピン軌道、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 を宣言します。次に、1電子積分と2電子積分をそれぞれ表す h1e と h2e を宣言します。これらはすべて、後でSQDにも使用されます。

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "assets/sqdrift/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 = "assets/sqdrift/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種類あります。

定数エネルギーオフセット、すなわちゼロ個の数演算子の積で、その時間発展は大域的な位相のみに寄与するもの。
個々の数演算子 $n_i$ で、その時間発展は単一量子ビットの $Z$ 回転に帰着するもの。
$n_i n_j$ のような高次の積。

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

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

# 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))

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

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

フェルミオン回路の作成

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

アンザッツの準備

InitializeModes クラスを使用してHartree-Fock状態を準備します。窒素の場合、このプロセスは単純に、最初の num_elec_a 個の量子ビットにXゲートを適用し、次に num_elec_b 個の量子ビットに適用するというものです。窒素ではどちらも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次を超える項を含む場合でも、量子ビット接続性が限られているにもかかわらず、ハードウェア上でより効率的に実行できるより浅い回路を作成できます。項のグループ化の後、演算子はその重みに基づいてサンプリングされます。各演算子 $h_i$ について、重み $W_{h_i}$ は以下のように定義されます。

W_{h_i} = |c_i| / \lambda

フェルミオンおよびハードウェアネイティブな最適化

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

QDriftTrotterization パスは、内部で重み計算とサンプリングを使用して、サンプリングに使用する回路を生成します
RelabelModes パスは、フェルミオンモードを並べ替えて量子ビット間の接続性を最適化し、ゲートの深さを削減するために使用できる別の最適化パスです。詳細はAPIリファレンスを参照してください

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

F2QLayout: プリセットのパスマネージャーは TrivialF2QLayout パスを適用し、 $n$ 個のフェルミオンビットを $n$ 個の量子ビットに自明にマッピングします。
F2QSynth: フェルミオンベースの回路命令を量子ビットベースの命令にマッピングするトランスパイルパスです。

qdrift = QDriftTrotterization(num_exc, rng=19)

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

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

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

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

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

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

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

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

ステップ4: 後処理を行い、望ましい古典的な形式で結果を返す

SQDへのビット列の使用

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

ハードウェアの例

この例では20量子ビット(10個の空間軌道)を使用します。この選択は、素早く実行できるべきチュートリアルのための便宜的なものであり、この手法における厳格な上限ではありません。

古典的なステップのコストは、量子ビット数によって直接決まるものではありません。SQDはサンプリングされた配置が張る部分空間に射影されたハミルトニアンを対角化するため、古典的なコストを左右するのは、その選択された部分空間の次元です。これはここでは samples_per_batch、num_batches、および回路が実際に生成する異なる配置の数によって決まり、それに加えて射影されたハミルトニアンを適用するために必要な疎な線形代数の計算も影響します。完全なCI空間は軌道と電子の数に応じて組み合わせ論的に増大しますが、選択される部分空間はそのうちの小さく調整可能な一部分であり、その大きさは直接制御できます。その結果、量子ビット数と古典的な難しさはある程度独立に変化させることができます。より広い軌道空間を適度な部分空間にサンプリングする方が、非常に大きな部分空間で対角化されるより小さな系よりも安価になる場合があります。

したがって実際には、実行可能な系のサイズは、望む精度に必要な部分空間の次元と、固有値ソルバーが利用できるメモリおよびコア数に依存します。一般に、より大きな軌道空間では化学的精度に到達するためにより大きな部分空間が必要となり、これが最終的に分散リソースを必要とする動機となります。このステップをスケールアウトする方法についてはqiskit-addon-sqd-hpcを参照してください。固定のカットオフを想定するのではなく、実際的なアプローチは、報告される部分空間の次元と反復にわたるエネルギーの収束を観察し、エネルギーの改善が止まるか利用可能なメモリを使い尽くすまで部分空間のサイズを増やすことです。

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

偽のビット列の刈り込み

ここでは、追加のステップを実行するかどうかを選択できます。回路の実行からすべてのビット列が得られたら、SQDを実行する前に無効なビット列をフィルタリングして除去することも、刈り込みを行わずに進めることもできます。刈り込みを省略することは、一般にハードウェア実行では望ましいとされます。これは、対称性が破れたショットを配置の復元に利用できる状態のまま残しておくことで、それらを有効な配置に修復し、それらのショットを単に破棄するのではなく部分空間を広げることができるためです。

窒素は7個の $\alpha$ 電子と7個の $\beta$ 電子しか持てないため、出力の前半と後半においてそれぞれ1の数が7個より多いか少ないビット列は破棄できます。ビット列が有効かどうかをチェックし、無効であれば破棄する関数を定義します。偽のビット列をフィルタリングして除去した後、残りは対角化スキームに送られます。以下の PRUNE フラグを使用して、この2つの動作を切り替えます。

刈り込みは、回路の数、発展時間の集合、対角項のフィルタリングと並んで、最終的な部分空間を形作るいくつかの選択肢のうちの一つに過ぎないことに留意してください。刈り込みを行った実行と行わなかった実行を比較することは、他のすべてが固定されている場合にのみ有意義です。C++版の関連資料では、復元ではなく事後選択を行い、他のパラメータも異なることから、この点についてより詳しく議論されています。

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

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

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

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

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

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

次のステップ

推奨事項

この研究に興味を持たれた方は、以下の資料にも興味があるかもしれません。

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

学習成果​

前提条件​

背景​

要件​

セットアップ​

シミュレーター例​

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

ステップ2: 量子ハードウェア実行のための問題の最適化​

ステップ3: Qiskitプリミティブを使用した実行​

ステップ4: 後処理を行い、望ましい古典的な形式で結果を返す​

ハードウェアの例​

次のステップ​