{ "cells": [ { "cell_type": "markdown", "id": "1682996d", "metadata": {}, "source": [ "---\n", "title: \"Qiskit Serverless を使用した暗黙的溶媒計算\"\n", "description: \"SQDおよびIEF-PCMを用いて、量子ハードウェアにおける暗黙的溶媒効果の計算方法を学びます。\"\n", "---\n", "\n" ] }, { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "\n", "\n", "# Qiskit Serverless を使用した暗黙的溶媒計算\n", "\n", "*所要時間の目安：Heron r2 プロセッサで約2分（注：これはあくまで目安です。（実行時間は環境によって異なる場合があります。）*\n", "\n", "{/* cspell:ignore avas AVAS hcore textit TRIC dmas mocore ncore ncas mocas fermilevel ecore orbts iiter IITER edup textcoords xytext fontsize fontweight frameon */}\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "\n", "\n", "## 学習成果\n", "\n", "このチュートリアルを終えた後、ユーザーは以下の点を理解できるようになります：\n", "\n", "* [Qiskit Serverless](/docs/guides/serverless) を使用したリモートワークフローの設定と実行方法\n", "* 量子コンピュータを用いて暗黙的溶媒効果を計算する方法\n", "\n", "\n", "\n", "## 前提条件\n", "\n", "このチュートリアルを進める前に、以下のトピックについて理解しておいていただくことをお勧めします：\n", "\n", "* [サンプルベース量子対角化](/learning/courses/quantum-diagonalization-algorithms/sqd-overview)\n", "\n", "\n", "\n", "## 背景\n", "\n", "計算生物物理学では、暗黙的溶媒法による計算が頻繁に行われている。これらのモデルは、溶媒系を直接モデル化することなく、溶質化合物が溶媒とどのように相互作用するかを記述するものである。その代わりに、溶質系モデルを、経験的に特徴づけられた誘電媒体の数学的表現で包み込むという近似が行われる。この誘電体近似は、直接モデル化された溶質と相互作用する。誘電体は、溶質系の電子場と相互作用することで、その基底状態エネルギーなどの特性に影響を与える。化合物は誘電環境によって挙動が異なるため、これは創薬などで用いられる生物物理モデルにとって重要な点である。空気中（真空）で化合物をモデル化すると、水中でモデル化した場合とは異なる挙動が得られる。医薬品成分は、その大部分が水で構成されている人体に取り込まれる必要があるため、真空状態ではなく、水のような溶液中でその成分をモデル化することが有用である。暗黙的溶媒モデルを用いれば、この挙動を低コストで実現できますが、最終的な結果は、溶質分子と溶媒分子の両方を直接的に表現する、計算コストが高い明示的溶媒モデルに比べて、一般的に近似度が高くなります。\n", "\n", "このチュートリアルでは、量子アルゴリズムであるサンプルベース量子対角化（SQD）を、計算負荷が比較的低い陰的溶媒モデルにどのように組み込むことができるかを解説します。この例では、メチルアミンが水に溶ける際の挙動について説明します。我々は、この量子アルゴリズムをCASCIと呼ばれる最先端の古典的比較手法と比較し、両者の計算結果が極めて一致していることを示す。本稿では、量子コンピューティングを中心とした超小型スーパーコンピューティングアーキテクチャを紹介する。このアーキテクチャでは、ルーチンの量子サンプリング部分における計算負荷の高い従来の後処理を、 Qiskit Serverless 内のクラウド環境にオフロードしている。また、このコードでは、リモートで利用可能なCPUコア間で並列処理を行うことで、計算時間を短縮しています。\n", "\n", "Qiskit Serverless インフラの管理を必要とせずに、分散型の量子および古典的なワークロードを実行するためのフレームワークです。サーバーのプロビジョニング（ EC2s、クラスター、 Docker コンテナの起動）は不要であり、オーケストレーションツール（ Kubernetes、 Docker Swarm）も不要で、監視やメンテナンスも必要ありません。各Serverlessジョブはクリーンなコンテナ内で実行され、コードを実行した後、シャットダウンします。ジョブ間のデータは保持されません。コードを書いて、ジョブを送信するだけです。 Serverlessジョブ内では、プログラムが IBM Quantum® のバックエンドにシームレスにアクセスし、そこで結果の後処理を行うことができます。 Qiskit Serverless を利用することで、ユーザーは常時稼働のリモート CPU コアとメモリにアクセスでき、これにより、従来型のワークロードの一部をリモートリソースに分散させることが可能になります。また、ユーザーはプログラムの並列処理において一定の利点を得られるだけでなく、実行中のデバイスシャットダウンによるよくあるトラブルも回避できます。 Qiskit Serverless の詳細については、 [同](/docs/guides/serverless)製品のドキュメントおよび [GitHub](https://qiskit.github.io/qiskit-serverless/index.html) に掲載されている関連資料をご覧ください。\n", "\n", "このチュートリアルでは、以下の内容の具体的な活用例を紹介します：\n", "\n", "* サンプルベース量子対角化\n", "* 量子コンピューティングのためのクライアント・サーバー型計算モデル\n", "\n", "このチュートリアルは、以下に記載されているクリーブランド・クリニックでの研究に着想を得て、それに基づいて作成されています [Kaliakin, Danil 他「溶媒サンプルに基づく暗黙的量子対角化」『Journal of Physical Chemistry B』 129.23 （2025）：5788-5796](https://pubs.acs.org/doi/10.1021/acs.jpcb.5c01030)。本論文では、暗黙的溶媒計算の全ワークフローを明らかにするとともに、反復的な溶媒自己一貫性法（「ハートウッド・アルゴリズム」、M. Motta, T. Pellegrini (2025)、幾何学的最適化、および量子ビット配置の自動選択。暗黙的溶媒計算を実行するための簡素化されたブラックボックス・インターフェースについては、クリーブランド・クリニックの研究に基づき、同クリニックと IBM® が共同で開発した「 [SQD IEF-PCM Qiskit Function](/docs/guides/function-template-chemistry-workflow) 」テンプレートをご参照ください。\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "metadata": {}, "source": [ "\n", "\n", "## 要件\n", "\n", "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n", "\n" ] }, { "cell_type": "markdown", "id": "9456ea67", "metadata": {}, "source": [ "* Qiskit SDK v2.0 またはそれ以降（ [可視化](/docs/api/qiskit/visualization)機能対応）\n", "* Qiskit Runtime v0.40 またはそれ以降 (`pip install qiskit-ibm-runtime`)\n", "* Qiskit「 IBM 」カタログ `pip install qiskit_ibm_catalog`\n", "* Qiskit、 IBM、サーバーレス `pip install qiskit_serverless`\n", "* Qiskit アドオン：サンプルベース量子対角化 (SQD) v0.12.0 `pip install qiskit_addon_sqd`\n", "* PySCF `pip install pyscf`\n", "* FFSIM `pip install ffsim`\n", "* Matplotlib `pip install matplotlib`\n", "* 幾何学的な `pip install geometric`\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "7db2e559", "metadata": {}, "source": [ "\n", "\n", "## セットアップ\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2c88b910", "metadata": {}, "outputs": [], "source": [ "# Establish Quantum Resource connection\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "service = QiskitRuntimeService()\n", "\n", "backend = service.least_busy()\n", "print(f\"Using backend {backend.name}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "af9286cf", "metadata": {}, "outputs": [], "source": [ "# Establish Classical HPC Resource connection\n", "from qiskit_ibm_catalog import QiskitFunction, QiskitServerless\n", "\n", "client = QiskitServerless()" ] }, { "cell_type": "markdown", "id": "aee647af", "metadata": {}, "source": [ "メインのノートブックプログラムのすぐ隣に、 `source_files`という名前のディレクトリを作成します。リモート計算環境と共有したい Python ファイルを、このディレクトリに配置してください。次の2つのファイルを作成する必要があります：\n", "\n", "* `source_files\\diagonalization_engine.py`\n", "* `source_files\\classical_simulation.py`\n", "\n", "以下の各スクリプトのテキストをクリックして展開し、その内容を以下のパス名でローカルファイルにコピーして貼り付けてください。\n", "\n", "\n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " 詳細については、前述の「 [SQD IEF-PCM Qiskit 関数テンプレートガイド](/docs/guides/function-template-chemistry-workflow) 」（クリーブランド・クリニックと IBM が共同で開発）をご参照ください。 [`qiskit_addon_sqd`](/docs/addons/qiskit-addon-sqd) ライブラリも参照してください。\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7347c7e5", "metadata": { "tags": [ "id-diagonalization" ] }, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "import numpy as np\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "from functools import partial\n", "import os\n", "\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "from qiskit_serverless import distribute_task, get_arguments, get, save_result\n", "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", " solve_sci,\n", ")\n", "\n", "\n", "### Argument retrieval\n", "args = get_arguments()\n", "\n", "data = args[\"data\"] # Chemistry Data\n", "energy_tol = args[\"energy_tol\"] # SQD option\n", "occupancies_tol = args[\"occupancies_tol\"] # SQD option\n", "max_iterations = args[\"max_iterations\"] # SQD option\n", "symmetrize_spin = args[\"symmetrize_spin\"] # Eigenstate solver option\n", "carryover_threshold = args[\"carryover_threshold\"] # Eigenstate solver option\n", "num_batches = args[\"num_batches\"] # Eigenstate solver option\n", "samples_per_batch = args[\"samples_per_batch\"] # Eigenstate solver option\n", "max_cycle = args[\"max_cycle\"] # Eigenstate solver option\n", "mem = args[\"mem\"] # Memory per Worker\n", "\n", "\n", "# --- fan‑out target: 1 CPU + mem GB RAM per call -------------\n", "@distribute_task(target={\"cpu\": 1, \"mem\": mem * 1024**3})\n", "def _solve_sci_worker(\n", " ix, ci_strs, one_body_tensor, two_body_tensor, norb, nelec, spin_sq\n", "):\n", " print(f\">>>>> WORKER {ix} INITIATED\")\n", " res = solve_sci(\n", " ci_strs,\n", " one_body_tensor,\n", " two_body_tensor,\n", " norb=norb,\n", " nelec=nelec,\n", " spin_sq=spin_sq,\n", " )\n", "\n", " print(f\">>>>> WORKER {ix} COMPLETE\")\n", " return res\n", "\n", "\n", "def distribute_solve_sci_batch(\n", " ci_strings: list[tuple[np.ndarray, np.ndarray]],\n", " one_body_tensor: np.ndarray,\n", " two_body_tensor: np.ndarray,\n", " norb: int,\n", " nelec: tuple[int, int],\n", " *,\n", " spin_sq: float | None = None,\n", " **kwargs,\n", ") -> list[SCIResult]:\n", " \"\"\"Diagonalize Hamiltonian in subspaces, parallelizing across\n", " vCPUs in the Serverless environment.\n", "\n", " Args:\n", " ci_strings: List of pairs (strings_a, strings_b) of arrays of\n", " spin-alpha CI strings and spin-beta CI strings whose Cartesian\n", " product gives the basis of the subspace in which to perform a\n", " diagonalization.\n", " one_body_tensor: The one-body tensor of the Hamiltonian.\n", " two_body_tensor: The two-body tensor of the Hamiltonian.\n", " norb: The number of spatial orbitals.\n", " nelec: The numbers of alpha and beta electrons.\n", " spin_sq: Target value for the total spin squared for the ground state.\n", " If ``None``, no spin will be imposed.\n", " **kwargs: Keyword arguments to pass to\n", " `pyscf.fci.selected_ci.kernel_fixed_space`\n", " (https://pyscf.org/pyscf_api_docs/pyscf.fci.html#pyscf.fci.selected_ci.kernel_fixed_space\n", "\n", " Returns:\n", " The results of the diagonalizations in the subspaces given by ci_strings.\n", " \"\"\"\n", " inputs = [\n", " (ix, ci_strs, one_body_tensor, two_body_tensor, norb, nelec, spin_sq)\n", " for ix, ci_strs in enumerate(ci_strings)\n", " ]\n", "\n", " # fan‑out: spawn one worker per input tuple\n", " print(\">>>>> ENTERING WORKER FAN-OUT\")\n", " refs = [_solve_sci_worker(*input_) for input_ in inputs]\n", " print(\">>>>> WAITING ON WORKERS TO FINISH TASKS\")\n", "\n", " # fan‑in: block until every worker finishes\n", " results = get(refs)\n", " print(\">>>>> DISTRIBUTED JOBS COMPLETED\")\n", "\n", " return results\n", "\n", "\n", "# A caveat of executing a Python program remotely is\n", "# that the inputs to the remote program must be passed\n", "# over an internet network. Similarly, the outputs\n", "# must be passed back to the local program via the same\n", "# structure. Python objects are not always able to be\n", "# passed over a network, and must be encoded in a\n", "# JSON serializable format.\n", "i_data = JSONDecoder().decode(data)\n", "\n", "# i_data has all of the information needed from the\n", "# local program to pick up where the computation left off\n", "# after its submission to the remote environment.\n", "[\n", " job_id,\n", " hcore,\n", " eri,\n", " num_orbitals,\n", " nuclear_repulsion_energy,\n", " num_elec_a,\n", " num_elec_b,\n", "] = i_data\n", "\n", "# Re-convert data back into numpy format, after serialization\n", "hcore = np.array(hcore)\n", "eri = np.array(eri)\n", "nuclear_repulsion_energy = np.float64(nuclear_repulsion_energy)\n", "\n", "# Instantiate Runtime Service to retrieve the\n", "# bitstrings from the QPU job. We provided these\n", "# credentials upon Serverless setup.\n", "service = QiskitRuntimeService(\n", " channel=os.environ.get(\"QISKIT_IBM_CHANNEL\"),\n", " token=os.environ.get(\"QISKIT_IBM_TOKEN\"),\n", " instance=os.environ.get(\"QISKIT_IBM_INSTANCE\"),\n", ")\n", "\n", "# retrieving the QPU job data from the Serverless side\n", "job = service.job(job_id)\n", "primitive_result = job.result()\n", "pub_result = primitive_result[0]\n", "bit_array = pub_result.data.meas # Getting the bitstrings\n", "\n", "# Pass options to the built-in eigensolver\n", "sci_solver = partial(\n", " distribute_solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle\n", ")\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\">>>>> SQD ITERATION {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\">>>>> SUBSAMPLE {i}\")\n", " print(f\">>>>> \\tENERGY: {result.energy + nuclear_repulsion_energy}\")\n", " print(\n", " f\">>>>> \\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=num_orbitals,\n", " nelec=(num_elec_a, num_elec_b),\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", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=12345,\n", ")\n", "\n", "print(\">>>>> EXACT DIAGONALIZATION COMPLETE. CLEANING UP, SERIALIZING DATA.\")\n", "# Numpy arrays are not JSON serializable.\n", "# Convert them to List objects before using the JSONEncoder\n", "o_data = JSONEncoder().encode(\n", " [\n", " result.energy + nuclear_repulsion_energy,\n", " result.energy,\n", " result.rdm1.tolist(),\n", " result.rdm2.tolist(),\n", " [x.tolist() for x in result.orbital_occupancies],\n", " [\n", " result.sci_state.nelec,\n", " result.sci_state.norb,\n", " [x.tolist() for x in result.sci_state.orbital_occupancies()],\n", " [x.tolist() for x in result.sci_state.rdm()],\n", " ],\n", " ]\n", ")\n", "\n", "# JSON-safe package\n", "save_result({\"outputs\": o_data}) # single JSON blob returned to client" ] }, { "cell_type": "code", "execution_count": null, "id": "41807b3a", "metadata": { "tags": [ "id-classical" ] }, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "\n", "from qiskit_serverless import get_arguments, save_result\n", "\n", "import pyscf\n", "from pyscf import gto, scf\n", "from pyscf.solvent import pcm\n", "from pyscf.mcscf import avas\n", "\n", "import psutil\n", "\n", "mem_info = (\n", " psutil.virtual_memory()\n", ") # Get information about virtual memory (RAM)\n", "total_ram_gb = mem_info.total / (1024**3) # Convert bytes to GB\n", "print(f\">>>>> SERVERLESS TOTAL RAM: {total_ram_gb:.2f} GB\")\n", "\n", "### Argument retrieval\n", "args = get_arguments()\n", "data = args[\"data\"] # Chemistry Data\n", "\n", "i_data = JSONDecoder().decode(data)\n", "[mol_geo, eps, ao_labels] = i_data\n", "\n", "print(\">>>>> DEFINING MOLECULE\")\n", "mol = gto.M()\n", "mol.atom = mol_geo\n", "mol.basis = \"cc-pVDZ\"\n", "mol.unit = \"Ang\"\n", "mol.charge = 0\n", "mol.spin = 0\n", "mol.verbose = 0\n", "\n", "print(\">>>>> BUILDING MOLECULE\")\n", "mol.build()\n", "\n", "print(\">>>>> DEFINING PCM\")\n", "cm = pcm.PCM(mol)\n", "cm.eps = eps # for water\n", "cm.method = \"IEF-PCM\"\n", "\n", "print(\">>>>> BUILDING RESTRICTED HARTREE FOCK\")\n", "mf = scf.RHF(mol).PCM(cm) # This is the Final SCF object\n", "mf.kernel(verbose=0)\n", "\n", "print(\">>>>> RUNNING AVAS\")\n", "avas_ = avas.AVAS(mf, ao_labels, with_iao=True, canonicalize=True, verbose=0)\n", "avas_.kernel()\n", "norb, ne_act, mo_avas = avas_.ncas, avas_.nelecas, avas_.mo_coeff\n", "\n", "print(\">>>>> STARTING CASCI\")\n", "mc_pcm = pyscf.mcscf.CASCI(mf, norb, ne_act).PCM(\n", " cm\n", ") # Make sure to decorate the CASCI object with PCM\n", "mc_pcm.mo_coeff = mo_avas\n", "# mc_pcm.max_memory = 140000\n", "\n", "(CASCI_E, _, _, _, _) = mc_pcm.kernel(verbose=0)\n", "\n", "print(f\">>>>> CASCI_E: {CASCI_E}\")\n", "o_data = JSONEncoder().encode([float(CASCI_E)])\n", "\n", "# JSON-safe package\n", "save_result({\"outputs\": o_data}) # single JSON blob returned to client" ] }, { "cell_type": "markdown", "id": "c94aebd8", "metadata": {}, "source": [ "クラウド環境で実行するプログラムを共有し、ソースコードを変更するたびに再アップロードする必要があります：\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cc00669f", "metadata": {}, "outputs": [], "source": [ "client.upload(\n", " QiskitFunction(\n", " title=\"diagonalization_engine\",\n", " entrypoint=\"diagonalization_engine.py\", # lives in ./source_files\n", " working_dir=\"source_files\",\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "ef726eee", "metadata": {}, "outputs": [], "source": [ "client.upload(\n", " QiskitFunction(\n", " title=\"classical_simulation\",\n", " entrypoint=\"classical_simulation.py\", # lives in ./source_files\n", " working_dir=\"source_files\",\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "0cd9185e", "metadata": {}, "source": [ "\n", "\n", "## 小規模シミュレータの例\n", "\n" ] }, { "cell_type": "markdown", "id": "b0db693d", "metadata": {}, "source": [ "このチュートリアルでは、シミュレータによる探索の範囲を超えた、拡張性のある量子アプリケーションを実証することを目的としているため、小規模なシミュレータは使用しません。その代わりに、後ほど、CASCIと呼ばれる従来の最先端の比較手法を用いて、この手法をどのように実装できるかを示す。\n", "\n" ] }, { "cell_type": "markdown", "id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8", "metadata": {}, "source": [ "\n", "\n", "## 大規模なハードウェアの例\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5df67bda", "metadata": {}, "outputs": [], "source": [ "# This is a useful helper function that displays\n", "# remote job execution details to the user's local machine\n", "def feedback_serverless(serverless_job):\n", " import time\n", "\n", " # Wait for the job to execute\n", " print(f\">>>>> Serverless status: {serverless_job.job_id}\")\n", " timer = 0\n", " while timer < 10000:\n", " if (\n", " serverless_job.status() == \"QUEUED\"\n", " or serverless_job.status() == \"INITIALIZING\"\n", " or serverless_job.status() == \"RUNNING\"\n", " ):\n", " print(f\">>>>> [{timer}s] Serverless job {serverless_job.job_id}: \\\n", " {serverless_job.status()}\")\n", " time.sleep(10)\n", " timer += 10\n", "\n", " elif serverless_job.status() == \"ERROR\":\n", " print(\n", " f\">>>>> Serverless job {serverless_job.job_id}: {serverless_job.status()}\"\n", " )\n", " print(\">>>>> Logs:\")\n", " print(serverless_job.logs())\n", " break\n", "\n", " elif serverless_job.status() == \"DONE\":\n", " print(\n", " f\">>>>> Serverless job {serverless_job.job_id}: {serverless_job.status()}\"\n", " )\n", " break\n", "\n", " else:\n", " break\n", "\n", " return" ] }, { "attachments": {}, "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "\n", "\n", "## ステップ1：古典的な入力を量子問題に写像する\n", "\n" ] }, { "cell_type": "markdown", "id": "cf4bbd59", "metadata": {}, "source": [ "\n", "\n", "### 1.1: 既知の $\\textit{a priori}$ 分子構造を用いて分子オブジェクトを初期化する\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "566e06b4", "metadata": {}, "outputs": [], "source": [ "# Reference guide for building molecule structures:\n", "# https://pyscf.org/user/gto.html\n", "# Video tutorial on building molecular objects in PySCF:\n", "# https://www.youtube.com/watch?v=cNC2cY9E9j0\n", "\n", "molecule_name = \"Methylamine\"\n", "\n", "methylamine_geo = \"\"\"\n", " N -0.7154 0.0000 0.0000;\n", " C 0.7154 0.0000 0.0000;\n", " H 1.1069 0.0916 1.0174;\n", " H 1.0996 0.8349 -0.5930;\n", " H 1.0996 -0.9274 -0.4345;\n", " H -1.0625 0.8564 0.4294;\n", " H -1.0625 -0.7661 0.5753;\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": null, "id": "4a7aec02", "metadata": {}, "outputs": [], "source": [ "# Imports\n", "import pyscf\n", "from pyscf import gto # Deals with molecular initialization\n", "from pyscf import scf # Solvation methods\n", "\n", "# Explicitly defining the Methylamine molecule\n", "mol = gto.M()\n", "mol.atom = methylamine_geo\n", "mol.basis = \"cc-pVDZ\"\n", "mol.unit = \"Ang\"\n", "mol.charge = 0\n", "mol.spin = 0\n", "mol.verbose = 0\n", "\n", "mol.build()" ] }, { "cell_type": "markdown", "id": "99359d80", "metadata": {}, "source": [ "\n", "\n", "### 1.2: 極化連続体モデル（PCM）を用いて溶媒和効果を定義する\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9294b528", "metadata": {}, "outputs": [], "source": [ "# You can explore other solvents (such as methanol) by\n", "# retrieving other dielectric parameters from:\n", "# https://gaussian.com/scrf/\n", "from pyscf.solvent import pcm\n", "\n", "eps_water = 78.3553 # If solvating in a different medium,\n", "# set this constant appropriately using a known value" ] }, { "cell_type": "code", "execution_count": null, "id": "ffccfb92", "metadata": {}, "outputs": [], "source": [ "cm = pcm.PCM(mol)\n", "cm.eps = eps_water # PySCF defaults to water solvation,\n", "# but here we show this solvation parameter explicitly\n", "\n", "cm.method = (\n", " \"IEF-PCM\" # Alternative solvation models include C-PCM, SS(V)PE, COSMO\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "bfef3c9b", "metadata": {}, "outputs": [], "source": [ "# Create a \"Restricted Hartree-Fock\" object for the solute,\n", "# then wrap the SCF object with a Polarizable Continuum Model\n", "mf_pcm0 = scf.RHF(mol).PCM(\n", " cm\n", ") # Restricted Hartree-Fock misses instantaneous correlations,\n", "# post-HF methods like CCSD, CI, MP2 might be worth exploring" ] }, { "cell_type": "markdown", "id": "b5e6c20a", "metadata": {}, "source": [ "\n", "\n", "### 1.3: TRICを用いた幾何最適化\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "493402b3", "metadata": {}, "outputs": [], "source": [ "# Geometry optimization with geomeTRIC\n", "from pyscf.geomopt.geometric_solver import (\n", " optimize,\n", ") # GeomeTRIC under the hood, for geometry optimization\n", "\n", "mol_opt = optimize(\n", " mf_pcm0, tol_grad=3e-4, verbose=0\n", ") # Use geomeTRIC/TRIC under the hood" ] }, { "cell_type": "markdown", "id": "58c168bf", "metadata": {}, "source": [ "\n", "\n", "### 1.4: 関連する変数を含む連続体モデルと平均場オブジェクトを準備する\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1834cb22", "metadata": {}, "outputs": [], "source": [ "from pyscf.mcscf import avas\n", "\n", "# Re-define PCM\n", "cm = pcm.PCM(mol_opt)\n", "cm.eps = eps_water # for water\n", "cm.method = \"IEF-PCM\"\n", "\n", "# Re-build Restricted Hartree Fock object\n", "mf_opt = scf.RHF(mol_opt).PCM(cm)\n", "mf_opt.kernel(verbose=0)\n", "\n", "# Run AVAS\n", "ao_labels = [\"C 2s\", \"C 2p\", \"N 2s\", \"N 2p\", \"H 1s\"]\n", "avas_ = avas.AVAS(\n", " mf_opt, ao_labels, with_iao=True, canonicalize=True, verbose=0\n", ")\n", "avas_.kernel()\n", "norb, ne_act, mo_avas = avas_.ncas, avas_.nelecas, avas_.mo_coeff\n", "\n", "num_elec_a = (ne_act + mol_opt.spin) // 2\n", "num_elec_b = (ne_act - mol_opt.spin) // 2" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "\n", "\n", "## ステップ2：量子ハードウェアでの実行に向けて問題を最適化する\n", "\n" ] }, { "cell_type": "markdown", "id": "52506835", "metadata": {}, "source": [ "ここで紹介したヘルパー関数の詳細については、 [「化学ハミルトニアンのサンプルベース量子対角化」](/docs/tutorials/sample-based-quantum-diagonalization) チュートリアルを参照してください。\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dc27c42e", "metadata": {}, "outputs": [], "source": [ "# Standard SQD helper functions (From SQD Tutorial)\n", "\n", "from typing import Sequence\n", "import rustworkx\n", "from qiskit.providers import BackendV2\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "\n", "from rustworkx import NoEdgeBetweenNodes, PyGraph\n", "\n", "IBM_TWO_Q_GATES = {\"cx\", \"ecr\", \"cz\"}\n", "\n", "\n", "def create_linear_chains(num_orbitals: int) -> PyGraph:\n", " \"\"\"In zig-zag layout, there are two linear chains (with connecting\n", " qubits between the chains). This function creates those two linear\n", " chains: a rustworkx PyGraph with two disconnected linear chains.\n", " Each chain contains `num_orbitals` number of nodes, that is, in the\n", " final graph there are `2 * num_orbitals` number of nodes.\n", "\n", " Args:\n", " num_orbitals (int): Number orbitals or nodes in each linear chain.\n", " They are also known as alpha-alpha interaction qubits.\n", "\n", " Returns:\n", " A rustworkx.PyGraph with two disconnected linear chains each with\n", " `num_orbitals` number of nodes.\n", " \"\"\"\n", " G = rustworkx.PyGraph()\n", "\n", " for n in range(num_orbitals):\n", " G.add_node(n)\n", "\n", " for n in range(num_orbitals - 1):\n", " G.add_edge(n, n + 1, None)\n", "\n", " for n in range(num_orbitals, 2 * num_orbitals):\n", " G.add_node(n)\n", "\n", " for n in range(num_orbitals, 2 * num_orbitals - 1):\n", " G.add_edge(n, n + 1, None)\n", "\n", " return G\n", "\n", "\n", "def create_lucj_zigzag_layout(\n", " num_orbitals: int, backend_coupling_graph: PyGraph\n", ") -> tuple[PyGraph, int]:\n", " \"\"\"This function creates the complete zigzag graph that 'can be mapped'\n", " to an IBM QPU with heavy-hex connectivity (the zigzag must be an\n", " isomorphic sub-graph to the QPU/backend coupling graph for it to be\n", " mapped). The zigzag pattern includes both linear chains (alpha-alpha\n", " interactions) and connecting qubits between the linear chains\n", " (alpha-beta interactions).\n", "\n", " Args:\n", " num_orbitals (int): Number of orbitals, that is, number of nodes in\n", " each alpha-alpha linear chain.\n", " backend_coupling_graph (PyGraph): The coupling graph of the backend\n", " on which the LUCJ ansatz will be mapped and run. This function takes\n", " the coupling graph as a undirected `rustworkx.PyGraph` where there\n", " is only one 'undirected' edge between two nodes, that is, qubits.\n", " Usually, the coupling graph of an IBM backend is directed (for\n", " example, Eagle devices such as ibm_brisbane) or may have two edges\n", " between two nodes (for example, Heron `ibm_torino`). A user\n", " needs to make such graphs undirected or remove duplicate edges\n", " (or do both) to make them compatible with this function.\n", "\n", " Returns:\n", " G_new (PyGraph): The graph with IBM backend compliant zigzag pattern.\n", " num_alpha_beta_qubits (int): Number of connecting qubits between the\n", " linear chains in the zigzag pattern. While we want as many\n", " connecting (alpha-beta) qubits between the linear (alpha-alpha)\n", " chains, we cannot accommodate all due to qubit and connectivity\n", " constraints of backends. This is the maximum number of connecting\n", " qubits the zigzag pattern can have while being backend compliant\n", " (that is, isomorphic to backend coupling graph).\n", " \"\"\"\n", " isomorphic = False\n", " G = create_linear_chains(num_orbitals=num_orbitals)\n", "\n", " num_iters = num_orbitals\n", " while not isomorphic:\n", " G_new = G.copy()\n", " num_alpha_beta_qubits = 0\n", " for n in range(num_iters):\n", " if n % 4 == 0:\n", " new_node = 2 * num_orbitals + num_alpha_beta_qubits\n", " G_new.add_node(new_node)\n", " G_new.add_edge(n, new_node, None)\n", " G_new.add_edge(new_node, n + num_orbitals, None)\n", " num_alpha_beta_qubits = num_alpha_beta_qubits + 1\n", " isomorphic = rustworkx.is_subgraph_isomorphic(\n", " backend_coupling_graph, G_new\n", " )\n", " num_iters -= 1\n", "\n", " return G_new, num_alpha_beta_qubits\n", "\n", "\n", "def lightweight_layout_error_scoring(\n", " backend: BackendV2,\n", " virtual_edges: Sequence[Sequence[int]],\n", " physical_layouts: Sequence[int],\n", " two_q_gate_name: str,\n", ") -> list[list[list[int], float]]:\n", " \"\"\"Lightweight and heuristic function to score isomorphic layouts. There\n", " can be many zigzag patterns, each with different set of physical qubits,\n", " that can be mapped to a backend. Some of them might include fewer noise\n", " qubits and couplings than others. This function computes a simple error\n", " score for each such layout. It sums up 2Q gate error for all couplings\n", " in the zigzag pattern (layout) and measurement of errors of physical\n", " qubits in the layout to compute the error score.\n", "\n", " Note:\n", " This lightweight scoring can be refined using concepts such as\n", " mapomatic.\n", "\n", " Args:\n", " backend (BackendV2): A backend.\n", " virtual_edges (Sequence[Sequence[int]]): Edges in the device-\n", " compliant zigzag pattern where nodes are numbered from 0 to (2 *\n", " num_orbitals + num_alpha_beta_qubits).\n", " physical_layouts (Sequence[int]): All physical layouts of the zigzag\n", " pattern that are isomorphic to each other and to the larger backend\n", " coupling map.\n", " two_q_gate_name (str): The name of the two-qubit gate of the\n", " backend. The name is used for fetching two-qubit gate error from\n", " backend properties.\n", "\n", " Returns:\n", " scores (list): A list of lists where each sublist contains two\n", " items. First item is the layout, and second item is a float\n", " representing error score of the layout. The layouts in the `scores`\n", " are sorted in the ascending order of error score.\n", " \"\"\"\n", " props = backend.properties()\n", " scores = []\n", " for layout in physical_layouts:\n", " total_2q_error = 0\n", " for edge in virtual_edges:\n", " physical_edge = (layout[edge[0]], layout[edge[1]])\n", " try:\n", " ge = props.gate_error(two_q_gate_name, physical_edge)\n", " except Exception:\n", " ge = props.gate_error(two_q_gate_name, physical_edge[::-1])\n", " total_2q_error += ge\n", " total_measurement_error = 0\n", " for qubit in layout:\n", " meas_error = props.readout_error(qubit)\n", " total_measurement_error += meas_error\n", " scores.append([layout, total_2q_error + total_measurement_error])\n", " return sorted(scores, key=lambda x: x[1])\n", "\n", "\n", "def _make_backend_cmap_pygraph(backend: BackendV2) -> PyGraph:\n", " graph = backend.coupling_map.graph\n", " if not graph.is_symmetric():\n", " graph.make_symmetric()\n", " backend_coupling_graph = graph.to_undirected()\n", "\n", " edge_list = backend_coupling_graph.edge_list()\n", " removed_edge = []\n", " for edge in edge_list:\n", " if set(edge) in removed_edge:\n", " continue\n", " try:\n", " backend_coupling_graph.remove_edge(edge[0], edge[1])\n", " removed_edge.append(set(edge))\n", " except NoEdgeBetweenNodes:\n", " pass\n", "\n", " return backend_coupling_graph\n", "\n", "\n", "def get_zigzag_physical_layout(\n", " num_orbitals: int, backend: BackendV2, score_layouts: bool = True\n", ") -> tuple[list[int], int]:\n", " \"\"\"The main function that generates the zigzag pattern\n", " with physical qubits that can be used as an `intial_layout` in a\n", " preset passmanager/transpiler.\n", "\n", " Args:\n", " num_orbitals (int): Number of orbitals.\n", " backend (BackendV2): A backend.\n", " score_layouts (bool): Optional. If `True`, it uses the\n", " `lightweight_layout_error_scoring` function to score the\n", " isomorphic layouts and returns the layout with\n", " fewer erroneous qubits.\n", " If `False`, returns the first isomorphic subgraph.\n", "\n", " Returns:\n", " A tuple of device compliant layout (list[int]) with zigzag pattern\n", " and an int representing number of alpha-beta-interactions.\n", " \"\"\"\n", " backend_coupling_graph = _make_backend_cmap_pygraph(backend=backend)\n", "\n", " G, num_alpha_beta_qubits = create_lucj_zigzag_layout(\n", " num_orbitals=num_orbitals,\n", " backend_coupling_graph=backend_coupling_graph,\n", " )\n", "\n", " isomorphic_mappings = rustworkx.vf2_mapping(\n", " backend_coupling_graph, G, subgraph=True\n", " )\n", " isomorphic_mappings = list(isomorphic_mappings)\n", "\n", " edges = list(G.edge_list())\n", "\n", " layouts = []\n", " for mapping in isomorphic_mappings:\n", " initial_layout = [None] * (2 * num_orbitals + num_alpha_beta_qubits)\n", " for key, value in mapping.items():\n", " initial_layout[value] = key\n", " layouts.append(initial_layout)\n", "\n", " two_q_gate_name = IBM_TWO_Q_GATES.intersection(\n", " backend.configuration().basis_gates\n", " ).pop()\n", "\n", " if score_layouts:\n", " scores = lightweight_layout_error_scoring(\n", " backend=backend,\n", " virtual_edges=edges,\n", " physical_layouts=layouts,\n", " two_q_gate_name=two_q_gate_name,\n", " )\n", "\n", " return scores[0][0][:-num_alpha_beta_qubits], num_alpha_beta_qubits\n", "\n", " return layouts[0][:-num_alpha_beta_qubits], num_alpha_beta_qubits" ] }, { "cell_type": "code", "execution_count": null, "id": "87f8e3ec", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler import generate_preset_pass_manager\n", "import ffsim\n", "\n", "# Initial LUCJ ansatz layout\n", "initial_layout, _ = get_zigzag_physical_layout(norb, backend=backend)\n", "\n", "# Initialize a pass manager\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend, initial_layout=initial_layout\n", ")\n", "\n", "pass_manager.pre_init = ffsim.qiskit.PRE_INIT" ] }, { "attachments": {}, "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "\n", "\n", "## 手順 3 および 4：Qiskit を使用して実行し、 Qiskit Serverless を使用して後処理を行う\n", "\n", "ここでは、ステップ3（実行）とステップ4（後処理）を統合しています。これは、暗黙的溶媒モデルの適用において、最終的な計算結果を改善するために、実行と後処理を数サイクル繰り返す反復的なフィードバックループが必要となるためです。\n", "\n" ] }, { "cell_type": "markdown", "id": "8715bb6b", "metadata": {}, "source": [ "\n", "\n", "### 3.1 制限付きハートリー・フォックエネルギーを計算する\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3aa879ad", "metadata": {}, "outputs": [], "source": [ "# Run the kernel to get the RHF energy\n", "mf_opt = scf.RHF(mol_opt).PCM(cm)\n", "hf_e = float(mf_opt.kernel())" ] }, { "cell_type": "code", "execution_count": null, "id": "480dd476", "metadata": {}, "outputs": [], "source": [ "print(f\"Restricted Hartree-Fock Energy: {hf_e}\")" ] }, { "cell_type": "markdown", "id": "fd0145c6", "metadata": {}, "source": [ "\n", "\n", "### 3.2: CASCIを用いて古典的な基準エネルギーを決定する\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f629119d", "metadata": {}, "outputs": [], "source": [ "# Setup a Serverless Client\n", "worker = client.load(\"classical_simulation\")" ] }, { "cell_type": "code", "execution_count": null, "id": "bb9d2407", "metadata": {}, "outputs": [], "source": [ "from json.encoder import JSONEncoder\n", "\n", "ao_labels = [\"C 2s\", \"C 2p\", \"N 2s\", \"N 2p\", \"H 1s\"]\n", "data_e = JSONEncoder().encode([mol_opt.tostring(), eps_water, ao_labels])" ] }, { "cell_type": "code", "execution_count": null, "id": "93874c58", "metadata": {}, "outputs": [], "source": [ "serverless_job = worker.run(data=data_e)" ] }, { "cell_type": "code", "execution_count": null, "id": "756f0fdf", "metadata": {}, "outputs": [], "source": [ "# Optionally, check the Serverless status feedback\n", "# Don't sit here and stare at the feedback unless debugging.\n", "# You can go develop something else while the Serverless job runs.\n", "_ = feedback_serverless(serverless_job)" ] }, { "cell_type": "code", "execution_count": null, "id": "d3ede43c", "metadata": {}, "outputs": [], "source": [ "# If you make a mistake and need to cancel something\n", "\n", "# for job in client.jobs():\n", "# job.cancel()" ] }, { "cell_type": "code", "execution_count": null, "id": "7379c49f", "metadata": {}, "outputs": [], "source": [ "from json.decoder import JSONDecoder\n", "\n", "CASCI_E = JSONDecoder().decode(serverless_job.result()[\"outputs\"])[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "0a5f4b88", "metadata": {}, "outputs": [], "source": [ "# We have approximated the red, classical baseline from\n", "# Figure 1 for Methanol (North-West panel)\n", "print(f\"CASCI/IEF-PCM(cc-pVDZ): E={CASCI_E}\")" ] }, { "cell_type": "markdown", "id": "b22a1b00", "metadata": {}, "source": [ "\n", "\n", "## アプリケーションパラメータの設定\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "618375af", "metadata": {}, "outputs": [], "source": [ "# Systematically vary these parameters to improve hardware results\n", "\n", "# Set to \"True\" to run on real hardware\n", "use_hardware = True\n", "\n", "# Error suppression/mitigation options\n", "# >> Configure within Sampler primitive\n", "\n", "# Transpiler Options\n", "optimization_level = 3\n", "\n", "# Heartwood algorithm options\n", "n_iter = 15 # How many update loops to run\n", "resample = 1 # (resample=1 -> resample the QPU after every\n", "# update loop; resample=n_iter -> sample QPU only once)\n", "shots = 10000\n", "\n", "# SQD options\n", "energy_tol = 1e-4\n", "occupancies_tol = 1e-3\n", "max_iterations = 12\n", "\n", "# Eigenstate solver options\n", "num_batches = 5\n", "samples_per_batch = 300\n", "symmetrize_spin = True\n", "carryover_threshold = 1e-5\n", "max_cycle = 200\n", "\n", "# Classical post-processing options\n", "local = (\n", " False # Remote, Serverless (False) versus Local Post-Processing (True)\n", ")\n", "mem = 16 # Memory allocated to each diagonalization worker (Gb)" ] }, { "cell_type": "code", "execution_count": null, "id": "ce75a126", "metadata": {}, "outputs": [], "source": [ "# Heartwood algorithm subroutines\n", "import os\n", "import numpy as np\n", "import pyscf\n", "from pyscf import ao2mo, cc\n", "from functools import reduce\n", "import ffsim\n", "\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "import time\n", "\n", "from functools import partial\n", "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", " solve_sci_batch,\n", ")\n", "\n", "\n", "def update_rdm(casci_object, dmas):\n", " \"\"\"\n", " Inputs:\n", " mc -> CASCI object\n", " dmas -> Spin-summed 1-particle reduced density matrix\n", "\n", " This function returns the CASCI/SQD one-body density matrix in\n", " the full basis of atomic orbitals, written as the sum (last line)\n", " of two terms:\n", " - a contribution from the core orbitals,\n", " np.dot(mocore, mocore.conj().T) * 2, (core = inactive and doubly-occupied)\n", " - a contribution from the active-space orbitals and electrons (dmas)\n", " rotated from the active-space to the AO basis (the reduce operation)\n", "\n", " Outputs:\n", " rho_approximation: The CASCI/SQD one-body density matrix\n", " in the full basis of atomic orbitals\n", "\n", " \"\"\"\n", " mo_coeff = casci_object.mo_coeff\n", " ncore = casci_object.ncore\n", " ncas = casci_object.ncas\n", " mocore = mo_coeff[:, :ncore]\n", " mocas = mo_coeff[:, ncore : ncore + ncas]\n", " dm1 = np.dot(mocore, mocore.conj().T) * 2\n", "\n", " rho_approximation = dm1 + reduce(np.dot, (mocas, dmas, mocas.conj().T))\n", " return rho_approximation\n", "\n", "\n", "def run_active_space_calculation(\n", " h1e_cas, h2e_cas, norb, ne_act, orbs, fermilevel, ecore\n", "):\n", " # ----- perform an HF and a CCSD calculation in the active space\n", " from pyscf import tools\n", " from datetime import datetime\n", "\n", " now = datetime.now().strftime(\"%H:%M:%S\")\n", " print(\">>>>> ACTIVE SPACE CALCULATIONS \")\n", " tools.fcidump.from_integrals(\n", " f\"as_fcidump_{now}.txt\",\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " ms=0,\n", " nuc=ecore,\n", " ) # Forcefully represents the active space in the correct structure\n", " mf_as = tools.fcidump.to_scf(f\"as_fcidump_{now}.txt\")\n", " os.remove(f\"as_fcidump_{now}.txt\")\n", " mf_as.kernel()\n", " print(\">>>>> RUNNING CCSD\")\n", "\n", " mf_cc = cc.CCSD(mf_as)\n", " mf_cc.kernel()\n", " orbts = mf_as.mo_coeff\n", " t1, t2 = mf_cc.t1, mf_cc.t2\n", " print(\">>>>> UPDATED t1, t2 PARAMETERS\")\n", "\n", " # ----- update the HF orbitals\n", " active = list(\n", " range(fermilevel - ne_act // 2, fermilevel - ne_act // 2 + norb)\n", " )\n", " orbs[:, active] = np.dot(orbs[:, active], orbts)\n", " return orbs, t1, t2\n", "\n", "\n", "def get_lucj(norb, num_elec_a, num_elec_b, t1, t2, n_reps=1):\n", " print(\">>>>> CONSTRUCTING LUCJ CIRCUIT\")\n", "\n", " alpha_alpha_indices = [(p, p + 1) for p in range(norb - 1)]\n", " alpha_beta_indices = [(p, p) for p in range(0, norb, 4)]\n", "\n", " ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n", " t1=t1, # <---- Update t1 each loop\n", " t2=t2, # <---- Update t2 each loop\n", " n_reps=n_reps,\n", " interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n", " )\n", " nelec = (num_elec_a, num_elec_b)\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\n", " # and append it 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", " return circuit\n", "\n", "\n", "# Classical diagonalization engine sent to HPC\n", "def classically_diagonalize(\n", " bit_array=None, # Bit string array (only needed if locally processing data)\n", " nuclear_repulsion_energy=None, # Electronic energy from the core orbitals\n", " hcore=None, # 1-electron hamiltonian integrals\n", " eri=None, # 2-electron hamiltonian integrals\n", " num_orbitals=None, # Number of spatial orbitals\n", " nelec=None, # Number of electrons\n", " num_elec_a=None, # Alpha orbitals\n", " num_elec_b=None, # Beta orbitals\n", " job_id=None, # QPU bitstring Job ID\n", " client=None, # Diagonalization engine worker\n", " energy_tol=1e-4, # SQD option\n", " occupancies_tol=1e-3, # SQD option\n", " max_iterations=12, # SQD option\n", " num_batches=8, # Eigenstate solver option\n", " samples_per_batch=300, # Eigenstate solver option\n", " symmetrize_spin=False, # Eigenstate solver option\n", " carryover_threshold=1e-5, # Eigenstate solver option\n", " max_cycle=200, # Eigenstate solver option\n", " local=True, # Remote vs Local Diagonalization\n", " mem=16.0, # Memory per Serverless Worker (Gb)\n", "):\n", " print(\">>>>> STARTING DIAGONALIZATION ENGINE \")\n", " # Pass options to the built-in eigensolver. If you just want to use\n", " # the defaults, you can omit this step, in which case you would not\n", " # specify the sci_solver argument in the call to\n", " # diagonalize_fermionic_hamiltonian below.\n", " if local:\n", " sci_solver = partial(\n", " solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle\n", " )\n", "\n", " # List to capture intermediate results\n", " result_history = []\n", "\n", " def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\">>>>> SQD ITERATION {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\">>>>> SUBSAMPLE {i}\")\n", " print(\n", " f\">>>>> \\tENERGY: {result.energy + nuclear_repulsion_energy}\"\n", " )\n", " print(\n", " f\">>>>> \\tSUBSPACE DIMENSION: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", " result = diagonalize_fermionic_hamiltonian(\n", " hcore,\n", " eri,\n", " bit_array,\n", " samples_per_batch=samples_per_batch,\n", " norb=num_orbitals,\n", " nelec=(nelec // 2, nelec // 2),\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", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=12345,\n", " )\n", "\n", " result = (result.energy, result.rdm1, result.rdm2)\n", "\n", " else:\n", " # Serverless Logic\n", " print(\n", " f\">>>>> SENDING QISKIT RUNTIME JOB {job_id} TO QISKIT SERVERLESS\"\n", " )\n", "\n", " data = [\n", " job_id,\n", " hcore.tolist(),\n", " eri.tolist(),\n", " int(num_orbitals),\n", " float(nuclear_repulsion_energy),\n", " int(num_elec_a),\n", " int(num_elec_b),\n", " ]\n", "\n", " # Encode the execution dependencies with the JSONEncoder\n", " data_e = JSONEncoder().encode(data)\n", "\n", " # Send to Serverless\n", " worker = client.load(\"diagonalization_engine\")\n", " serverless_job = worker.run(\n", " data=data_e,\n", " energy_tol=energy_tol, # SQD option\n", " occupancies_tol=occupancies_tol, # SQD option\n", " max_iterations=max_iterations, # SQD option\n", " symmetrize_spin=symmetrize_spin, # Eigenstate solver option\n", " carryover_threshold=carryover_threshold, # Eigenstate solver option\n", " num_batches=num_batches, # Eigenstate solver option\n", " samples_per_batch=samples_per_batch, # Eigenstate solver option\n", " max_cycle=max_cycle, # Eigenstate solver option\n", " mem=mem, # Memory per Worker (Gb)\n", " )\n", "\n", " # Wait for the job to execute\n", " _ = feedback_serverless(serverless_job)\n", "\n", " o_data = JSONDecoder().decode(serverless_job.result()[\"outputs\"])\n", " result = (o_data[1], np.array(o_data[2]), np.array(o_data[3]))\n", "\n", " print(f\">>>>>>>>>> Active Space Energy: {o_data[1]}\")\n", " print(f\">>>>>>>>>> rdm1: {o_data[2]}\")\n", " print(f\">>>>>>>>>> rdm2: {o_data[3]}\")\n", "\n", " return result" ] }, { "cell_type": "code", "execution_count": null, "id": "d1e23956", "metadata": {}, "outputs": [], "source": [ "# The Heartwood algorithm\n", "import numpy as np\n", "\n", "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n", "\n", "mc = pyscf.mcscf.CASCI(mf_opt, ncas=norb, nelecas=ne_act).PCM(cm)\n", "mc.with_solvent.method = mf_opt.with_solvent.method # Here we make sure\n", "# that mc is also using the same solvent method defined earlier (IEF-PCM)\n", "mc.with_solvent.eps = mf_opt.with_solvent.eps # Set the dielectric parameters\n", "mc.mo_coeff = mo_avas.copy() # Update the molecular orbitals to include\n", "# those computed in the presence of the solvent\n", "\n", "h1e_cas, ecore = (\n", " mc.get_h1eff()\n", ") # <-- h1eff is the 1-electron hamiltonian integrals. h1e_cas is\n", "# a common alias. ecore is the electronic energy from the core orbitals.\n", "h2e_cas = ao2mo.restore(\n", " 1, mc.get_h2eff(), norb\n", ") # <-- get the 2-electron hamiltonian integrals\n", "\n", "mc.mo_coeff, t1, t2 = run_active_space_calculation(\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " mo_avas.copy(),\n", " mf_opt.mol.nelectron // 2,\n", " ecore,\n", ")\n", "\n", "# Sampler primitive options\n", "sampler = Sampler(mode=backend)\n", "\n", "# Explore error suppression techniques and see if they can improve result quality\n", "sampler.options.dynamical_decoupling.enable = True\n", "sampler.options.dynamical_decoupling.sequence_type = \"XY4\"\n", "sampler.options.twirling.enable_measure = True\n", "sampler.options.environment.job_tags = [\"TUT_ISC\"]\n", "# sampler.options.twirling.enable_gates = False\n", "# sampler.options.twirling.num_randomizations = 10\n", "# sampler.options.twirling.shots_per_randomization = 1024\n", "\n", "# initial approximation for rdm1\n", "with_solvent_e, with_solvent_v = None, None # Don't touch\n", "data = []\n", "for iiter in range(n_iter):\n", " print(f\">>>>> IMPLICIT SOLVENT ITERATION {iiter+1}/{n_iter}\")\n", " if with_solvent_v is not None:\n", " # Subsequent update loops enter here\n", " mc.get_hcore = lambda *args: mc._scf.get_hcore() + with_solvent_v\n", " else:\n", " # First update loop starts here\n", " # hcore is the CAS space (classically computed) 1-electron\n", " # hamiltonian, which we default to at the start of the routine.\n", " mc.get_hcore = (\n", " lambda *args: mc._scf.get_hcore()\n", " ) # REF: https://pyscf.org/pyscf_api_docs/pyscf.mcscf.html#pyscf.scf.hf.CASBase.get_h1cas\n", "\n", " # Alias mapping\n", " # hcore : h1e_cas : h1e_eff\n", " # nuclear_repulsion_energy : ecore\n", " # eri : h2e_cas : h2e_eff\n", "\n", " h1e_cas, ecore = (\n", " mc.get_h1eff()\n", " ) # <-- h1eff is the 1-electron hamiltonian integrals. h1e_cas is a\n", " # common alias. ecore is the electronic energy from the core orbitals.\n", " h2e_cas = ao2mo.restore(\n", " 1, mc.get_h2eff(), norb\n", " ) # <-- get the 2-electron hamiltonian integrals\n", "\n", " mc.mo_coeff, t1, t2 = run_active_space_calculation(\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " mo_avas.copy(),\n", " mf_opt.mol.nelectron // 2,\n", " ecore,\n", " )\n", "\n", " if use_hardware:\n", " if (\n", " iiter % resample == 0\n", " ): # <-- Toggle how often you refresh your bitstrings here. The\n", " # developer suggests that you do it every time, but benevolently\n", " # provides the freedom to disagree with him via the resample\n", " # control variable.\n", " # The \"Quantum-Centric\" part\n", " print(\">>>>> GENERATING BITSTRINGS USING QUANTUM HARDWARE\")\n", " # LUCJ Ansatz construction\n", " circuit = get_lucj(norb, num_elec_a, num_elec_b, t1, t2, n_reps=1)\n", "\n", " print(f\">>>>> TRANSPILING LUCJ TO {backend.name}\")\n", " isa_circuit = pass_manager.run(circuit)\n", " print(f\">>>>> SUBMITTING ISA_CIRCUIT TO {backend.name}\")\n", " job = sampler.run(\n", " [isa_circuit], shots=shots\n", " ) # <----- Error Suppression/Mitigation configured performed upstream\n", " job_id = str(job.job_id())\n", "\n", " timer = 0\n", " while job.status() != \"DONE\":\n", " timer += 10\n", " print(\n", " f\">>>>> [{timer}s] RUNTIME JOB {job_id}: {job.status()}\"\n", " )\n", " time.sleep(10)\n", "\n", " primitive_result = job.result()\n", " print(f\">>>>> RETRIEVED {job_id} FROM {backend.name}\")\n", "\n", " pub_result = primitive_result[0]\n", " bit_array = pub_result.data.meas\n", "\n", " else:\n", " print(\">>>>> GENERATING BITSTRINGS CLASSICALLY\")\n", " rng = np.random.default_rng(24)\n", " bit_array = generate_bit_array_uniform(\n", " 100_000, 2 * norb, rand_seed=rng\n", " ) # <-- Sample bitstrings from a uniform distribution. This is\n", " # useful for debug, but runs out of steam on large systems\n", " job_id = float(\n", " \"nan\"\n", " ) # <-- we will check that valid job_id's were passed during grading\n", " local = True\n", "\n", " # The \"Classical Post-processing\" part\n", " result = classically_diagonalize(\n", " bit_array=bit_array,\n", " nuclear_repulsion_energy=ecore, # Electronic energy from the core orbitals\n", " hcore=h1e_cas, # 1-electron hamiltonian integrals\n", " eri=h2e_cas, # 2-electron hamiltonian integrals\n", " num_orbitals=norb, # Number of spatial orbitals\n", " nelec=ne_act, # Number of electrons\n", " num_elec_a=ne_act // 2, # Alpha orbitals\n", " num_elec_b=ne_act // 2, # Beta orbitals\n", " job_id=job_id, # QPU bitstring Job ID\n", " client=client, # Diagonalization engine worker\n", " energy_tol=energy_tol, # SQD option\n", " occupancies_tol=occupancies_tol, # SQD option\n", " max_iterations=max_iterations, # SQD option\n", " num_batches=num_batches, # Eigenstate solver option\n", " samples_per_batch=samples_per_batch, # Eigenstate solver option\n", " symmetrize_spin=symmetrize_spin, # Eigenstate solver option\n", " carryover_threshold=carryover_threshold, # Eigenstate solver option\n", " max_cycle=max_cycle, # Eigenstate solver option\n", " local=local, # Remote vs Local Diagonalization\n", " mem=mem, # Memory per Worker (Gb)\n", " )\n", "\n", " # e : SQD-based estimate of the energy\n", " # rdm1: Spin-summed 1-particle reduced density matrix\n", " e, rdm1 = result[0], result[1]\n", " rho_approximation = update_rdm(\n", " mc, rdm1\n", " ).copy() # <--- Reconstruct the one-body density matrix in the\n", " # atomic orbital basis to update the external potential due to\n", " # the solvent\n", "\n", " if with_solvent_e is not None:\n", " # Subsequent update loops enter here\n", " edup = np.einsum(\n", " \"ij,ji->\", with_solvent_v, rho_approximation\n", " ) # <-- edup: Incrementing the energy calculation with\n", " # subsequent iterations\n", " e += ecore + with_solvent_e - edup\n", "\n", " else:\n", " # First update loop enters here\n", " e += (\n", " ecore # Pulled from the CAS space object (molecule's core energy)\n", " )\n", "\n", " # Outputs:\n", " # with_solvent_e : scalar energy correction due to solvent polarization\n", " # with_solvent_v : Fock-like matrix to be added to the core Hamiltonian in SCF\n", " with_solvent_e, with_solvent_v = mc.with_solvent._get_vind(\n", " rho_approximation\n", " )\n", " data.append((iiter, float(e), job_id))\n", " print(f\">>>>> END IITER {iiter}\")\n", " print(f\">>>>> TOTAL ENERGY: {e}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "id": "3905922f", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib.ticker import ScalarFormatter\n", "\n", "\n", "def plot_data(data, baseline=0, name=None, save=False):\n", " x_vals, y_vals, job_ids = zip(*data)\n", " fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", " # Plot line + markers\n", " ax.plot(\n", " x_vals,\n", " y_vals,\n", " color=\"navy\",\n", " linewidth=2,\n", " marker=\"o\",\n", " markersize=5,\n", " label=\"Energy trajectory\",\n", " )\n", " ax.axhline(\n", " baseline,\n", " color=\"red\",\n", " linestyle=\"--\",\n", " linewidth=1.5,\n", " label=\"Reference energy\",\n", " )\n", "\n", " # Force plain formatting\n", " ax.yaxis.set_major_formatter(ScalarFormatter(useMathText=True))\n", " ax.ticklabel_format(style=\"plain\", axis=\"y\")\n", "\n", " # Annotate each point with its exact value\n", " for x, y, job_id in zip(x_vals, y_vals, job_ids):\n", " ax.annotate(\n", " f\"{y:.8f}, ID: {job_id}\",\n", " (x, y),\n", " textcoords=\"offset points\",\n", " xytext=(0, 8), # vertical offset\n", " ha=\"center\",\n", " fontsize=8,\n", " rotation=25,\n", " color=\"navy\",\n", " )\n", "\n", " # Annotate the Classical Reference line\n", " for x, y in zip([0.5], [baseline]):\n", " ax.annotate(\n", " f\"{y:.5f}\",\n", " (x, y),\n", " textcoords=\"offset points\",\n", " xytext=(0, 8), # vertical offset\n", " ha=\"center\",\n", " fontsize=8,\n", " rotation=25,\n", " color=\"red\",\n", " )\n", "\n", " # Titles, labels, etc\n", " ax.set_title(\n", " f\"SQD/IEF-PCM(cc-pVDZ) - {name}\\nEnergy Convergence\",\n", " fontsize=14,\n", " fontweight=\"bold\",\n", " pad=15,\n", " )\n", " ax.set_xlabel(\"Update Iterations\", fontsize=12)\n", " ax.set_ylabel(\"Total Energy (Hartrees)\", fontsize=12)\n", " ax.grid(True, linestyle=\"--\", linewidth=0.6, alpha=0.7)\n", " ax.legend(frameon=True, loc=\"best\")\n", " plt.tight_layout()\n", "\n", " if save:\n", " plt.savefig(f\"./results/{name}_energy_convergence.png\")\n", " return fig, ax" ] }, { "cell_type": "code", "execution_count": null, "id": "6fdb639b", "metadata": {}, "outputs": [], "source": [ "# Plot your data\n", "\n", "fig, ax = plot_data(data, baseline=CASCI_E, name=molecule_name, save=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "75f48e6a-c7e4-46f3-9d39-a7a877427a04", "metadata": {}, "source": [ "\n", "\n", "## 次のステップ\n", "\n", "\n", " この作品に興味を持たれた方は、以下のリソースもご参照ください：\n", "\n", " * [Quantum Resource Management Interface（QRMI）](https://github.com/qiskit-community/qrmi) — 量子リソースと従来型リソースをSlurmなどのHPCワークロードマネージャーに統合し、本チュートリアルで示した分散コンピューティングのパターンを拡張します\n", " * [『 Qiskit Serverless 』ガイド](/docs/guides/serverless)\n", " * [Qiskit Serverless GitHub](https://qiskit.github.io/qiskit-serverless/index.html)\n", " * [化学ハミルトニアンのサンプルベース量子対角化](/docs/tutorials/sample-based-quantum-diagonalization)チュートリアル\n", " * [サンプルベース量子対角化（SQD）の概要](/docs/addons/qiskit-addon-sqd)\n", " * [暗黙的溶媒モデルを用いた電子構造シミュレーション用のテンプレートをデプロイして実行する](/docs/guides/function-template-chemistry-workflow) （SQD IEF-PCM Qiskit Functionテンプレート、クリーブランド・クリニックと IBM が共同開発）\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": 120 }, "nbformat": 4, "nbformat_minor": 5 }