{ "cells": [ { "cell_type": "markdown", "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e", "metadata": {}, "source": [ "---\n", "title: \"化学ハミルトニアンのサンプル・ベースの量子対角化\"\n", "description: \"ノイズの多い量子ハードウェアを用いて窒素分子をシミュレートするために、サンプルベース量子対角化アルゴリズムを使用する。\"\n", "---\n", "\n", "{/* cspell:ignore fontdict fontsize milli LUCJ CCSD ccsd hcore pvdz */}\n", "\n", "\n", "\n", "# 化学ハミルトニアンのサンプル・ベースの量子対角化\n", "\n", "*使用時間の見積もり：Heron r2 プロセッサーで1分未満（注：これはあくまでも見積もりです。 ランタイムは異なるかもしれない)。*\n", "\n" ] }, { "cell_type": "markdown", "id": "9701917b", "metadata": {}, "source": [ "\n", "\n", "## 学習成果\n", "\n", "このチュートリアルを終えた後、ユーザーは以下の点を理解できるようになります：\n", "\n", "* [SQD Qiskit アドオン](/docs/addons/qiskit-addon-sqd)を使用して、量子処理ユニット（QPU）からサンプリングしたビット列を用いて、分子系の基底状態エネルギーを近似する方法。\n", "* [ffsim](https://github.com/qiskit-community/ffsim) を使用して、量子化学シミュレーション用の局所ユニタリークラスター・ジャストロウ（LUCJ）回路を構築する方法。\n", "\n", "\n", "\n", "## 前提条件\n", "\n", "このチュートリアルを進める前に、以下のトピックについて事前に確認しておくことをお勧めします：\n", "\n", "* 量子化学と第二量子化\n", "* サンプラープリミティブを使用した量子回路からのサンプリング\n", "\n", "\n", "\n", "## 背景\n", "\n", "このチュートリアルでは、 [SQD Qiskit アドオン](https://github.com/Qiskit/qiskit-addon-sqd)を使用して[サンプルベース量子対角化（SQD）アルゴリズム](https://arxiv.org/abs/2405.05068)を実装し、ノイズを含む量子サンプルを後処理して、平衡結合長における窒素分子 $\\text{N}_2$ の基底状態を近似する方法を解説します。 このソフトウェアの詳細については、対応する[ドキュメント](/docs/addons/qiskit-addon-sqd)に記載されています。そこには、使い始めるための[簡単な例](/docs/addons/qiskit-addon-sqd/guides/quickstart)も掲載されています。\n", "\n", "このチュートリアルは、量子化学、特に分子の基底状態のエネルギーを求める方法に精通しているユーザーにおすすめです。 ワークフローの詳細な手順については、 [「量子対角化アルゴリズム」コース](/learning/courses/quantum-diagonalization-algorithms)を参照してください。\n", "\n", "SQDは、量子計算と分散型古典計算を組み合わせて、量子系のハミルトニアンなどの量子演算子の固有値と固有ベクトルを求める手法である。 従来の分散計算は、量子プロセッサから得られたサンプルを処理し、それらが張る部分空間において対象のハミルトニアンを射影および対角化するために用いられる。 SQD ベースのワークフローには、以下の手順が含まれます：\n", "\n", "1. 回路アナザッツを選び、量子コンピュータ上で参照状態（この場合は[ハートリーフォック](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method)状態）に適用する。\n", "2. 出来上がった量子状態からビット列をサンプリングする。\n", "3. ビット列に対して*自己一貫性のある構成復元*手順を実行し、基底状態の近似を求める。\n", "\n", "SQDは、対象となる固有状態がスパースである場合にうまく機能することが知られています。波動関数は基底状態のセット $\\mathcal{S} = \\{|x\\rangle \\}$、そのサイズは問題の大きさに対して指数関数的に増加することはありません。\n", "\n", "\n", "\n", "### 量子化学\n", "\n", "分子系のハミルトニアンは次のように書ける\n", "\n", "$$\n", "\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n", "+ \\frac12\n", "\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n", "h_{prqs} \\,\n", "\\hat{a}^\\dagger_{p\\sigma}\n", "\\hat{a}^\\dagger_{q\\tau}\n", "\\hat{a}_{s\\tau}\n", "\\hat{a}_{r\\sigma},\n", "$$\n", "\n", "ここで、 $h_{pr}$ と $h_{prqs}$ は分子積分と呼ばれる複素数で、コンピュータ・プログラムを使って分子の仕様から計算できる。 このチュートリアルでは [PySCF](https://pyscf.org/) ソフトウェアパッケージを使用して積分を計算します。\n", "\n", "分子ハミルトニアンの導出方法の詳細については、量子化学の教科書（例えば、SzaboとOstlund著『 *Modern Quantum Chemistry* 』）を参照されたい。 量子化学の問題がどのように量子コンピュータにマッピングされるかについての高レベルの説明は、Qiskit Global Summer School 2024の講義 [*Mapping Problems to Qubitsを*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900)ご覧ください。\n", "\n", "\n", "\n", "### 局所ユニタリークラスター・ジャストロー（LUCJ）近似\n", "\n", "SQDでは、サンプルを抽出するための量子回路のアンザッツが必要です。 このチュートリアルでは、物理的な根拠とハードウェアへの適応性を兼ね備えていることから、 [ローカルユニタリークラスター・ジャストロウ（LUCJ）](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) アンザッツを用います。 [ffsim](https://qiskit-community.github.io/ffsim/) を使用して、近似回路を構築します。\n", "\n", "LUCJアプローチは、量子ビット間の接続性が制限されているQPUにも適応する。 スピン軌道は、SWAPゲートによるルーティングを必要としないような形で量子ビットにマッピングされる。 IBM® このハードウェアは、ヘビー・ヘックス格子型の量子ビットトポロジーを採用しており、その場合は、下図に示すような「ジグザグ」パターンを採用することができます。 このパターンでは、同じスピンを持つ軌道は線形トポロジーを持つ量子ビット（赤と青の円）にマッピングされ、異なるスピンを持つ軌道間の接続は空間上の4つおきの軌道に存在し、その接続は補助量子ビット（紫の円）によって実現されている。\n", "\n", "![重ヘックス格子上のLUCJアサッツのQubit写像図](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/7e0ee7e1-2d24-417f-ac59-25c58db79aa9.avif)\n", "\n", "\n", "\n", "### 自己整合的な構成復旧\n", "\n", "自己無撞着構成回復手順は、ノイズの多い量子サンプルから可能な限り多くの信号を抽出するように設計されている。 分子ハミルトニアンは粒子数とスピンZを保存するので、これらの対称性も保存する回路アサッツを選ぶことは理にかなっている。 ハートリーフォック状態に適用すると、結果として得られる状態は、ノイズのない設定において固定された粒子数とスピンZを持つ。 したがって、この状態からサンプリングされたビット列のスピン $\\alpha$ とスピン $\\beta$ の半分は、ハートリーフォック状態と同じ[ハミングウェイトを](https://en.wikipedia.org/wiki/Hamming_weight)持つはずである。 現在の量子プロセッサーにはノイズが存在するため、測定されたビット列の中にはこの性質に反するものがある。 単純なポストセレクションでは、これらのビット列は破棄されるが、これは無駄である。 自己無撞着な回復手順は、後処理でその信号の一部を回復しようとする。 この手順は反復的であり、基底状態における各軌道の平均占有率の推定値を入力として必要とする。 この手順はループで実行され、各反復には以下のステップがある：\n", "\n", "1. 指定された対称性に違反する各ビット列について、ビット列を平均軌道占有率の現在の推定値に近づけるように設計された確率的手順でビットを反転させ、新しいビット列を得る。\n", "2. 対称性を満たす新旧のビット列をすべて集め、あらかじめ選んだ固定サイズの部分集合をサブサンプルする。\n", "3. ビット列の各サブセットについて、ハミルトニアンを対応する基底ベクトル（これらの基底ベクトルの説明については[前節を](#quantum-chemistry)参照）によってスパンされる部分空間に射影し、古典コンピュータ上で射影されたハミルトニアンの基底状態推定値を計算する。\n", "4. 平均軌道占有率の推定値を、最もエネルギーの低い基底状態の推定値で更新する。\n", "\n", "\n", "\n", "### SQDワークフロー図\n", "\n", "SQDのワークフローは以下の図のようになっている：\n", "\n", "![SQDアルゴリズムのワークフロー図](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/fd7e816f-4e2e-4dd7-a7da-f71afb9ca68d.avif)\n", "\n" ] }, { "cell_type": "markdown", "id": "88422c4b", "metadata": {}, "source": [ "\n", "\n", "## 要件\n", "\n", "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n", "\n", "* Qiskit SDK v1.0 またはそれ以降、 [可視化](/docs/api/qiskit/visualization)サポート付き\n", "* Qiskit Runtime v0.22 またはそれ以降 (`pip install qiskit-ibm-runtime`)\n", "* SQD Qiskit アドオン v0.11 以降 (`pip install qiskit-addon-sqd`)\n", "* ffsim v0.0.75 以降 (`pip install ffsim`)\n", "\n" ] }, { "cell_type": "markdown", "id": "c6e44a31", "metadata": {}, "source": [ "\n", "\n", "## セットアップ\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "6e51c3d8", "metadata": {}, "outputs": [], "source": [ "import math\n", "\n", "import ffsim\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pyscf\n", "import pyscf.cc\n", "import pyscf.mcscf\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "from qiskit.primitives import StatevectorSampler\n", "from qiskit.providers.fake_provider import GenericBackendV2\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "from qiskit_ibm_runtime import SamplerV2 as Sampler" ] }, { "cell_type": "markdown", "id": "4bc6ee26-4371-4cd2-80a7-60752bf8775d", "metadata": {}, "source": [ "\n", "\n", "## 小規模シミュレータの例\n", "\n", "このチュートリアルでは、平衡結合距離付近にある窒素分子の基底状態の近似解を求めます。 まず、実験をシミュレーションして正常に動作することを確認するため、小規模な STO-6G 基底セットを使用します。\n", "\n" ] }, { "cell_type": "markdown", "id": "afeb054c", "metadata": {}, "source": [ "\n", "\n", "### ステップ1：古典的な入力を量子問題にマッピングする\n", "\n", "まず、分子とその特性を特定します。\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "b821e660", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -108.464957764796\n", "CASCI E = -108.595987350986 E(CI) = -32.4115475088426 S^2 = 0.0000000\n", "norb = 8\n", "nelec = (5, 5)\n" ] } ], "source": [ "# Specify molecule properties\n", "spin_sq = 0\n", "\n", "# Build N2 molecule\n", "mol = pyscf.gto.Mole()\n", "mol.build(\n", " atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n", " basis=\"sto-6g\",\n", " symmetry=\"Dooh\",\n", ")\n", "\n", "# Define active space\n", "n_frozen = 2\n", "active_space = range(n_frozen, mol.nao_nr())\n", "\n", "# Get molecular integrals\n", "scf = pyscf.scf.RHF(mol).run()\n", "norb = len(active_space)\n", "n_electrons = int(sum(scf.mo_occ[active_space]))\n", "n_alpha = (n_electrons + mol.spin) // 2\n", "n_beta = (n_electrons - mol.spin) // 2\n", "nelec = (n_alpha, n_beta)\n", "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n", "mo = cas.sort_mo(active_space, base=0)\n", "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n", "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n", "\n", "# Compute exact energy using FCI\n", "reference_energy = cas.run().e_tot\n", "\n", "print(f\"norb = {norb}\")\n", "print(f\"nelec = {nelec}\")" ] }, { "cell_type": "markdown", "id": "96bfe018", "metadata": {}, "source": [ "LUCJアンサッツ回路を構成する前に、まず以下のコードセルでCCSD計算を行う。 この計算で得られた [$t_1$ と $t_2$ の振幅は](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator)、ansatzのパラメータを初期化するのに使われる。\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "efe83d98", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "E(CCSD) = -108.5933309085008 E_corr = -0.1283731437052354\n" ] } ], "source": [ "# Get CCSD t2 amplitudes for initializing the ansatz\n", "ccsd = pyscf.cc.CCSD(\n", " scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n", ").run()\n", "t1 = ccsd.t1\n", "t2 = ccsd.t2" ] }, { "attachments": {}, "cell_type": "markdown", "id": "f4d882fa", "metadata": {}, "source": [ "ここで、 [ffsim](https://github.com/qiskit-community/ffsim) を使用してアンザッツ回路を作成します。 我々の分子は閉殻のハートリー・フォック状態を持つため、UCJアンザッツのスピン平衡型、すなわち [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced) を用いる。 $t_2$ の振幅の「圧縮」二重因数分解を有効にするため、`method` `from_t_amplitudes` メソッドに `compressed` を設定 `optimize=True` しました（詳細については、ffsimのドキュメントにある [「The local unitary cluster Jastrow (LUCJ) ansatz」](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) を参照してください）。\n", "\n", "LUCJアンザッツはQPUの利用可能な接続性に適応するため、アンザッツを作成する前にQPUバックエンドを初期化する必要があります。 とりあえず、ヘビー・ヘックス結合マップと、LUCJアンザッツによって自然に分解されるゲートセットを持つ汎用的なバックエンドを作成します。 次に、 [LUCJアプローチに関する背景セクション](#local-unitary-cluster-jastrow-lucj-ansatz)で説明した「ジグザグ」レイアウトに従って、LUCJアプローチを指定されたバックエンドへトランスパイルするための専用パスマネージャーを作成 `ffsim.qiskit.generate_lucj_pass_manager` します。 この関数は、選択されたレイアウトに関連する誤差を最小化するためにスコアリング・ヒューリスティックを使用します。これは、バックエンドが実際のQPUである場合や、ノイズモデルを備えたシミュレータである場合に重要です。 この関数は、パスマネージャーを返すだけでなく、ハードウェア上で実装可能なアルファ・ベータ結合ペアも返します。 すべてのペアを実装できない場合、警告が出力されます。\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dd69a86c", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "from qiskit.transpiler import CouplingMap\n", "\n", "warnings.formatwarning = lambda msg, *args, **kwargs: f\"Warning: {msg}\\n\"\n", "\n", "# Set ansatz properties\n", "n_reps = 1\n", "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n", "\n", "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n", "pairs_ab = None\n", "\n", "# Initialize backend\n", "coupling_map = CouplingMap.from_heavy_hex(3)\n", "backend = GenericBackendV2(\n", " coupling_map.size(),\n", " coupling_map=coupling_map,\n", " basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\", \"swap\"],\n", ")\n", "\n", "# Create pass manager\n", "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n", " backend=backend,\n", " norb=norb,\n", " connectivity=\"heavy-hex\",\n", " interaction_pairs=(pairs_aa, pairs_ab),\n", " optimization_level=3,\n", ")\n", "\n", "# Create the LUCJ ansatz operator\n", "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n", " t2=t2,\n", " t1=t1,\n", " n_reps=n_reps,\n", " interaction_pairs=(pairs_aa, pairs_ab),\n", " # Setting optimize=True enables the \"compressed\" factorization\n", " optimize=True,\n", " # Limit the number of optimization iterations to prevent the code cell\n", " # from running too long. Removing this line may improve results.\n", " options=dict(maxiter=1000),\n", ")\n", "\n", "# create an empty quantum circuit\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "\n", "# prepare Hartree-Fock state as the reference state and append it\n", "# to the quantum circuit\n", "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n", "\n", "# apply the UCJ operator to the reference state\n", "circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n", "circuit.measure_all()" ] }, { "cell_type": "markdown", "id": "db11bf6d", "metadata": {}, "source": [ "\n", "\n", "### ステップ2：量子ハードウェアでの実行に向けた最適化\n", "\n", "次に、対象のハードウェアに合わせて回路を最適化します。 通常、このステップでは、ハードウェアバックエンドとそのバックエンド用のパスマネージャーを初期化します。 ただし、LUCJアプローチはハードウェアの接続性に適合しているため、これらの処理は前のステップですでに実施済みです。 あとは、回路上でパスマネージャーを実行し、QPU上で直接実行可能なISA回路にトランスパイルするだけです。\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "7d554aa5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, 'barrier': 1})\n" ] } ], "source": [ "isa_circuit = pass_manager.run(circuit)\n", "print(f\"Gate counts: {isa_circuit.count_ops()}\")" ] }, { "cell_type": "markdown", "id": "0cc1edef", "metadata": {}, "source": [ "\n", "\n", "### ステップ3: `Qiskit primitives`を使用して実行する\n", "\n", "ハードウェア実行向けに回路を最適化した後、ターゲットハードウェア上で実行し、基底状態のエネルギー推定のためのサンプルを収集する準備が整いました。 回路は1つしかないため、 IBM Quantum の「Compute Service」 [ジョブ実行モード](/docs/guides/execution-modes)を使用して、この回路を実行します。\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout\n" ] } ], "source": [ "rng = np.random.default_rng()\n", "sampler = StatevectorSampler(seed=rng)\n", "job = sampler.run([isa_circuit], shots=100_000)" ] }, { "cell_type": "code", "execution_count": 7, "id": "332ecab3-77e6-473f-b0e7-af30f983393a", "metadata": {}, "outputs": [], "source": [ "primitive_result = job.result()\n", "pub_result = primitive_result[0]" ] }, { "cell_type": "markdown", "id": "6df05b6e", "metadata": {}, "source": [ "\n", "\n", "### ステップ4：後処理を行い、結果を希望の古典形式で返す\n", "\n", "QPU出力の品質を判断する有用な指標は、返される有効な構成の数である。 有効な構成は正しい粒子数とスピンZを持つ。これはビット列の右半分がスピンアップ電子の数に等しいハミング重量を持ち、左半分がスピンダウン電子の数に等しいハミング重量を持つことを意味する。 以下のセルは、サンプリングされた構成のうち有効なものの割合を計算します。\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "718f8517", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fraction of sampled configurations that are valid: 1.0\n" ] } ], "source": [ "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", "bit_array = pub_result.data.meas\n", "num_valid = sum(\n", " is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n", ")\n", "valid_fraction = num_valid / bit_array.num_shots\n", "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")" ] }, { "cell_type": "markdown", "id": "7b126f3a", "metadata": {}, "source": [ "ノイズのないシミュレータで回路をサンプリングしているため、すべてのビット列は有効です。 ノイズの多いQPU上で実行した場合、この比率は1未満になりますが、ビット列が均一にランダムにサンプリングされた場合に予想される比率よりは大きくなることを期待しています。その値は、次のセルで計算されます。\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6b3e4bca", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625\n" ] } ], "source": [ "expected_fraction_random = (\n", " math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n", ")\n", "print(\n", " f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n", " f\"{expected_fraction_random}\"\n", ")" ] }, { "cell_type": "markdown", "id": "eb704101-0fe8-4d12-b572-b1d844e35a90", "metadata": {}, "source": [ "ここで、 `diagonalize_fermionic_hamiltonian` 関数を用いてハミルトニアンの基底状態エネルギーを推定する。 この関数は、エネルギー推定値を改善するために、ノイズの多い量子サンプルを繰り返し精緻化する自己無撞着構成回復手順を実行する。 コールバック関数を渡すことで、後の分析のために中間結果を保存することができる。 `diagonalize_fermionic_hamiltonian` への引数の説明については、 [APIドキュメントを](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian)参照のこと。\n", "\n", "ここでは、基底状態における軌道占有の初期推定値としてハートリー・フォック配置を指定するために、引数 を用いて `diagonalize_fermionic_hamiltonian``initial_occupancies` を指定する。 この手法は、基底状態がハートリー・フォック配置に顕著な支持を持つシステムでは合理的であるが、他の状況では適切でない可能性がある。ただし、より高度な計算手法を用いれば、そうしたケースでもより良い初期推定値が得られるかもしれない。 Specifying `initial_occupancies` は、ノイズの多い QPU 上で大規模な回路をサンプリングする場合など、有効な構成がサンプリングされなかった場合でも構成復元を実行できるようにします。 この引数がなければ、有効な構成が提供されていない場合、構成の復元は失敗しエラーが発生します。\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "2f32a352", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "\tSubsample 1\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "\tSubsample 2\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "\tSubsample 1\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "\tSubsample 2\n", "\t\tEnergy: -108.59275573641656\n", "\t\tSubspace dimension: 900\n", "Final energy: -108.59275573641656\n", "Final energy error: 0.0032316145694579745\n" ] } ], "source": [ "from functools import partial\n", "\n", "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", " solve_sci_batch,\n", ")\n", "\n", "# SQD options\n", "energy_tol = 1e-3\n", "occupancies_tol = 1e-3\n", "max_iterations = 5\n", "\n", "# Eigenstate solver options\n", "num_batches = 3\n", "samples_per_batch = 300\n", "symmetrize_spin = True\n", "carryover_threshold = 1e-4\n", "max_cycle = 200\n", "\n", "# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n", "initial_occupancies = (\n", " np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n", " np.array([1] * n_beta + [0] * (norb - n_beta)),\n", ")\n", "\n", "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n", "# you can omit this step, in which case you would not specify the sci_solver argument\n", "# in the call to diagonalize_fermionic_hamiltonian below.\n", "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n", "\n", "# List to capture intermediate results\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", "result = diagonalize_fermionic_hamiltonian(\n", " hcore,\n", " eri,\n", " bit_array,\n", " samples_per_batch=samples_per_batch,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=num_batches,\n", " energy_tol=energy_tol,\n", " occupancies_tol=occupancies_tol,\n", " max_iterations=max_iterations,\n", " sci_solver=sci_solver,\n", " symmetrize_spin=symmetrize_spin,\n", " initial_occupancies=initial_occupancies,\n", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=rng,\n", ")\n", "\n", "final_energy = result.energy + nuclear_repulsion_energy\n", "energy_error = final_energy - reference_energy\n", "print(f\"Final energy: {final_energy}\")\n", "print(f\"Final energy error: {energy_error}\")" ] }, { "cell_type": "markdown", "id": "9d78906b-4759-4506-9c69-85d4e67766b3", "metadata": {}, "source": [ "\n", "\n", "#### 結果を可視化する\n", "\n", "最初のグラフは、このシミュレーションにおいて、最初の反復計算の後ですでに正確な解 `1 mH` に極めて近い値を得ていることを示しています（化学的精度は通常、 `1 kcal/mol`$\\approx$ 程度と `1.6 mH`みなされます）。ただし、これは小規模な系であり、サンプルにノイズが含まれていないため、構成の復元は必要ありません。 ノイズの多いQPU上で動作する大規模なシステムでは、構成の復元を複数回繰り返す必要が生じ、最終的な精度が低下する可能性があります。 一般的に、構成の復元反復回数を増やすか、バッチあたりのサンプル数を増やすことで、エネルギーを改善することができます。\n", "\n", "2番目のプロットは、最終反復後の各空間軌道の平均占有率を示している。 スピン・アップ電子とスピン・ダウン電子の両方が、解の最初の5つの軌道を高い確率で占めていることがわかる。\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Data for energies plot\n", "x1 = range(len(result_history))\n", "min_e = [\n", " min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n", " for result in result_history\n", "]\n", "e_diff = [abs(e - reference_energy) for e in min_e]\n", "yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n", "\n", "# Chemical accuracy (+/- 1 milli-Hartree)\n", "chem_accuracy = 0.001\n", "\n", "# Data for avg spatial orbital occupancy\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "\n", "# Plot energies\n", "axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].set_yticks(yt1)\n", "axs[0].set_yticklabels(yt1)\n", "axs[0].set_yscale(\"log\")\n", "axs[0].set_ylim(1e-4)\n", "axs[0].axhline(\n", " y=chem_accuracy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"chemical accuracy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "\n", "# Plot orbital occupancy\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ce0eecb3-8a23-4118-aa1e-a28afcec6334", "metadata": {}, "source": [ "\n", "\n", "## 大規模なハードウェアの例\n", "\n", "それでは、実際の量子ハードウェアを使って、より大規模な例を実行してみましょう。 ここでは、 cc-pVDZ 基底関数セットを用いて、窒素分子のアクティブ空間を導出する。\n", "\n" ] }, { "cell_type": "markdown", "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062", "metadata": {}, "source": [ "\n", "\n", "### 手順 1～4\n", "\n", "ここでは、これらすべての手順を統合して大規模な単一のワークフローとし、それを実際の量子ハードウェア上で実行します。\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3858949c-a55d-4ff8-a0fc-54fb53e131b5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -108.929838385609\n", "norb = 26\n", "nelec = (5, 5)\n", "E(CCSD) = -109.2177884185544 E_corr = -0.2879500329450045\n", "Using backend ibm_boston\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].\n", "Removing interaction (24, 24) from the end.\n", "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].\n", "Removing interaction (20, 20) from the end.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Gate counts: OrderedDict({'sx': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})\n", "Fraction of sampled configurations that are valid: 0.02124\n", "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n", "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -109.13889134249762\n", "\t\tSubspace dimension: 120409\n", "\tSubsample 1\n", "\t\tEnergy: -109.11785470455858\n", "\t\tSubspace dimension: 110889\n", "\tSubsample 2\n", "\t\tEnergy: -109.13234360554011\n", "\t\tSubspace dimension: 130321\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -109.16392179579177\n", "\t\tSubspace dimension: 223729\n", "\tSubsample 1\n", "\t\tEnergy: -109.16281938332986\n", "\t\tSubspace dimension: 223729\n", "\tSubsample 2\n", "\t\tEnergy: -109.16955816711932\n", "\t\tSubspace dimension: 233289\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -109.17905772999075\n", "\t\tSubspace dimension: 324900\n", "\tSubsample 1\n", "\t\tEnergy: -109.17532445048462\n", "\t\tSubspace dimension: 357604\n", "\tSubsample 2\n", "\t\tEnergy: -109.1733168689756\n", "\t\tSubspace dimension: 348100\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -109.18437778820451\n", "\t\tSubspace dimension: 474721\n", "\tSubsample 1\n", "\t\tEnergy: -109.18450164209159\n", "\t\tSubspace dimension: 476100\n", "\tSubsample 2\n", "\t\tEnergy: -109.18493571190754\n", "\t\tSubspace dimension: 487204\n", "Iteration 5\n", "\tSubsample 0\n", "\t\tEnergy: -109.18616522497996\n", "\t\tSubspace dimension: 622521\n", "\tSubsample 1\n", "\t\tEnergy: -109.18652868888333\n", "\t\tSubspace dimension: 644809\n", "\tSubsample 2\n", "\t\tEnergy: -109.18753326484406\n", "\t\tSubspace dimension: 585225\n", "Final energy: -109.18753326484406\n", "Final energy error: 0.040495951813099396\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# ------------------------------ Step 1 ------------------------------\n", "# Build N2 molecule\n", "mol = pyscf.gto.Mole()\n", "mol.build(\n", " atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n", " basis=\"cc-pvdz\",\n", " symmetry=\"Dooh\",\n", ")\n", "\n", "# Define active space\n", "n_frozen = 2\n", "active_space = range(n_frozen, mol.nao_nr())\n", "\n", "# Get molecular integrals\n", "scf = pyscf.scf.RHF(mol).run()\n", "norb = len(active_space)\n", "n_electrons = int(sum(scf.mo_occ[active_space]))\n", "n_alpha = (n_electrons + mol.spin) // 2\n", "n_beta = (n_electrons - mol.spin) // 2\n", "nelec = (n_alpha, n_beta)\n", "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n", "mo = cas.sort_mo(active_space, base=0)\n", "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n", "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n", "\n", "# Store reference energy from SCI calculation performed separately\n", "reference_energy = -109.22802921665716\n", "\n", "print(f\"norb = {norb}\")\n", "print(f\"nelec = {nelec}\")\n", "\n", "# Get CCSD t2 amplitudes for initializing the ansatz\n", "ccsd = pyscf.cc.CCSD(\n", " scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n", ").run()\n", "t1 = ccsd.t1\n", "t2 = ccsd.t2\n", "\n", "# Set ansatz properties\n", "n_reps = 1\n", "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n", "\n", "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n", "pairs_ab = None\n", "\n", "# Initialize backend\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=133\n", ")\n", "print(f\"Using backend {backend.name}\")\n", "\n", "# Create pass manager\n", "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n", " backend=backend,\n", " norb=norb,\n", " connectivity=\"heavy-hex\",\n", " interaction_pairs=(pairs_aa, pairs_ab),\n", " optimization_level=3,\n", ")\n", "\n", "# Create the LUCJ ansatz operator\n", "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n", " t2=t2,\n", " t1=t1,\n", " n_reps=n_reps,\n", " interaction_pairs=(pairs_aa, pairs_ab),\n", " # Setting optimize=True enables the \"compressed\" factorization\n", " optimize=True,\n", " # Limit the number of optimization iterations to prevent the code cell\n", " # from running too long. Removing this line may improve results.\n", " options=dict(maxiter=1000),\n", ")\n", "\n", "# create an empty quantum circuit\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "\n", "# prepare Hartree-Fock state as the reference state and append it\n", "# to the quantum circuit\n", "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n", "\n", "# apply the UCJ operator to the reference state\n", "circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n", "circuit.measure_all()\n", "\n", "\n", "# ------------------------------ Step 2 ------------------------------\n", "\n", "isa_circuit = pass_manager.run(circuit)\n", "print(f\"Gate counts: {isa_circuit.count_ops()}\")\n", "\n", "\n", "# ------------------------------ Step 3 ------------------------------\n", "sampler = Sampler(mode=backend)\n", "sampler.options.environment.job_tags = [\"TUT_SQD\"]\n", "job = sampler.run([isa_circuit], shots=100_000)\n", "primitive_result = job.result()\n", "pub_result = primitive_result[0]\n", "\n", "\n", "# ------------------------------ Step 4 ------------------------------\n", "\n", "bit_array = pub_result.data.meas\n", "num_valid = sum(\n", " is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n", ")\n", "valid_fraction = num_valid / bit_array.num_shots\n", "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")\n", "expected_fraction_random = (\n", " math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n", ")\n", "print(\n", " f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n", " f\"{expected_fraction_random}\"\n", ")\n", "# SQD options\n", "energy_tol = 1e-3\n", "occupancies_tol = 1e-3\n", "max_iterations = 5\n", "\n", "# Eigenstate solver options\n", "num_batches = 3\n", "samples_per_batch = 300\n", "symmetrize_spin = True\n", "carryover_threshold = 1e-4\n", "max_cycle = 200\n", "\n", "# Use the Hartree-Fock configuration as an initial guess for the\n", "# orbital occupancies\n", "initial_occupancies = (\n", " np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n", " np.array([1] * n_beta + [0] * (norb - n_beta)),\n", ")\n", "\n", "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n", "# you can omit this step, in which case you would not specify the\n", "# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.\n", "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n", "\n", "# List to capture intermediate results\n", "result_history = []\n", "\n", "\n", "result = diagonalize_fermionic_hamiltonian(\n", " hcore,\n", " eri,\n", " bit_array,\n", " samples_per_batch=samples_per_batch,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=num_batches,\n", " energy_tol=energy_tol,\n", " occupancies_tol=occupancies_tol,\n", " max_iterations=max_iterations,\n", " sci_solver=sci_solver,\n", " symmetrize_spin=symmetrize_spin,\n", " initial_occupancies=initial_occupancies,\n", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=rng,\n", ")\n", "\n", "final_energy = result.energy + nuclear_repulsion_energy\n", "energy_error = final_energy - reference_energy\n", "print(f\"Final energy: {final_energy}\")\n", "print(f\"Final energy error: {energy_error}\")\n", "\n", "# Data for energies plot\n", "x1 = range(len(result_history))\n", "min_e = [\n", " min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n", " for result in result_history\n", "]\n", "e_diff = [abs(e - reference_energy) for e in min_e]\n", "yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n", "\n", "# Chemical accuracy (+/- 1 milli-Hartree)\n", "chem_accuracy = 0.001\n", "\n", "# Data for avg spatial orbital occupancy\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "\n", "# Plot energies\n", "axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].set_yticks(yt1)\n", "axs[0].set_yticklabels(yt1)\n", "axs[0].set_yscale(\"log\")\n", "axs[0].set_ylim(1e-4)\n", "axs[0].axhline(\n", " y=chem_accuracy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"chemical accuracy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "\n", "# Plot orbital occupancy\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "405e89ea-57da-4021-bb18-91e8d583d310", "metadata": {}, "source": [ "\n", "\n", "## 次のステップ\n", "\n", "\n", " この作品が興味深かったなら、以下の資料も気に入るかもしれません：\n", "\n", " * [フェルミオン格子モデルのサンプルベース・クリロフ量子対角化法](/docs/tutorials/sample-based-krylov-quantum-diagonalization) ― 変分アンサッツの代わりに時間発展回路を用いた関連チュートリアル\n", " * [DiceソルバーによるSQD化学ワークフローのスケールアップ](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) ― 対角化に、より効率的なDiceソフトウェアを活用する方法を紹介するページ\n", " * [SQD アドオン API ドキュメント](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) - 関 `diagonalize_fermionic_hamiltonian` 数のリファレンス\n", " * [*量子特化型スーパーコンピュータにおける厳密対角化のスケールを超えた化学計算*](https://www.science.org/doi/10.1126/sciadv.adu9991) ― 本チュートリアルが基づく論文\n", "\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" }, "hours": 1.5, "qpuSeconds": 60 }, "nbformat": 4, "nbformat_minor": 5 }