サンプルベース量子対角化による化学ハミルトニアンの解析

使用量の目安: Heron r2 プロセッサで1分未満（注: これはあくまで目安です。実際の実行時間は異なる場合があります。）

学習成果

• 量子処理ユニット（QPU）からサンプリングしたビット列を使用して、SQD Qiskit アドオンで分子系の基底状態エネルギーを近似する方法。

• ffsimを使用して量子化学シミュレーション用の局所ユニタリクラスターヤストロー（LUCJ）回路を構築する方法。

前提条件

• 量子化学と第二量子化

• 量子回路からサンプリングするためのSampler プリミティブの使用方法

背景

このチュートリアルでは、SQD Qiskit アドオンを使用してサンプルベース量子対角化（SQD）アルゴリズムを実装し、ノイズのある量子サンプルを後処理することで平衡結合長における窒素分子 $\text{N}_2$ の基底状態の近似を求める方法を示します。ソフトウェアの詳細は対応するドキュメントに記載されており、簡単な例から始めることもできます。

窒素化学は、肥料用アンモニアの生産を通じて持続可能性と関わっており、この例もその文脈に位置付けられます。触媒が窒素をどのように活性化するかを理解することは、よりエネルギー消費の少ない生産方法の研究の一部です。孤立した窒素分子の基底状態を計算することで、こうしたより広い課題に関連する電子構造法を紹介します。

このチュートリアルは、量子化学（特に分子の基底状態エネルギーの探索）に精通したユーザーを対象としています。ワークフローの詳細なウォークスルーについては、量子対角化アルゴリズムコースを参照してください。

SQDは、量子システムのハミルトニアンなどの量子演算子の固有値と固有ベクトルを、量子コンピューティングと分散古典コンピューティングを組み合わせて求める手法です。古典分散コンピューティングは、量子プロセッサから得られたサンプルを処理し、それらが張る部分空間にターゲットハミルトニアンを射影して対角化するために使用されます。SQDベースのワークフローは以下のステップで構成されます：

1. 回路アンザッツを選択し、参照状態（この場合は Hartree-Fock 状態）に対して量子コンピュータ上で適用します。

2. 得られた量子状態からビット列をサンプリングします。

3. ビット列に対して自己無撞着配置回復手順を実行し、基底状態の近似を取得します。

SQDは、ターゲット固有状態がスパースである場合にうまく機能することが知られています。つまり、波動関数が計算基底状態の集合 $\mathcal{S} = \{|x\rangle \}$ で支持されており、そのサイズが問題のサイズに対して指数的に増大しない場合です。

量子化学

分子系のハミルトニアンは次のように記述できます。

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},$$

ここで、$h_{pr}$ と $h_{prqs}$ は分子積分と呼ばれる複素数であり、コンピュータプログラムを使用して分子の仕様から計算することができます。このチュートリアルでは、PySCF ソフトウェアパッケージを使用して積分を計算します。

分子ハミルトニアンの導出の詳細については、量子化学の教科書（例えば、Szabo と Ostlund 著 Modern Quantum Chemistry）を参照してください。量子化学問題が量子コンピュータにどのようにマッピングされるかの概要については、Qiskit Global Summer School 2024 の講義 Mapping Problems to Qubits を参照してください。

局所ユニタリクラスターヤストロー（LUCJ）アンザッツ

SQDではサンプルを抽出するための量子回路アンザッツが必要です。このチュートリアルでは、物理的な動機付けとハードウェア親和性の両方を兼ね備えた局所ユニタリクラスターヤストロー（LUCJ）アンザッツを使用します。アンザッツ回路の構築には ffsim を使用します。

LUCJアンザッツは、制限されたQubit接続性を持つQPUに適応します。スピン軌道は、アンザッツがSWAPゲートによるルーティングを必要としないようにQubitsにマッピングされます。IBM® ハードウェアは、正方格子またはヘビーヘックス格子のいずれかの量子ビットトポロジーを持っています。ヘビーヘックス格子の場合、以下に示す「ジグザグ」パターンを採用できます。このパターンでは、同じスピンを持つ軌道は線形トポロジーの量子ビットにマッピングされ（赤と青の円）、異なるスピンの軌道間の接続は4番目の空間軌道ごとに存在し、補助量子ビット（紫の円）を介して接続が行われます。

自己無撞着配置回復

自己無撞着配置回復手順は、ノイズのある量子サンプルからできる限り多くの信号を抽出するために設計されています。分子ハミルトニアンは粒子数とスピンZを保存するため、これらの対称性も保存する回路アンザッツを選択することが合理的です。Hartree-Fock状態に適用した場合、ノイズがない設定では、得られる状態は固定された粒子数とスピンZを持ちます。したがって、この状態からサンプリングされた任意のビット列のスピン-$\alpha$ 半分とスピン-$\beta$ 半分は、Hartree-Fock状態と同じハミング重みを持つはずです。現在の量子プロセッサにはノイズが存在するため、測定されたビット列の一部はこの性質に違反します。単純なポストセレクションではこれらのビット列を破棄しますが、それらのビット列にはまだ何らかの信号が含まれている可能性があるため、これは無駄です。自己無撞着回復手順は、後処理でその信号の一部を回復しようとするものです。この手順は反復的であり、入力として基底状態における各軌道の平均占有数の推定値が必要で、まず生のサンプルから計算されます。手順はループ内で実行され、各反復は以下のステップで構成されます：

1. 指定された対称性に違反する各ビット列について、現在の平均軌道占有数の推定値にビット列を近づけるように設計された確率的手順でビットを反転させ、新しいビット列を取得します。

2. 対称性を満たすすべての古いビット列と新しいビット列を収集し、事前に選択した固定サイズの部分集合をサブサンプリングします。

3. 各ビット列の部分集合について、対応する基底ベクトルが張る部分空間にハミルトニアンを射影し（これらの基底ベクトルの説明は前のセクションを参照）、古典コンピュータ上で射影ハミルトニアンの基底状態推定を計算します。

4. 最もエネルギーが低い基底状態推定で、平均軌道占有数の推定値を更新します。

SQDワークフロー図

SQDワークフローは以下の図に示されています：

要件

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

• Qiskit SDK v1.0 以降、visualization サポート付き

• Qiskit Runtime v0.22 以降（pip install qiskit-ibm-runtime）

• SQD Qiskit アドオン v0.12 以降（pip install qiskit-addon-sqd）

• ffsim v0.0.75 以降（pip install ffsim）

セットアップ

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

_needed = {"ffsim": "ffsim", "matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "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")

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

小規模シミュレーターの例

このチュートリアルでは、平衡結合距離付近の窒素分子の基底状態の近似を求めます。まず、実験をシミュレートして動作確認ができるよう、小規模なSTO-6G基底セットを使用します。

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

まず、分子とその性質を指定します。

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="sto-6g",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

converged SCF energy = -108.464957764796
CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000
norb = 8
nelec = (5, 5)

LUCJアンザッツ回路を構築する前に、まず以下のコードセルでCCSD計算を実行します。この計算から得られる $t_1$ および $t_2$ 振幅 は、アンザッツのパラメータの初期化に使用されます。

# 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]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052352

次に、ffsim を使用してアンザッツ回路を作成します。この分子は閉殻Hartree-Fock状態を持つため、UCJアンザッツのスピンバランス型変種である UCJOpSpinBalanced を使用します。from_t_amplitudes メソッドで optimize=True を設定して、$t_2$ 振幅の「圧縮」二重因子分解を有効にします（詳細はffismのドキュメントの The local unitary cluster Jastrow (LUCJ) ansatz を参照してください）。

LUCJアンザッツはQPUの利用可能な接続性に適応するため、アンザッツを作成する前にQPUバックエンドを初期化する必要があります。ここでは、正方格子カップリングマップとLUCJアンザッツが自然に分解されるゲートセットを持つ汎用バックエンドを作成します。次に、ffsim.qiskit.generate_lucj_pass_manager を使用して、「ラダー」レイアウト（詳細はffsimのドキュメントの説明を参照）に従って、LUCJアンザッツを指定バックエンドにトランスパイルするための特化パスマネージャを作成します。この関数は、選択したレイアウトに関連するエラーを最小化するスコアリングヒューリスティックを使用しており、バックエンドが実際のQPUやノイズモデルを持つシミュレーターの場合に重要です。パスマネージャを返すほか、ハードウェア上で実装できるアルファ・ベータカップリングペアも返します。すべてのペアを実装できない場合は警告が出力されます。

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 4
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
coupling_map = CouplingMap.from_grid(norb, 2)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, 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(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

ステップ 2: 量子ハードウェア実行のために最適化する

次に、ターゲットハードウェアに対して回路を最適化します。通常、このステップではハードウェアバックエンドとそのバックエンド用のパスマネージャを初期化しますが、LUCJアンザッツはハードウェアの接続性に適応しているため、これらの処理はすでに前のステップで行いました。残るのは、パスマネージャを回路に対して実行して、QPU上で直接実行できるISA回路にトランスパイルすることだけです。

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})

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

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

rng = np.random.default_rng(12345)
sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

primitive_result = job.result()
pub_result = primitive_result[0]

ステップ 4: 後処理を行い、結果を所望の古典フォーマットで返す

QPU出力の品質を判断する有用な指標は、返された有効な配置の数です。有効な配置とは、正しい粒子数とスピンZを持つもので、ビット列の右半分のハミング重みがスピンアップ電子数に等しく、左半分のハミング重みがスピンダウン電子数に等しいことを意味します。次のセルでは、サンプリングされた配置のうち有効なものの割合を計算します。

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
    )

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")

Fraction of sampled configurations that are valid: 1.0

すべてのビット列が有効なのは、ノイズのないシミュレーター上で回路をサンプリングしているためです。ノイズのあるQPU上で実行する場合、この割合は1未満になりますが、以下のセルで計算されるようにビット列が一様ランダムにサンプリングされた場合に期待される割合よりも大きくなることが望まれます。

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

次に、diagonalize_fermionic_hamiltonian 関数を使用してハミルトニアンの基底状態エネルギーを推定します。この関数は、自己無撞着配置回復（self-consistent configuration recovery）手続きを実行し、ノイズを含む量子サンプルを反復的に改善してエネルギー推定値の精度を向上させます。中間結果を後の解析のために保存できるよう、コールバック関数を渡します。diagonalize_fermionic_hamiltonian の引数の説明については、API ドキュメント を参照してください。

ここでは、diagonalize_fermionic_hamiltonian の initial_occupancies 引数を使用して、基底状態における軌道占有数の初期推定値としてHartree-Fock配置を指定します。このアプローチは、基底状態がHartree-Fock配置に有意なサポートを持つ系には適切ですが、他の状況では適切でない場合もあります。そのような場合は、より高度な計算手法によってより良い初期推定値が得られる場合があります。initial_occupancies を指定することで、有効な配置がサンプリングされなかった場合でも配置回復を実行できます（大規模な回路をノイズのあるQPUでサンプリングする際にこのようなことが起こり得ます）。この引数がない場合、有効な配置が提供されなければ配置回復は失敗してエラーを発生させます。

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 50
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
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)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

Iteration 1
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Iteration 2
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Final energy: -108.59589540253083
Final energy error: 9.19484552071026e-05

結果の可視化

最初のプロットは、このシミュレーションでは最初の反復後に既に 1 mH 以内の精度に達していることを示しています（化学精度は一般に 1 kcal/mol $\approx$ 1.6 mH とされています）。ただし、これは小規模な系であり、サンプルがノイズレスであるため配置回復は必要ありません。より大規模な系をノイズのあるQPU上で実行した場合は、複数回の配置回復反復が必要になることがあり、最終的な精度が悪化する場合があります。一般的に、配置回復の反復回数を増やすか、バッチあたりのサンプル数を増やすことで、エネルギーの精度をさらに向上させることができます。

2番目のプロットは、最終反復後の各空間軌道の平均占有数を示しています。スピンアップ電子とスピンダウン電子の両方が、解において高い確率で最初の5つの軌道を占有していることがわかります。

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

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_yscale("log")
axs[0].set_ylim(1e-5, 1e-2)
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})

plt.tight_layout()
plt.show()

大規模ハードウェアの例

次に、実際の量子ハードウェア上でより大規模な例を実行します。ここでは、cc-pVDZ基底関数系から窒素分子のアクティブスペースを導出します。また、以下の変更も行います：

• ゲートエラーによるノイズを減らすため、n_reps を 1 に下げる
• バックエンドのプロセッサファミリーに基づいてアンザッツの接続性を設定する（Heronの場合はヘビーヘックス、Nighthawkの場合は正方格子）
• より大きな部分空間次元を可能にするため、max_dim を 1000 に増やす。

ステップ 1〜4

ここでは、すべてのステップを1つのワークフローにまとめてより大規模なスケールで実行します。これを実際の量子ハードウェア上で実行します。

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="cc-pvdz",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# 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]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(f"Using backend {backend.name}")

# Create pass manager
connectivity_map = {"Heron": "heavy-hex", "Nighthawk": "square"}
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity_map[
        backend.configuration().processor_type["family"]
    ],
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, 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(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 1000
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the
# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

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_yscale("log")
axs[0].set_ylim(1e-3, 1e-1)
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})

plt.tight_layout()
plt.show()

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450041
Using backend ibm_boston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.01677
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.17922180754668
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.17713541195594
		Subspace dimension: 962361
	Subsample 2
		Energy: -109.17407767345108
		Subspace dimension: 968256
Iteration 2
	Subsample 0
		Energy: -109.1873108379026
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.1876970140745
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.18255107152576
		Subspace dimension: 1000000
Iteration 3
	Subsample 0
		Energy: -109.19052148187927
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19316375313544
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19087348337145
		Subspace dimension: 1000000
Iteration 4
	Subsample 0
		Energy: -109.19374847699348
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19517039328333
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19422416032815
		Subspace dimension: 1000000
Iteration 5
	Subsample 0
		Energy: -109.19580768221108
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19608553581007
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19548131833729
		Subspace dimension: 1000000
Final energy: -109.19608553581007
Final energy error: 0.03194368084709254

次のステップ

推奨事項

この研究が興味深いと感じた方は、以下の資料も参照してください：

• フェルミオン格子モデルのサンプルベース量子Krylov対角化 - 変分アンザッツの代わりに時間発展回路を使用した関連チュートリアル
• Diceソルバーを使用したSQD化学ワークフローのスケールアップ - より効率的なDiceソフトウェアを対角化に使用する方法を示すページ
• SQDアドオンAPIドキュメント - diagonalize_fermionic_hamiltonian 関数のリファレンス
• Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer - このチュートリアルの基となった論文