{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "048b37e0",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"SqDRIFT 基底状態推定のためのアルゴリズム\"\n",
        "description: \"SqDRIFT 基底状態推定の問題に対して、 qDRIFT とSKQDという2つの著名なアルゴリズムを組み合わせると同時に、回路の深さを低減する。\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore cisolver ECORE combinatorially multiset */}\n",
        "\n",
        "<span id=\"sqdrift-algorithm-for-ground-state-estimation\" />\n",
        "\n",
        "# SqDRIFT 基底状態推定のためのアルゴリズム\n",
        "\n",
        "推定*所要時間：Heron r3 プロセッサで 180 秒（注：これはあくまで推定値です。 （実行時間は状況によって異なる場合があります。）*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "sqdrift-cpp-note",
      "metadata": {},
      "source": [
        "<Admonition type=\"note\" title=\"C++版をお探しですか？\">\n",
        "  このチュートリアルでは、 Python を使用しています。 C++による実装（ソースコードおよびビルド手順を含む）については、「 [C++ SqDRIFT チュートリアル](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/cpp) 」を参照してください。\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5edbd059",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## 学習成果\n",
        "\n",
        "* トロッター化に比べて、奥行きが浅い回路を作成する方法について学びましょう\n",
        "* qDRIFT とSQDを用いた基底状態推定のエンドツーエンドのワークフローを順を追って解説します\n",
        "* 他のQiskitアドオンと組み合わせて `qiskit-fermions` 、このようなワークフローを実装する方法について学びましょう\n",
        "\n",
        "このチュートリアルは、教育目的で Python ノートブックとして提供されています。\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## 前提条件\n",
        "\n",
        "* 「 [サンプルベース量子対角化（SQD）](/docs/addons/qiskit-addon-sqd) 」の概要を読む\n",
        "* 「 [サンプルベースのクリロフ量子対角化（SKQD）](/learning/courses/quantum-diagonalization-algorithms/skqd) 」のレッスンをご覧ください\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## 背景\n",
        "\n",
        "[SqDRIFT](https://arxiv.org/abs/2508.02578) これはSKQDの変種であり、ビット列をサンプリングするためのアンザッツを選択する必要性を、対象のハミルトニアンから直接構築された時間発展回路の集合体に置き換えたものである。 これは、ハミルトニアンの係数に基づいて、より小さな時間発展演算子をハミルトニアンからサブサンプリングすることで実現され、これは「 qDRIFT トロッター化法」として知られている。\n",
        "\n",
        "このチュートリアルでは、 [Qiskit Fermions](/docs/addons/qiskit-fermions) を利用して、 [qDRIFT](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.070503) アルゴリズム向けのより自然なフェルミオン回路を作成します。その後、フェルミオン向けのレイアウトおよび合成処理を適用し、最終的に回路を従来のQiskitパイプラインに組み込んでハードウェア上で実行します。\n",
        "\n",
        "ハミルトニアンを次のような形とする：\n",
        "\n",
        "$$\n",
        "H = \\sum_{i=1}^{N} c_i h_i\n",
        "$$\n",
        "\n",
        "ここで、一般性を失うことなく、 $c_i > 0$ を満たし、かつ $h_i$ の最大固有値の絶対値が $1$ に等しいものと仮定する。符号付きまたは複素数の前因子はいずれも $h_i$ に吸収されるため、係数 $c_i$ は厳密に正の重みとなり、 $h_i$ は各項の方向を表す。 ここで、 $N$ はハミルトニアンの項の数（あるいは、グループ分け後のグループの数）であり、これはハミルトニアンの性質である。これは、以下で $n$ と表記される、単一の回路にサンプリングされる演算子の数とは異なる。\n",
        "\n",
        "qDRIFT アルゴリズムは、目標時間 $t$ に対して、ある演算子 $V_k$ を実現する。ここで、 $k$ は $1 \\cdots K$ から始まり、次のように定義される $k_{th}$SqDRIFT 回路を表す：\n",
        "\n",
        "$$\n",
        "V_k = \\prod_{j=1}^{n} e^{-i h_{k_j} \\lambda t / n }\n",
        "$$\n",
        "\n",
        "ここで、 $n$ は回路あたりのサンプリングされた演算子の数、 $K$ はアンサンブル内の回路の数である。 この製品は、すべての $N$ ハミルトニアン項ではなく、 $n$ の抽出結果に基づいて動作します。また、項は重複ありで抽出されるため、単一の $V_k$ 内に同じ $h_i$ が複数回出現する可能性があります。\n",
        "\n",
        "数量：\n",
        "\n",
        "$$\n",
        "\\lambda = \\sum_{i=1}^{N} c_i\n",
        "$$\n",
        "\n",
        "は係数の $L_1$ ノルムであるため、 $n$ の各ステップは、どの項が抽出されたかに関係なく、同じ時間 $\\lambda t / n$ だけ進行する。 ステップ角の均一性は、 qDRIFT: の特徴であり、係数は、その項がどれだけ回転するかではなく、 *その項がどれだけ頻繁に*引き出されるかによって結果に影響を与えます。 各指標は、以下の分布からサンプリングされます：\n",
        "\n",
        "$$\n",
        "P[k_i] = \\frac{c_i}{\\lambda}\n",
        "$$\n",
        "\n",
        "したがって、 $(k_1, \\ldots, k_n)$ という数列は、この分布から抽出された項の添字のランダムな列である。 $c_i$ は正の値であり、その和は $\\lambda$ となるため、これは正規化された確率分布であり、ランダムな抽出に基づく結果のチャネルの期待値は、 $H$ の下での進化を近似する。この誤差は、 $n$ が大きくなるにつれて減少する。 なお、近似誤差は、項の数 $N$ ではなく、 $\\lambda$ に依存することに注意してください。\n",
        "\n",
        "（『 SqDRIFT 』の論文では、項の数を $\\mathcal{N}$、数列の長さを $N$ と表記しているが、ここでは両者を明確に区別するために、 $N$ および $n$ という表記を用いる。）\n",
        "\n",
        "このチュートリアルでは、このようなランダム化された回路のアンサンブルを生成する方法について説明します。 これらの回路を作成した後、異なる演算子に対してクリロフ部分空間を作成する場合と同様に、異なる時間パラメータを持つ複数のそのような演算子からビット列をサンプリングする。 これにより、基底状態ベクトルとサンプリングされたビット列との間の重なりがより高くなる。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cda4ef2a",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## 要件\n",
        "\n",
        "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください\n",
        "\n",
        "* Python （ 3.10 以上）の仮想環境\n",
        "* pip>=25.1\n",
        "* qiskit ≈ 2.5\n",
        "* qiskit-fermions==0.1.0 （この名称は複数形であることに注意してください）\n",
        "* numpy\n",
        "* pyscf\n",
        "* qiskit-aer\n",
        "* qiskit-ibm-runtime\n",
        "* qiskit-addon-sqd\n",
        "\n",
        "以下のコマンドで、必要なパッケージをすべてインストールできます：\n",
        "\n",
        "```\n",
        "pip install \"qiskit~=2.5\" \"qiskit-fermions==0.1.0\" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e24facd1",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## セットアップ\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "1e5d4294",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Third-party scientific computing\n",
        "import numpy as np\n",
        "\n",
        "# PySCF\n",
        "from pyscf import tools, ao2mo, fci\n",
        "\n",
        "# Qiskit core\n",
        "from qiskit import transpile\n",
        "from qiskit.primitives import BitArray\n",
        "\n",
        "# Qiskit Aer\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "# IBM Quantum Compute Service\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler\n",
        "\n",
        "# Qiskit Fermions\n",
        "from qiskit_fermions.operators.library import FCIDump\n",
        "from qiskit_fermions.operators import FermionOperator\n",
        "from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms\n",
        "from qiskit_fermions.operators.terms.grouping import (\n",
        "    group_terms_by_electronic_structure,\n",
        ")\n",
        "from qiskit_fermions.operators.terms.ordering import canonical_order\n",
        "from qiskit_fermions.circuit import FermionicCircuit\n",
        "from qiskit_fermions.circuit.library import Evolution\n",
        "from qiskit_fermions.transpiler import FermionicPassManager\n",
        "from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager\n",
        "from qiskit_fermions.transpiler.passes import QDriftTrotterization\n",
        "from qiskit_fermions.circuit.library import InitializeModes\n",
        "\n",
        "# Qiskit addon SQD\n",
        "from qiskit_addon_sqd.fermion import (\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        "    SCIResult,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "65478ff1",
      "metadata": {},
      "source": [
        "<span id=\"simulator-example\" />\n",
        "\n",
        "## シミュレータの例\n",
        "\n",
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### ステップ1：古典的な入力を量子問題に写像する\n",
        "\n",
        "**FCIDump の読み込みと準備**\n",
        "\n",
        "このチュートリアルでは、窒素の電子構造ハミルトニアンを読み込みます（ N2 ）。 フェルミオン演算子を作る方法は他にもあります。  のドキュメントを参照してください [`qiskit_fermions.operators.library`](https://qiskit.github.io/qiskit-fermions/stable/0.1/pydoc/qiskit_fermions.operators.library.html#module-qiskit_fermions.operators.library)。\n",
        "\n",
        "**このFCIDumpについて。** このファイルは、最小 STO-3G 基底関数系における窒素分子（ $N_2$ ）について記述 `N2_sto_3g` しており、原子間距離は 1.09$\\AA$ に設定されている。これは、実験的に得られた平衡結合長である。 そのヘッダーには `NORB=10`、 `NELEC=14`、、およびが宣言されている `MS2=0`。すなわち、10個の空間軌道（したがって、20個のスピン軌道、およびジョーダン・ウィグナー表現の下では20個の量子ビット）、スピンシングレット状態にある14個の電子、つまり7個の $\\alpha$ 電子と7個の $\\beta$ 電子である。 すべての軌道には対称性ラベル「1」が割り当てられており、つまり、点群対称性は利用されていない。 これは全空間の STO-3G ダンプであるため、軌道は凍結されておらず、相関空間も十分に小さいため、次のセルに示すように、比較のために古典的に正確なFCI基準エネルギーを計算することができます。\n",
        "\n",
        "PySCF: を実行することで、同等のファイルを再生成できます\n",
        "\n",
        "```python\n",
        "from pyscf import gto, scf, tools\n",
        "\n",
        "mol = gto.M(atom=\"N 0 0 0; N 0 0 1.09\", basis=\"sto-3g\", symmetry=False)\n",
        "mf = scf.RHF(mol).run()\n",
        "tools.fcidump.from_scf(mf, \"N2_sto_3g\")\n",
        "```\n",
        "\n",
        "積分は収束したSCF軌道に依存するため、再生成されたファイルは、軌道の位相や順序において出荷時のファイルと異なる場合がありますが、総エネルギーには影響しません。\n",
        "\n",
        "**ファイルの取得。** この [GitHub リポジトリ](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/fcidump_files)で FCIDump を見つけてください。 以下のセルを実行すると、チュートリアルの残りの部分で想定されている場所にそのデータを取得できます。\n",
        "\n",
        "まず、pyscf が提供する `cisolver` を使用して、基準エネルギーを求めます。 これが、我々が扱っている分子の真の基底状態エネルギーである。 このため、まず `nelec`、それぞれ軌道数と電子数を表す `norb` と を定義する。 次に、それぞれ1電子積分および2電子積分を表す および `h1e` を定義する `h2e`。 これらはすべて、後でSQDにも使用されることになる。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "f15f2d83",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g\n"
          ]
        }
      ],
      "source": [
        "import os\n",
        "from urllib.request import urlopen\n",
        "\n",
        "# The FCIDump is stored with this tutorial in the Qiskit documentation repository.\n",
        "FCIDUMP_URL = \"https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "FCIDUMP_PATH = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "if not os.path.exists(FCIDUMP_PATH):\n",
        "    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)\n",
        "    with urlopen(FCIDUMP_URL) as response:\n",
        "        contents = response.read()\n",
        "    with open(FCIDUMP_PATH, \"wb\") as f:\n",
        "        f.write(contents)\n",
        "    print(f\"Downloaded FCIDump to {FCIDUMP_PATH}\")\n",
        "else:\n",
        "    print(f\"Using existing FCIDump at {FCIDUMP_PATH}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "41b178cf",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Parsing assets/sqdrift/fcidump_files/N2_sto_3g\n",
            "Reference FCI Energy  = -107.6481842917 Ha\n",
            "Nuclear Repulsion Energy = 23.7887003074 Ha\n"
          ]
        }
      ],
      "source": [
        "name = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "fcidump = tools.fcidump.read(name)\n",
        "\n",
        "# Extract metadata from the FCIDump header\n",
        "norb = fcidump[\"NORB\"]  # number of spatial orbitals\n",
        "nelec = fcidump[\"NELEC\"]  # total number of electrons\n",
        "e_nuc = fcidump[\"ECORE\"]  # nuclear repulsion / core energy\n",
        "ms2 = fcidump[\"MS2\"]  # 2S (spin)\n",
        "\n",
        "num_elec_a = (nelec + ms2) // 2  # alpha electrons\n",
        "num_elec_b = (nelec - ms2) // 2  # beta  electrons\n",
        "\n",
        "# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)\n",
        "h1e = fcidump[\"H1\"]  # shape (norb, norb)\n",
        "h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)\n",
        "    1, fcidump[\"H2\"], norb\n",
        ")\n",
        "\n",
        "cisolver = fci.direct_spin1.FCI()\n",
        "cisolver.max_cycle = 200\n",
        "cisolver.conv_tol = 1e-12\n",
        "\n",
        "e_fci, _ = cisolver.kernel(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    norb,\n",
        "    (num_elec_a, num_elec_b),\n",
        "    ecore=e_nuc,  # adds nuclear repulsion to the final energy\n",
        ")\n",
        "\n",
        "reference_energy = e_fci\n",
        "\n",
        "print(f\"Reference FCI Energy  = {reference_energy:.10f} Ha\")\n",
        "\n",
        "nuclear_repulsion_energy = fcidump[\"ECORE\"]\n",
        "print(f\"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7818e18a",
      "metadata": {},
      "source": [
        "**ハミルトニアンの読み込み**\n",
        "\n",
        "必要なデータが準備できたので、FCIファイルから、以下の形式と互換性のあるハミルトニアンを読み込みます。 `qiskit-fermions`\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "565b92fc",
      "metadata": {},
      "outputs": [],
      "source": [
        "fcidump = FCIDump.from_file(name)\n",
        "hamiltonian = FermionOperator.from_fcidump(fcidump)\n",
        "num_modes = 2 * fcidump.norb"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d88bb917",
      "metadata": {},
      "source": [
        "**フェルミオンを用いたワークフロー `qiskit-fermions`**\n",
        "\n",
        "まず、トランスパイラーのパスやフェルミオン回路特有のゲートを提供する を用いて `qiskit-fermions`、ハミルトニアンをフェルミオン回路モデルにマッピングします。 これらは、このワークフローにおいて、Qiskitの従来のトランスパイラーによる処理が行われる前に使用されます。\n",
        "\n",
        "**項目のグループ化**\n",
        "\n",
        "結果の再現性を確保するため、まず を用いて、項をその構造のみに基づいて並べ `canonical_order` 替えます。 したがって、リスト `canon` 内の演算子の順序は固定されています。 これにより、作成された演算子の再現性が確保されます。というのも、今後使用する` `QDriftTrotterization` pass`が、ランダムなインデックスをサンプリングして qDRIFT 演算子を作成するからです。\n",
        "\n",
        "このステップでは、電子構造ハミルトニアンに存在する多くの対称性を活用し、係数が同じである関連する項をグループ化します。 そうすることで、 qDRIFT プロトコルがサンプリングを行う演算子係数の分布は変化しますが、その収束保証には影響しません。 重要なことに、対称性によって関連する項をグループ化することで、パウリ項が好都合に相殺され、それらの作用下で状態を時間発展させる際の回路深度が全体的に短くなる。\n",
        "\n",
        "`qiskit-fermions` このグループ化を自動的に行ってくれる関 `group_terms_by_electronic_structure` 数が用意されています。\n",
        "\n",
        "なお、この式では、項が[通常の順序で](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.operators.FermionOperator.html#qiskit_fermions.operators.FermionOperator.normal_ordered)並んでいることを前提 `group_terms_by_electronic_structure` としていることに注意してください。\n",
        "\n",
        "**対角項のフィルタリング**\n",
        "\n",
        "回路を生成するために使用するハミルトニアンから対角成分を除去することで、 $n$ qDRIFT サンプリングスロットが、構成間の状態分布を変化させる項に割り当てられるようにします。 こうした項は、次のステップでゲート `Evolution` が構築される前のこの段階で、ハミルトニアンから除外しておくのが最善である。\n",
        "\n",
        "ここで問題となる項は、職業-数基底において対角にあるもの、すなわち、数演算子の積である $a^\\dagger_i a_i$ である。この定義に該当する項には、以下の3種類がある：\n",
        "\n",
        "* **定数エネルギーオフセット**。これはゼロ数演算子の積であり、その時間発展はグローバルな位相のみに寄与する；\n",
        "* **個々の数演算子**$n_i$。これらは、時間発展により単一量子ビットの $Z$ 回転に還元される；\n",
        "* $n_i n_j$ などの**高階**積。\n",
        "\n",
        "これらはいずれも、それ自体では「占領数構成」間の住民の移動を引き起こすことはなく、すでに存在する構成の「フェーズ」にのみ作用する。 しかし、それらは決して無関係なものではない。これらの相対位相は、回路の後半にある励起項によって生じる干渉に影響を与えるため、それらを除去すると、実際に生成される挙動が変化し、サンプリング分布も変化する可能性がある。 これは、回路生成段階における意図的な近似であり、サンプリングを励起項に集中させるために行われるものであって、サンプリングされた分布をそのまま残す段階ではない。 上記の対称性グループ化とは異なり、このフィルタは、 qDRIFT の収束保証を損なうことなく、進化の対象となる演算子を変更します。 したがって、これらの回路はもはや完全ハミルトニアン下での進化を近似するものではなく、 qDRIFT の誤差上限は、元の演算子ではなく、フィルタリングされた演算子に適用されることになる。 ここではこれが許容されるのは、回路が構成案を提案するために用いられる単なるサンプリング手法に過ぎないためである。フィルタは回路の構築に使用されるハミルトニアンのみに適用されるのに対し、その後の古典的な対角化では対角項を含む完全なハミルトニアンが使用されるため、エネルギー推定そのものから項が失われることはない。 SQDの精度は、その古典的なステップに依存しており、このステップは、構成がどのように提案されたかに関わらず、サンプリングされた部分空間において変分的なままである。\n",
        "\n",
        "この関 `filter_diagonal_terms()` 数は、演算子からそのような項をその場で削除します。 これは、それらの「正規順序構造」――すなわち、生成モードの多集合が消滅モードの多集合と一致すること――によってそれらを特定するため、すでに正規順序付けられている演算子に対してのみ有効である。 この仮定は実行時に検証されません。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "464a4470",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "5060\n"
          ]
        }
      ],
      "source": [
        "# Apply automatic grouping\n",
        "canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))\n",
        "exit_code = group_terms_by_electronic_structure(\n",
        "    canon, num_modes, two_body_physicist_order=False\n",
        ")\n",
        "filter_diagonal_terms(canon)\n",
        "\n",
        "print(len(canon.groups))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fe33ada4",
      "metadata": {},
      "source": [
        "ハミルトニアンの項をグループ分けしたので、回路のアンサンブルを生成するために、以下のパラメータを決定します：\n",
        "\n",
        "* 生成する回路の数： `num_circuits`\n",
        "* 励起群に基づく各回路の長さ： `num_exc`\n",
        "* 進化時間が異なる要因： `times`\n",
        "\n",
        "**フェルミオン回路の構築**\n",
        "\n",
        "それでは、各時間ステップごとにフェルミオン回路を作成していきます。 各回路は、先ほど定義した進化時間を用いた単一の進化ゲートで構成されます。 進化演算子はハミルトニアンである。 その後、これらの回路に対してトランスパイラー処理を行い、 qDRIFT 回路を生成します。\n",
        "\n",
        "**アンツァッツの準備**\n",
        "\n",
        "class `InitializeModes` を使用して、ハートリー・フォック状態を生成します。 窒素の場合、この処理は、まず最初の 量子 `num_elec_a` ビットに X ゲートを適用し、次に 量子 `num_elec_b` ビットに X ゲートを適用するという単純なもので、窒素の場合、これらはいずれも 7 個ずつです。 この状態は、窒素の7つの $\\alpha$ 電子と7つの $\\beta$ 電子を表しています。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "140abcc6",
      "metadata": {},
      "outputs": [],
      "source": [
        "# SqDRIFT parameters\n",
        "times = [1.0, 10.0]  # Total evolution times used for the subspace creation\n",
        "num_exc = 10  # Number of excitation groups per circuit\n",
        "num_circuits = 200  # Number of circuits to generate\n",
        "\n",
        "\n",
        "init_circuits = []\n",
        "\n",
        "hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))\n",
        "\n",
        "for time in times:\n",
        "    evo_gate = Evolution(num_modes, canon, time)\n",
        "    circ = FermionicCircuit(num_modes)\n",
        "    circ.append(hf_gate, circ.modes)\n",
        "    circ.append(evo_gate, circ.modes)\n",
        "    init_circuits.append(circ)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c18ed9a5",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### ステップ2：量子ハードウェアでの実行に向けて問題を最適化する\n",
        "\n",
        "回路が完成したので、まずは で利用可能なパスを使用してフェルミオンレベルの最適化を行い `qiskit-fermions` 、その後、選択したバックエンド向けに回路をトランスパイルします。 これはシミュレータ実験であるため、まずは AerSimulator についてこの作業を行います。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7de434b9",
      "metadata": {},
      "source": [
        "**各グループごとの重量の算出**\n",
        "\n",
        "このステップでは、ハミルトニアンの係数に比例する確率に従って、項の qDRIFT サンプリングを確率的に実行します。 qDRIFT トランスパイラ・パスが、この処理を代行してくれます。 これにより、量子ビット間の接続性が限られている場合でも、ハミルトニアンに長距離結合や二次以上の項が含まれている場合でも、ハードウェア上でより効率的に実行できる、より浅い回路を作成できるようになりました。\n",
        "項のグループ化を行った後、重みに基づいて演算子をサンプリングします。 各演算子 $h_i$ について、重み $W_{h_i}$ は次のように定義される：\n",
        "\n",
        "$$\n",
        "W_{h_i} = |c_i| / \\lambda\n",
        "$$\n",
        "\n",
        "**フェルミオンおよびハードウェアネイティブの最適化**\n",
        "\n",
        "この関数は、を受け取り `FermionicCircuit` 、ハードウェア上で実行できるようトランスパイルできる最適化された最終回路を生成する `MultiStagePassManager` を返します `generate_preset_jw_pass_manager()` 。 そのデフォルトの最適化ステージを、私たちのパ `QDriftTrotterization` スを含む `FermionicPassManager` ものに置き換えます：\n",
        "\n",
        "* この `QDriftTrotterization` パスでは、内部で重み計算とサンプリングを行い、サンプリングに使用する回路を生成します\n",
        "* このパス `RelabelModes` は、フェルミオンモードの順列を変更して量子ビット間の接続性を最適化し、ゲートの深さを削減するために使用できる、もう1つの最適化パスです。詳細については、 [APIリファレンス](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.RelabelModes.html#qiskit_fermions.transpiler.passes.RelabelModes)をご覧ください\n",
        "\n",
        "残りの段階は自動的に実行 `MultiStagePassManager` され、フェルミオンから量子ビットへのマッピングをすべて処理します：\n",
        "\n",
        "* [F2QLayout](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.TrivialF2QLayout.html) ：プリセット・パス・マネージャーは、 $n$ のフェルミオン・ビットを $n$ の量子ビットに単純にマッピングするパスを適用します `TrivialF2QLayout` 。\n",
        "* [F2QSynth](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.F2QSynthesis.html) ：フェルミオンベースの回路命令をキュービットベースの命令に変換するトランスパイル処理。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "34dbca48",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "400\n"
          ]
        }
      ],
      "source": [
        "qdrift = QDriftTrotterization(num_exc, rng=19)\n",
        "\n",
        "pm = generate_preset_jw_pass_manager()\n",
        "pm.optimization = FermionicPassManager([qdrift])\n",
        "\n",
        "sqdrift_circuits = []\n",
        "for circ in init_circuits:\n",
        "    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))\n",
        "\n",
        "for circ in sqdrift_circuits:\n",
        "    circ.measure_all()\n",
        "\n",
        "print(len(sqdrift_circuits))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6f71b1cb",
      "metadata": {},
      "source": [
        "フェルミオンレベルの最適化が完了したので、シミュレータ上で実行するために回路をトランスパイルすることができます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "de130514",
      "metadata": {},
      "outputs": [],
      "source": [
        "simulator = AerSimulator()\n",
        "shots = 100\n",
        "\n",
        "transpiled_circuits = transpile(sqdrift_circuits, simulator)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "98db4d80",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### ステップ 3: `Qiskit primitives` を使用して実行する\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "74496e28",
      "metadata": {},
      "source": [
        "回路が完成したので、 AerSimulator 上で Qiskit primitives を使用して実行することができます。 異なる回路からの集計結果をすべて合算します。 これらをブールベクトルに変換してから、最終的にSQDを用いて後処理を行います。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "24b9df39",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Executing 400 circuits with 100 shots each...\n",
            "400 length before post processing\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    f\"Executing {len(transpiled_circuits)} circuits with {shots} shots each...\"\n",
        ")\n",
        "\n",
        "job = simulator.run(transpiled_circuits, shots=shots)\n",
        "result = job.result()\n",
        "\n",
        "all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]\n",
        "\n",
        "print(len(all_counts), \"length before post processing\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b797c220",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### ステップ4：後処理を行い、結果を所望の従来の形式で出力する\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2129d0ac",
      "metadata": {},
      "source": [
        "**SQDにおけるビット文字列の使用**\n",
        "\n",
        "これで、選択したビット列に対して対角化アルゴリズムを実行し、分子の基底状態のエネルギーに対応する最小の固有値を求めることができます。 コールバック関数を作成し、初期占有率を宣言し、パラメータを設定してから、最終的に対角化スキームを実行します。 コールバック関数は、各反復ごとに、現在の反復番号と現在の固有値の推定値を出力するために使用されます。\n",
        "\n",
        "最後に、基底状態の推定値を求めるために、得られたエネルギー `nuclear_repulsion_energy` に を加えます。\n",
        "\n",
        "**注** ：サブスペースの次元は、ノイズのないシミュレータであっても、反復ごとに一定ではありません。各サブサンプルごとに異なる構成セットが抽出され、復元ステップによって反復ごとにプールの形状が変更されるため、報告される次元はサブサンプルごとに異なります。 ノイズレスサンプリングだけでは、選択された部分空間の次元を決定づけることはできない。 しかし、ハードウェア実行では、ノイズの混入したショットによって粒子数の対称性が破られ、構成の復元によってそれらが追加の基底ベクトルとなるため、体系的により大きな部分空間が得られる傾向がある。 そのため、ハードウェアのセクションでは、ビット文字列を剪定するための別の手順についても紹介します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "f7c5b2ed",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "40000\n",
            "  Alpha electrons: 7\n",
            "  Beta electrons: 7\n",
            "  Number of orbitals: 10\n",
            "  Number of spin orbitals (qubits): 20\n",
            "Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)\n",
            "\n",
            "Running SQD with configuration recovery...\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64767025226178\n",
            "\t\tSubspace dimension: 5538\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64772799119115\n",
            "\t\tSubspace dimension: 5670\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64765512281548\n",
            "\t\tSubspace dimension: 5767\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64795948524682\n",
            "\t\tSubspace dimension: 6080\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64806617355072\n",
            "\t\tSubspace dimension: 6300\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64802260640258\n",
            "\t\tSubspace dimension: 6308\n",
            "FINAL SQD RESULTS\n",
            "Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712\n",
            " 0.99301927 0.0373282  0.0373266  0.00944508]\n",
            "Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737\n",
            " 0.99302145 0.03733536 0.03733399 0.0094585 ]\n",
            "Reference Energy: -107.6481842917 Ha\n",
            "Computed Energy:  -107.6480661736 Ha\n",
            "Error:            1.1811817564e-04 Ha\n"
          ]
        }
      ],
      "source": [
        "combined_counts = {}\n",
        "for counts in all_counts:\n",
        "    for bitstring, count in counts.items():\n",
        "        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count\n",
        "\n",
        "bit_array = BitArray.from_counts(combined_counts)\n",
        "print(bit_array.num_shots)\n",
        "\n",
        "print(f\"  Alpha electrons: {num_elec_a}\")\n",
        "print(f\"  Beta electrons: {num_elec_b}\")\n",
        "print(f\"  Number of orbitals: {norb}\")\n",
        "print(f\"  Number of spin orbitals (qubits): {2*norb}\")\n",
        "print(f\"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}\")\n",
        "\n",
        "# SQD parameters\n",
        "samples_per_batch = 300\n",
        "num_batches = 3\n",
        "max_iterations = 5\n",
        "\n",
        "initial_occupancies = (\n",
        "    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha\n",
        "    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta\n",
        ")\n",
        "\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "# Run SQD with configuration recovery\n",
        "print(\"\\nRunning SQD with configuration recovery...\")\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=1e-3,\n",
        "    occupancies_tol=1e-3,\n",
        "    max_iterations=max_iterations,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    seed=42,\n",
        "    callback=callback,\n",
        ")\n",
        "\n",
        "computed_energy = result.energy + nuclear_repulsion_energy\n",
        "\n",
        "print(\"FINAL SQD RESULTS\")\n",
        "print(f\"Orbital occupancies (alpha): {result.orbital_occupancies[0]}\")\n",
        "print(f\"Orbital occupancies (beta): {result.orbital_occupancies[1]}\")\n",
        "\n",
        "\n",
        "energy_error = abs(computed_energy - reference_energy)\n",
        "print(f\"Reference Energy: {reference_energy:.10f} Ha\")\n",
        "print(f\"Computed Energy:  {computed_energy:.10f} Ha\")\n",
        "print(f\"Error:            {energy_error:.10e} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "84705007",
      "metadata": {},
      "source": [
        "<span id=\"hardware-example\" />\n",
        "\n",
        "## ハードウェアの例\n",
        "\n",
        "この例では、20キュービット（10個の空間軌道）を使用しています。 その選択は、チュートリアルをスムーズに実行するための便宜上の措置であり、この手法に対する絶対的な制限というわけではありません。\n",
        "\n",
        "古典ステップのコストは、量子ビット数によって直接決定されるわけではない。 SQDは、 *サンプリングされた*構成が張る部分空間に射影されたハミルトニアンを対角化するため、古典的な計算コストを左右するのは、その選択された部分空間の次元（ここでは `samples_per_batch`、 `num_batches`、および回路が実際に生成する異なる構成の数によって決まる）と、射影されたハミルトニアンを適用するために必要な疎な線形代数の計算量である。 CI空間全体は、軌道数や電子数に応じて組み合わせ的に増加しますが、選択された部分空間はそのごく一部であり、調整可能な範囲に限定されており、その大きさを直接制御することができます。 したがって、量子ビットの数と古典的な計算難度は、ある程度独立して調整することが可能です。つまり、広範な軌道空間から適度な部分空間をサンプリングする方が、非常に広い空間にわたって対角化される小規模な系よりも計算コストが低くなる場合があります。\n",
        "\n",
        "したがって、実際には、実現可能なシステムの規模は、求める精度に必要な部分空間の次元と、固有値ソルバーが利用できるメモリおよびコア数によって決まります。 通常、軌道空間が大きくなると、化学的精度を達成するためにより大きなサブスペースが必要となります。これが、最終的には分散リソースの導入につながる要因となります。このステップをスケールアウトする方法については、 [qiskit-addon-sqd-hpc を](https://qiskit.github.io/qiskit-addon-sqd-hpc/)参照してください。 固定のカットオフ値を想定するよりも、反復計算の過程で報告される部分空間の次元とエネルギーの収束状況を注視し、エネルギーの改善が止まるか、利用可能なメモリが尽きるまで部分空間のサイズを拡大していくのが実用的なアプローチである。\n",
        "\n",
        "*注：* ハードウェアのノイズによるサンプリング誤差のため、ハードウェア実行時に対角化のために生成される部分空間は、シミュレータを使用した場合に得られるものよりも大きくなります。 これにより、対角化したい部分空間の次元は増えますが、SQDはノイズに対して頑健であるため、このワークフローでも正確な答えが得られます。\n",
        "\n",
        "**不要な文字列の削除**\n",
        "\n",
        "ここでは、追加の手順を実行するかどうかの選択ができます。 回路の実行結果からすべてのビット文字列が得られたら、SQDを実行する前に無効なビット文字列をフィルタリングするか、あるいは剪定を行わずに処理を進めるか、どちらかを選択できます。 ハードウェア実行においては、一般的に剪定を省略する方が望ましい。そうすることで、対称性が破れたショットを構成回復の対象として残すことができ、それらを有効な構成に修復して部分空間を拡大することができるためである。これに対し、剪定を行うと、そうしたショットは即座に破棄されてしまう。\n",
        "\n",
        "窒素は $\\alpha$ と $\\beta$ の電子をそれぞれ7個しか持つことができないため、出力の前半と後半で 1s が7個より多い、あるいは少ないビット列はすべて除外することができます。 ビット文字列が有効かどうかを検証し、有効でない場合はそれらを破棄する関数を定義します。 誤ったビット列を除去した後、残りのビット列は対角化処理に送られます。 以下のフラグ `PRUNE` を使用して、2つの動作を切り替えてください。\n",
        "\n",
        "剪定は、回路数、進化時間のセット、対角項フィルタリングと並んで、最終的な部分空間を形作るいくつかの選択肢のうちの1つにすぎないことを念頭に置いておいてください。 プリニングを施した実行と施していない実行を比較しても、他のすべての条件が固定されている場合にのみ意味があります。 [C++の解説書](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/cpp)では、この点についてさらに詳しく論じています。なぜなら、C++の解説書ではリカバリではなくポストセレクションを採用しており、またその他のパラメータにおいても異なるからです。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8276be0a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Parsing assets/sqdrift/fcidump_files/N2_sto_3g\n",
            "Reference FCI Energy  = -107.6481842917 Ha\n",
            "Nuclear Repulsion Energy = 23.7887003074 Ha\n",
            "5060\n",
            "400\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": []
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Selected backend: ibm_aachen (156 qubits)\n",
            "40000\n",
            "Electron configuration:\n",
            "  Total electrons: 14\n",
            "  Alpha electrons: 7\n",
            "  Beta electrons: 7\n",
            "  Number of orbitals: 10\n",
            "  Number of spin orbitals (qubits): 20\n",
            "Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)\n",
            "\n",
            "Running SQD with configuration recovery...\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64593072647523\n",
            "\t\tSubspace dimension: 7221\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.6458270048177\n",
            "\t\tSubspace dimension: 7209\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64007673117075\n",
            "\t\tSubspace dimension: 7138\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64757372124944\n",
            "\t\tSubspace dimension: 9009\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64674060104392\n",
            "\t\tSubspace dimension: 8245\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64731360491942\n",
            "\t\tSubspace dimension: 8178\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64765518770588\n",
            "\t\tSubspace dimension: 8835\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64767975712016\n",
            "\t\tSubspace dimension: 8649\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64761634415606\n",
            "\t\tSubspace dimension: 8648\n",
            "FINAL SQD RESULTS\n",
            "Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295\n",
            " 0.99298797 0.03728154 0.03728186 0.00938359]\n",
            "Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174\n",
            " 0.99298424 0.03728121 0.03728169 0.00939159]\n",
            "Reference Energy: -107.6481842917 Ha\n",
            "Computed Energy:  -107.6476797571 Ha\n",
            "Error:            5.0453460619e-04 Ha\n"
          ]
        }
      ],
      "source": [
        "name = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "fcidump = tools.fcidump.read(name)\n",
        "\n",
        "# Extract metadata from the FCIDump header\n",
        "norb = fcidump[\"NORB\"]  # number of spatial orbitals\n",
        "nelec = fcidump[\"NELEC\"]  # total number of electrons\n",
        "e_nuc = fcidump[\"ECORE\"]  # nuclear repulsion / core energy\n",
        "ms2 = fcidump[\"MS2\"]  # 2S (spin)\n",
        "\n",
        "num_elec_a = (nelec + ms2) // 2  # alpha electrons\n",
        "num_elec_b = (nelec - ms2) // 2  # beta  electrons\n",
        "\n",
        "# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)\n",
        "h1e = fcidump[\"H1\"]  # shape (norb, norb)\n",
        "h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)\n",
        "    1, fcidump[\"H2\"], norb\n",
        ")\n",
        "\n",
        "cisolver = fci.direct_spin1.FCI()\n",
        "cisolver.max_cycle = 200\n",
        "cisolver.conv_tol = 1e-12\n",
        "\n",
        "e_fci, _ = cisolver.kernel(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    norb,\n",
        "    (num_elec_a, num_elec_b),\n",
        "    ecore=e_nuc,  # adds nuclear repulsion to the final energy\n",
        ")\n",
        "\n",
        "reference_energy = e_fci\n",
        "\n",
        "print(f\"Reference FCI Energy  = {reference_energy:.10f} Ha\")\n",
        "\n",
        "nuclear_repulsion_energy = fcidump[\"ECORE\"]\n",
        "print(f\"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha\")\n",
        "\n",
        "fcidump = FCIDump.from_file(name)\n",
        "hamiltonian = FermionOperator.from_fcidump(fcidump)\n",
        "num_modes = 2 * fcidump.norb\n",
        "\n",
        "# Apply automatic grouping\n",
        "canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))\n",
        "exit_code = group_terms_by_electronic_structure(\n",
        "    canon, num_modes, two_body_physicist_order=False\n",
        ")\n",
        "filter_diagonal_terms(canon)\n",
        "\n",
        "print(len(canon.groups))\n",
        "\n",
        "# SqDRIFT parameters\n",
        "times = [1.0, 10.0]  # Total evolution times used for the subspace creation\n",
        "num_exc = 10  # Number of excitation groups per circuit\n",
        "num_circuits = 200  # Number of circuits to generate\n",
        "\n",
        "init_circuits = []\n",
        "hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))\n",
        "\n",
        "for time in times:\n",
        "    evo_gate = Evolution(num_modes, canon, time)\n",
        "    circ = FermionicCircuit(num_modes)\n",
        "    circ.append(hf_gate, circ.modes)\n",
        "    circ.append(evo_gate, circ.modes)\n",
        "    init_circuits.append(circ)\n",
        "\n",
        "# Calculate weights for sampling (one per group)\n",
        "qdrift = QDriftTrotterization(num_exc, rng=19)\n",
        "\n",
        "pm = generate_preset_jw_pass_manager()\n",
        "pm.optimization = FermionicPassManager([qdrift])\n",
        "\n",
        "sqdrift_circuits = []\n",
        "for circ in init_circuits:\n",
        "    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))\n",
        "\n",
        "for circ in sqdrift_circuits:\n",
        "    circ.measure_all()\n",
        "\n",
        "print(len(sqdrift_circuits))\n",
        "\n",
        "# This example assumes you have saved your IBM Quantum Platform account locally.\n",
        "service = QiskitRuntimeService(channel=\"ibm_quantum_platform\")\n",
        "\n",
        "# Select backend (choose based on qubit requirements)\n",
        "backend = service.least_busy(\n",
        "    operational=True,\n",
        "    simulator=False,\n",
        "    min_num_qubits=2 * norb,\n",
        ")\n",
        "\n",
        "print(f\"Selected backend: {backend.name} ({backend.num_qubits} qubits)\")\n",
        "\n",
        "# Transpile for hardware\n",
        "transpiled_circuits = transpile(\n",
        "    sqdrift_circuits,\n",
        "    backend=backend,\n",
        "    optimization_level=3,\n",
        "    seed_transpiler=42,\n",
        ")\n",
        "\n",
        "shots = 100\n",
        "\n",
        "sampler = Sampler(mode=backend)\n",
        "\n",
        "sampler.options.environment.job_tags = [\"TUT-SqDRIFT\"]\n",
        "\n",
        "job = sampler.run(transpiled_circuits, shots=shots)\n",
        "result = job.result()\n",
        "\n",
        "# Extract counts from SamplerV2 results\n",
        "all_counts = [pub_result.data.meas.get_counts() for pub_result in result]\n",
        "\n",
        "# Set to True to filter out bitstrings that violate electron-number conservation\n",
        "PRUNE = False\n",
        "\n",
        "\n",
        "def is_valid_bitstring(\n",
        "    bitstring: str, norb: int, nelec: tuple[int, int]\n",
        ") -> bool:\n",
        "    n_alpha, n_beta = nelec\n",
        "    return (\n",
        "        len(bitstring) == 2 * norb\n",
        "        and bitstring[norb:].count(\"1\") == n_alpha\n",
        "        and bitstring[:norb].count(\"1\") == n_beta\n",
        "    )\n",
        "\n",
        "\n",
        "if PRUNE:\n",
        "    all_counts_filtered = []\n",
        "    for counts in all_counts:\n",
        "        filtered_count = {}\n",
        "        for key in counts:\n",
        "            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):\n",
        "                continue\n",
        "            elif key not in filtered_count.keys():\n",
        "                filtered_count[key] = counts[key]\n",
        "            else:\n",
        "                filtered_count[key] += counts[key]\n",
        "        all_counts_filtered.append(filtered_count)\n",
        "    all_counts = all_counts_filtered\n",
        "\n",
        "combined_counts = {}\n",
        "for counts in all_counts:\n",
        "    for bitstring, count in counts.items():\n",
        "        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count\n",
        "\n",
        "bit_array = BitArray.from_counts(combined_counts)\n",
        "print(bit_array.num_shots)\n",
        "\n",
        "print(\"Electron configuration:\")\n",
        "print(f\"  Total electrons: {nelec}\")\n",
        "print(f\"  Alpha electrons: {num_elec_a}\")\n",
        "print(f\"  Beta electrons: {num_elec_b}\")\n",
        "print(f\"  Number of orbitals: {norb}\")\n",
        "print(f\"  Number of spin orbitals (qubits): {2*norb}\")\n",
        "print(f\"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}\")\n",
        "\n",
        "# SQD parameters\n",
        "samples_per_batch = 300\n",
        "num_batches = 3\n",
        "max_iterations = 5\n",
        "\n",
        "initial_occupancies = (\n",
        "    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha\n",
        "    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta\n",
        ")\n",
        "\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "# Run SQD with configuration recovery\n",
        "print(\"\\nRunning SQD with configuration recovery...\")\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=1e-3,\n",
        "    occupancies_tol=1e-3,\n",
        "    max_iterations=max_iterations,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    seed=42,\n",
        "    callback=callback,\n",
        ")\n",
        "\n",
        "computed_energy = result.energy + nuclear_repulsion_energy\n",
        "\n",
        "print(\"FINAL SQD RESULTS\")\n",
        "print(f\"Orbital occupancies (alpha): {result.orbital_occupancies[0]}\")\n",
        "print(f\"Orbital occupancies (beta): {result.orbital_occupancies[1]}\")\n",
        "\n",
        "\n",
        "energy_error = abs(computed_energy - reference_energy)\n",
        "print(f\"Reference Energy: {reference_energy:.10f} Ha\")\n",
        "print(f\"Computed Energy:  {computed_energy:.10f} Ha\")\n",
        "print(f\"Error:            {energy_error:.10e} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e18de7d4",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## 次のステップ\n",
        "\n",
        "<Admonition type=\"note\" title=\"推奨事項\">\n",
        "  この作品に興味を持たれた方は、以下の資料もご参照ください：\n",
        "\n",
        "  * [フェルミオン格子モデルのサンプルベース・クリロフ量子対角化](/docs/tutorials/sample-based-krylov-quantum-diagonalization) ― 変分アンザッツの代わりに時間発展回路を用いた関連チュートリアル。\n",
        "  * [化学ハミルトニアンのサンプルベース量子対角化](/docs/tutorials/sample-based-quantum-diagonalization) ― 量子化学シミュレーションのための局所ユニタリクラスター・ジャストロウ（LUCJ）回路の構築方法に関するチュートリアル。\n",
        "  * 論文『 [SqDRIFT](https://arxiv.org/abs/2508.02578) 』――本チュートリアルが基にしている文献。 （なお、本論文で取り上げている最適化の一部は現在開発中の段階であり、このチュートリアルは、使用しているライブラリの進化に伴い、将来変更される可能性があります。）\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}