{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "de193554-b271-4295-95e4-8904f0f6ee8a",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"分子幾何学\"\n",
        "description: \"このレッスンでは、単純な分子の構造を変化させ、各段階でエネルギーを最小化します。\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore pxxr prqs nelecas mcscf chmax Dmax vmax ecore ncas Excp disp */}\n",
        "\n",
        "<span id=\"determining-a-molecular-geometry\" />\n",
        "\n",
        "# 分子構造の決定\n",
        "\n",
        "前のセクションでは、分子の基底状態エネルギーを決定するためにVQEを実装した。 これは量子コンピューターの有効な利用法だが、それ以上に有用なのは分子の構造を決定することだろう。\n",
        "\n",
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "## ステップ1：古典的な入力を量子問題にマッピングする\n",
        "\n",
        "二原子水素の基本的な例に戻ると、変化する唯一の幾何学的パラメータは結合長である。 これを達成するために、前と同じように進めるが、最初の分子構築に変数（引数の結合長 *x* ）を使う。 これはかなり単純な変更だが、フェルミオンのハミルトニアンの構築から始まり、マッピングを経て最終的にコスト関数に伝搬するため、プロセス全体を通して変数を関数に含める必要がある。\n",
        "\n",
        "まず、前に使ったパッケージをいくつかロードし、コレスキー関数を定義する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "id": "0a5d39bd-8c26-404f-8057-c29e3af70df4",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.quantum_info import SparsePauliOp\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "#!pip install pyscf==2.4.0\n",
        "from pyscf import ao2mo, gto, mcscf, scf\n",
        "\n",
        "\n",
        "def cholesky(V, eps):\n",
        "    # see https://arxiv.org/pdf/1711.02242.pdf section B2\n",
        "    # see https://arxiv.org/abs/1808.02625\n",
        "    # see https://arxiv.org/abs/2104.08957\n",
        "    no = V.shape[0]\n",
        "    chmax, ng = 20 * no, 0\n",
        "    W = V.reshape(no**2, no**2)\n",
        "    L = np.zeros((no**2, chmax))\n",
        "    Dmax = np.diagonal(W).copy()\n",
        "    nu_max = np.argmax(Dmax)\n",
        "    vmax = Dmax[nu_max]\n",
        "    while vmax > eps:\n",
        "        L[:, ng] = W[:, nu_max]\n",
        "        if ng > 0:\n",
        "            L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])\n",
        "        L[:, ng] /= np.sqrt(vmax)\n",
        "        Dmax[: no**2] -= L[: no**2, ng] ** 2\n",
        "        ng += 1\n",
        "        nu_max = np.argmax(Dmax)\n",
        "        vmax = Dmax[nu_max]\n",
        "    L = L[:, :ng].reshape((no, no, ng))\n",
        "    print(\n",
        "        \"accuracy of Cholesky decomposition \",\n",
        "        np.abs(np.einsum(\"prg,qsg->prqs\", L, L) - V).max(),\n",
        "    )\n",
        "    return L, ng\n",
        "\n",
        "\n",
        "def identity(n):\n",
        "    return SparsePauliOp.from_list([(\"I\" * n, 1)])\n",
        "\n",
        "\n",
        "def creators_destructors(n, mapping=\"jordan_wigner\"):\n",
        "    c_list = []\n",
        "    if mapping == \"jordan_wigner\":\n",
        "        for p in range(n):\n",
        "            if p == 0:\n",
        "                ell, r = \"I\" * (n - 1), \"\"\n",
        "            elif p == n - 1:\n",
        "                ell, r = \"\", \"Z\" * (n - 1)\n",
        "            else:\n",
        "                ell, r = \"I\" * (n - p - 1), \"Z\" * p\n",
        "            cp = SparsePauliOp.from_list([(ell + \"X\" + r, 0.5), (ell + \"Y\" + r, -0.5j)])\n",
        "            c_list.append(cp)\n",
        "    else:\n",
        "        raise ValueError(\"Unsupported mapping.\")\n",
        "    d_list = [cp.adjoint() for cp in c_list]\n",
        "    return c_list, d_list"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "40e25b54-7927-4dbb-a26f-1c6b33f7f349",
      "metadata": {},
      "source": [
        "ここでハミルトニアンを定義するために、前の例とまったく同じように PySCF を使いますが、今度は原子間距離の役割を果たす変数 `x` を入れます。 これにより、コアエネルギー、単電子エネルギー、2電子エネルギーが以前と同じように返される。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "dbd10d0c-feb1-4a86-9bd9-b61101a08b95",
      "metadata": {},
      "outputs": [],
      "source": [
        "def ham_terms(x: float):\n",
        "    distance = x\n",
        "    a = distance / 2\n",
        "    mol = gto.Mole()\n",
        "    mol.build(\n",
        "        verbose=0,\n",
        "        atom=[\n",
        "            [\"H\", (0, 0, -a)],\n",
        "            [\"H\", (0, 0, a)],\n",
        "        ],\n",
        "        basis=\"sto-6g\",\n",
        "        spin=0,\n",
        "        charge=0,\n",
        "        symmetry=\"Dooh\",\n",
        "    )\n",
        "\n",
        "    # mf = scf.RHF(mol)\n",
        "    # mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))\n",
        "    # mx.kernel()\n",
        "\n",
        "    mf = scf.RHF(mol)\n",
        "    mf.kernel()\n",
        "    if not mf.converged:\n",
        "        raise RuntimeError(f\"SCF did not converge for distance {x}\")\n",
        "\n",
        "    mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))\n",
        "    casci_energy = mx.kernel()\n",
        "    if casci_energy is None:\n",
        "        raise RuntimeError(f\"CASCI failed for distance {x}\")\n",
        "\n",
        "    # Other variables that might come in handy:\n",
        "    # active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)\n",
        "    #    E1 = mf.kernel()\n",
        "    # mo = mx.sort_mo(active_space, base=0)\n",
        "    #    E2 = mx.kernel(mo)[:2]\n",
        "\n",
        "    h1e, ecore = mx.get_h1eff()\n",
        "    h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)\n",
        "    return ecore, h1e, h2e"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "db49a702-e60e-4c9b-8ff2-94aa2ada022c",
      "metadata": {},
      "source": [
        "上記の構成は、原子種、幾何学、電子軌道に基づいてフェルミオンのハミルトニアンを作っていることを思い出してほしい。 以下では、このフェルミオン的ハミルトニアンをパウリ作用素に写像する。 この `build_hamiltonian` 関数には、幾何学変数も引数として含まれる。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "84e6a56b-eaea-4c0a-8502-567cfc5140a2",
      "metadata": {},
      "outputs": [],
      "source": [
        "def build_hamiltonian(distx: float) -> SparsePauliOp:\n",
        "    ecore = ham_terms(distx)[0]\n",
        "    h1e = ham_terms(distx)[1]\n",
        "    h2e = ham_terms(distx)[2]\n",
        "\n",
        "    ncas, _ = h1e.shape\n",
        "\n",
        "    C, D = creators_destructors(2 * ncas, mapping=\"jordan_wigner\")\n",
        "    Exc = []\n",
        "    for p in range(ncas):\n",
        "        Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]\n",
        "        for r in range(p + 1, ncas):\n",
        "            Excp.append(\n",
        "                C[p] @ D[r]\n",
        "                + C[ncas + p] @ D[ncas + r]\n",
        "                + C[r] @ D[p]\n",
        "                + C[ncas + r] @ D[ncas + p]\n",
        "            )\n",
        "        Exc.append(Excp)\n",
        "\n",
        "    # low-rank decomposition of the Hamiltonian\n",
        "    Lop, ng = cholesky(h2e, 1e-6)\n",
        "    t1e = h1e - 0.5 * np.einsum(\"pxxr->pr\", h2e)\n",
        "\n",
        "    H = ecore * identity(2 * ncas)\n",
        "    # one-body term\n",
        "    for p in range(ncas):\n",
        "        for r in range(p, ncas):\n",
        "            H += t1e[p, r] * Exc[p][r - p]\n",
        "    # two-body term\n",
        "    for g in range(ng):\n",
        "        Lg = 0 * identity(2 * ncas)\n",
        "        for p in range(ncas):\n",
        "            for r in range(p, ncas):\n",
        "                Lg += Lop[p, r, g] * Exc[p][r - p]\n",
        "        H += 0.5 * Lg @ Lg\n",
        "\n",
        "    return H.chop().simplify()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ffcb2569-f832-4fa9-a8bb-d7626c2a233d",
      "metadata": {},
      "source": [
        "efficient\\_su2 ansatzや SciPy minimizersなど、VQE自体を実行するための残りのパッケージをロードする：\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "5f8dd8e0-6dfb-4be7-af15-1591f569201c",
      "metadata": {},
      "outputs": [],
      "source": [
        "# General imports\n",
        "\n",
        "# Pre-defined ansatz circuit and operator class for Hamiltonian\n",
        "from qiskit.circuit.library import efficient_su2\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "\n",
        "# SciPy minimizer routine\n",
        "from scipy.optimize import minimize\n",
        "\n",
        "# Plotting functions\n",
        "\n",
        "# qiskit-ibm-runtime tools\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "\n",
        "service = QiskitRuntimeService()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f94c1680-9c2a-4585-a63e-a49da4eb02f3",
      "metadata": {},
      "source": [
        "再びコスト関数を定義するが、これは常に完全に構築されマッピングされたハミルトニアンを引数として取るので、この関数については何も変わらない。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "0a80ad8c-a9cb-4cab-835d-f65131b99c87",
      "metadata": {},
      "outputs": [],
      "source": [
        "def cost_func(params, ansatz, H, estimator):\n",
        "    pub = (ansatz, [H], [params])\n",
        "    result = estimator.run(pubs=[pub]).result()\n",
        "    energy = result[0].data.evs[0]\n",
        "    return energy\n",
        "\n",
        "\n",
        "# def cost_func_sim(params, ansatz, H, estimator):\n",
        "#    energy = estimator.run(ansatz, H, parameter_values=params).result().values[0]\n",
        "#    return energy"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a068e645-8dfb-4e9c-b1b1-7dd936808188",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-execution\" />\n",
        "\n",
        "## ステップ2：量子実行のための問題最適化\n",
        "\n",
        "ハミルトニアンは新しいジオメトリごとに変わるので、演算子の転置はステップごとに変わる。 とはいえ、各ステップで適用される一般的なパスマネージャーを、使用したいハードウェアに合わせて定義することはできる。\n",
        "\n",
        "ここでは、利用可能なバックエンドの中で最もビジーでないものを使用する。 そのバックエンドを我々の AerSimulator, のモデルとして使い、シミュレーターが例えば実際のバックエンドのノイズの挙動を模倣できるようにする。 これらのノイズモデルは完璧なものではないが、実際のハードウェアに何を期待すればよいかを知る一助にはなるだろう。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "fa07518f-a2c0-4f7a-b344-8d9576427478",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Here, we select the least busy backend available:\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "print(backend)\n",
        "# Or to select a specific real backend use the line below, and substitute 'ibm_strasbourg'\n",
        "# for your chosen device. backend = service.get_backend('ibm_strasbourg')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "id": "e67fc84a-431b-4efd-937f-49c8b7ac3abb",
      "metadata": {},
      "outputs": [],
      "source": [
        "# To run on a simulator:\n",
        "# -----------\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "backend_sim = AerSimulator.from_backend(backend)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0636c55f-f46f-46ca-acaf-72bcc2f5f663",
      "metadata": {},
      "source": [
        "パス・マネージャーと関連パッケージをインポートし、サーキットの最適化に役立てる。 このステップとその上のステップは、ハミルトニアンから独立しているので、前回のレッスンから変更はない。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "id": "8202332e-69ca-4049-af27-e98a77f15a5d",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.transpiler import PassManager\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit.transpiler.passes import (\n",
        "    ALAPScheduleAnalysis,\n",
        "    PadDynamicalDecoupling,\n",
        "    ConstrainedReschedule,\n",
        ")\n",
        "from qiskit.circuit.library import XGate\n",
        "\n",
        "target = backend.target\n",
        "pm = generate_preset_pass_manager(target=target, optimization_level=3)\n",
        "pm.scheduling = PassManager(\n",
        "    [\n",
        "        ALAPScheduleAnalysis(target=target),\n",
        "        ConstrainedReschedule(\n",
        "            acquire_alignment=target.acquire_alignment,\n",
        "            pulse_alignment=target.pulse_alignment,\n",
        "            target=target,\n",
        "        ),\n",
        "        PadDynamicalDecoupling(\n",
        "            target=target,\n",
        "            dd_sequence=[XGate(), XGate()],\n",
        "            pulse_alignment=target.pulse_alignment,\n",
        "        ),\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4932ef5e-3c81-46a1-92cb-3082399734a0",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-ibm-quantum-primitives\" />\n",
        "\n",
        "## ステップ 3： IBM Quantum プリミティブを使用して実行します。\n",
        "\n",
        "以下のコードブロックでは、原子間距離の $x$ における各ステップの出力値を格納するための配列を設定しています。 $x$ の範囲は、平衡結合長の実験値である 0.74 オングストロームという知見に基づいて選択しました。 まずはシミュレータ上で実行するため、. `qiskit.primitives`からEstimator（ BackendEstimator ）をインポートします。 各幾何学的ステップについて、ハミルトニアンを構築し、最適化ツール「cobyla」を用いて一定回数（ここでは500回）の最適化ステップを実行する。 各幾何学的ステップにおいて、総エネルギーと電子エネルギーの両方を保存する。 オプティマイザの処理ステップが多いため、1時間以上かかる可能性があります。 所要時間を短縮するために、以下の入力内容を変更することをお勧めします。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "id": "4c41a221-8a02-4932-882b-5afcc98d1d8d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "accuracy of Cholesky decomposition  1.1102230246251565e-15\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/scipy/_lib/pyprima/common/preproc.py:68: UserWarning: COBYLA: Invalid MAXFUN; it should be at least num_vars + 2; it is set to 34\n",
            "  warn(f'{solver}: Invalid MAXFUN; it should be at least {min_maxfun_str}; it is set to {maxfun}')\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 1.316011435623847\n",
            "The corresponding X is:\n",
            "[2.32948769 5.39918229 3.03787975 4.11789904 4.97130735 2.68662232\n",
            " 1.76573151 2.48982571 5.40431972 3.65780829 1.33792786 5.48472494\n",
            " 6.18738702 1.78741883 0.78195251 2.96658955 1.35827677 5.599321\n",
            " 4.54850148 1.0939048  4.26158726 0.52100721 0.82318    4.76796961\n",
            " 3.75795507 3.8526447  5.51100375 5.91023075 2.61494836 1.79908918\n",
            " 2.65937756 5.53964148]\n",
            "\n",
            "-0.44791260077615314\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 1.316011435623847\n",
            "       x: [ 2.329e+00  5.399e+00 ...  2.659e+00  5.540e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.7235003672327549\n",
            "The corresponding X is:\n",
            "[2.56282915 5.63369524 5.58059887 4.049643   4.2021266  3.06866011\n",
            " 6.01619635 1.52520776 4.35403161 0.33673958 0.32623161 1.2179545\n",
            " 2.84001371 3.98956684 4.89632562 1.38303588 1.96194695 2.13182089\n",
            " 0.29739166 1.77895165 3.29151585 3.54355374 4.49626674 0.95756626\n",
            " 0.87103927 4.53068385 1.31051302 0.37103108 1.02961355 3.13342311\n",
            " 5.65815319 2.24770604]\n",
            "\n",
            "-0.5994426600672451\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.7235003672327549\n",
            "       x: [ 2.563e+00  5.634e+00 ...  5.658e+00  2.248e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.34960914928810116\n",
            "The corresponding X is:\n",
            "[5.44143165 6.75955835 1.56836472 3.09522093 4.67873235 1.67071481\n",
            " 0.3056494  0.65998337 1.02197668 5.21162959 0.43690354 3.56522934\n",
            " 4.56033119 1.90736037 0.40863891 2.87007312 3.2516952  5.90360196\n",
            " 1.99057799 5.20726456 0.74710237 6.03179202 3.80685028 0.03844391\n",
            " 5.88580196 3.62233258 3.98723567 2.50591888 5.44020267 2.2792993\n",
            " 5.57102303 4.46548617]\n",
            "\n",
            "-0.7087452725518989\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.34960914928810116\n",
            "       x: [ 5.441e+00  6.760e+00 ...  5.571e+00  4.465e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  2.220446049250313e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.10594558882184543\n",
            "The corresponding X is:\n",
            "[5.35675483 2.26629567 1.45430546 5.56758296 5.76309509 0.73239338\n",
            " 5.1216998  3.03258872 4.33624828 1.93197674 0.5292902  3.32274987\n",
            " 3.43247633 0.81490741 0.48060245 1.9944799  5.67519646 5.12534057\n",
            " 0.06510627 2.52989834 6.1699519  0.94828957 5.91634548 1.5994961\n",
            " 4.27902164 2.3129213  1.82353095 2.10634209 1.43740426 4.06988733\n",
            " 0.59624074 4.93925418]\n",
            "\n",
            "-0.7760164293781545\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.10594558882184543\n",
            "       x: [ 5.357e+00  2.266e+00 ...  5.962e-01  4.939e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.06473600797229297\n",
            "The corresponding X is:\n",
            "[6.07735568 0.18019501 0.20743128 4.15445985 3.59388894 5.10047555\n",
            " 6.09938474 6.54707528 3.36251167 2.05475223 3.67078456 5.96010605\n",
            " 2.58589996 5.2723619  3.26352977 2.47432334 3.50289983 2.06620525\n",
            " 6.0946056  1.22751903 0.97320057 2.19564095 5.73174941 2.05127682\n",
            " 5.73805165 3.84046105 1.84816963 2.1247504  3.11106736 2.44136052\n",
            " 3.39002685 0.81596991]\n",
            "\n",
            "-0.8207034521437214\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.06473600797229297\n",
            "       x: [ 6.077e+00  1.802e-01 ...  3.390e+00  8.160e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-17\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.19562982094782935\n",
            "The corresponding X is:\n",
            "[-0.02184462  3.67041038  7.25918653  5.89799546  0.63583624  1.84214506\n",
            "  2.84059837  5.31485182  1.6053784   0.04556618  0.32018993 -0.03884066\n",
            "  0.69131496  0.24203727  1.97397262  3.59723495  0.43355775  2.30131056\n",
            "  4.63482292  3.9857415   4.32320753  4.55388437  2.18753433  5.99034987\n",
            "  2.50489913  0.90650534  4.82518088  2.32954849  2.29901832  5.33658863\n",
            "  5.91246716  3.2405013 ]\n",
            "\n",
            "-0.8571013345978292\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.19562982094782935\n",
            "       x: [-2.184e-02  3.670e+00 ...  5.912e+00  3.241e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.2833766309947055\n",
            "The corresponding X is:\n",
            "[ 3.1700088   5.05055456  1.2545611   4.28751811  0.6255103   1.67526577\n",
            "  5.48201473  4.83820497  7.34880059  5.99705431  4.2502643   0.32066274\n",
            "  0.41001404  0.27271241  4.15682546  4.22393693  4.35148115  0.64538137\n",
            "  5.26288622  5.03810489  4.62426621  4.74997689  1.09603919  0.34752466\n",
            "  1.8116275   0.7474807   5.31754143  4.11181763  1.58797998  5.6299796\n",
            "  3.0109383  -0.19062772]\n",
            "\n",
            "-0.8713513097947054\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.2833766309947055\n",
            "       x: [ 3.170e+00  5.051e+00 ...  3.011e+00 -1.906e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.3527503628484244\n",
            "The corresponding X is:\n",
            "[3.90513622 4.61398739 5.92552705 1.99953405 4.82157369 1.35702441\n",
            " 2.77701782 5.73612247 4.22710527 1.83463189 0.45796297 4.62509318\n",
            " 0.98998668 0.11666217 3.0234641  4.54298546 0.14034033 4.15635797\n",
            " 1.41257357 4.48719602 2.39365535 0.19672041 5.0763044  1.86357581\n",
            " 3.657757   4.60298344 2.49769577 1.88086199 3.00108725 1.84475841\n",
            " 5.24047385 4.91142914]\n",
            "\n",
            "-0.8819275737684243\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.3527503628484244\n",
            "       x: [ 3.905e+00  4.614e+00 ...  5.240e+00  4.911e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  2.7755575615628914e-17\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.4022181851996095\n",
            "The corresponding X is:\n",
            "[6.09453981 3.5109422  3.37216019 4.94732621 1.25662002 5.89645164\n",
            " 5.06403334 2.68073141 4.40385083 1.13638366 1.73347762 6.82932871\n",
            " 1.15265014 2.07145964 4.36520459 1.14960341 1.62288871 4.32315915\n",
            " 5.45622821 0.93554005 3.17418483 0.47230243 1.31535502 5.77698726\n",
            " 2.04927925 2.50663538 5.9706002  5.4984681  2.9421232  1.56636313\n",
            " 1.09394523 4.62582   ]\n",
            "\n",
            "-0.8832883769450639\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.4022181851996095\n",
            "       x: [ 6.095e+00  3.511e+00 ...  1.094e+00  4.626e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.44423031870708934\n",
            "The corresponding X is:\n",
            "[4.05765050e+00 3.99144950e+00 3.13287593e+00 3.28855137e+00\n",
            " 4.32613515e+00 4.91104512e+00 1.86521867e+00 2.18822879e+00\n",
            " 6.01336171e+00 1.82501276e+00 2.64830637e+00 5.53045823e+00\n",
            " 2.36110093e+00 3.98821703e+00 4.69013438e-01 4.38996815e+00\n",
            " 7.78103801e-04 1.72994378e+00 2.24970934e+00 1.11978200e+00\n",
            " 2.24846445e+00 4.90745512e+00 5.38474921e+00 5.03587994e+00\n",
            " 3.54297277e+00 4.78147533e+00 1.25990218e+00 1.99168068e+00\n",
            " 5.89203503e+00 1.77673987e+00 5.37848357e+00 5.60245198e-01]\n",
            "\n",
            "-0.8852113278070892\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.44423031870708934\n",
            "       x: [ 4.058e+00  3.991e+00 ...  5.378e+00  5.602e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "All energies have been calculated\n"
          ]
        }
      ],
      "source": [
        "from qiskit.primitives import BackendEstimatorV2\n",
        "\n",
        "estimator = BackendEstimatorV2(backend=backend_sim)\n",
        "\n",
        "distances_sim = np.arange(0.3, 1.3, 0.1)\n",
        "vqe_energies_sim = []\n",
        "vqe_elec_energies_sim = []\n",
        "\n",
        "for dist in distances_sim:\n",
        "    xx = dist\n",
        "\n",
        "    # Random initial state and efficient_su2 ansatz\n",
        "    H = build_hamiltonian(xx)\n",
        "    ansatz = efficient_su2(H.num_qubits)\n",
        "    ansatz_isa = pm.run(ansatz)\n",
        "    x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)\n",
        "    H_isa = H.apply_layout(ansatz_isa.layout)\n",
        "    nuclear_repulsion = ham_terms(xx)[0]\n",
        "\n",
        "    res = minimize(\n",
        "        cost_func,\n",
        "        x0,\n",
        "        args=(ansatz_isa, H_isa, estimator),\n",
        "        method=\"cobyla\",\n",
        "        options={\"maxiter\": 20, \"disp\": True},\n",
        "    )\n",
        "\n",
        "    # Note this returns the total energy, and we are often interested in the electronic energy\n",
        "    tot_energy = getattr(res, \"fun\")\n",
        "    electron_energy = getattr(res, \"fun\") - nuclear_repulsion\n",
        "    print(electron_energy)\n",
        "    vqe_energies_sim.append(tot_energy)\n",
        "    vqe_elec_energies_sim.append(electron_energy)\n",
        "\n",
        "    # Print all results\n",
        "    print(res)\n",
        "\n",
        "print(\"All energies have been calculated\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "id": "bc6164ca-6909-4780-8009-6dc274c66268",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(1.2000000000000004)"
            ]
          },
          "execution_count": 39,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "xx"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "633a2a89-950f-4d13-a56e-6adc079245ea",
      "metadata": {},
      "source": [
        "この出力の結果については、後述の「後処理」のセクションで詳しく説明しますが、ここではひとまず、シミュレーションが成功したことに留めておきます。 これで、実際のハードウェア上で実行する準備が整いました。 レジリエンスを に設定し `1`、TREXエラー緩和機能を使用することを指定します。 実際にハードウェアを使って作業することになったので、 IBM Quantum のCompute Serviceと IBM Quantum のプリミティブを使用します。 なお、ジオメトリに関連するforループと、複数の変分試行の両方が、このセッション内に含まれていることに注意してください。\n",
        "\n",
        "実際のハードウェアでの実行にはコストと時間の制限があるため、以下ではジオメトリステップとオプティマイザステップの数を減らしている。 精度の目標や制限時間に応じて、これらのステップを調整するようにしてください。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b09744f5-87a9-4744-b495-cf993e5ffcb3",
      "metadata": {},
      "outputs": [],
      "source": [
        "# To continue running on real hardware use\n",
        "from qiskit_ibm_runtime import Session\n",
        "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n",
        "from qiskit_ibm_runtime import EstimatorOptions\n",
        "\n",
        "estimator_options = EstimatorOptions(resilience_level=1, default_shots=2000)\n",
        "\n",
        "distances = np.arange(0.5, 0.9, 0.1)\n",
        "vqe_energies = []\n",
        "vqe_elec_energies = []\n",
        "\n",
        "with Session(backend=backend) as session:\n",
        "    estimator = Estimator(mode=session, options=estimator_options)\n",
        "\n",
        "    for dist in distances:\n",
        "        xx = dist\n",
        "\n",
        "        # Random initial state and efficient_su2 ansatz\n",
        "\n",
        "        H = build_hamiltonian(xx)\n",
        "        ansatz = efficient_su2(H.num_qubits)\n",
        "        ansatz_isa = pm.run(ansatz)\n",
        "        H_isa = H.apply_layout(ansatz_isa.layout)\n",
        "        nuclear_repulsion = ham_terms(xx)[0]\n",
        "        x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)\n",
        "\n",
        "        res = minimize(\n",
        "            cost_func,\n",
        "            x0,\n",
        "            args=(ansatz_isa, H_isa, estimator),\n",
        "            method=\"cobyla\",\n",
        "            options={\"maxiter\": 50, \"disp\": True},\n",
        "        )\n",
        "\n",
        "        # Note this returns the total energy, and we are often interested in the electronic energy\n",
        "        tot_energy = getattr(res, \"fun\")\n",
        "        electron_energy = getattr(res, \"fun\") - nuclear_repulsion\n",
        "        print(electron_energy)\n",
        "        vqe_energies.append(tot_energy)\n",
        "        vqe_elec_energies.append(electron_energy)\n",
        "\n",
        "        # Print all results\n",
        "        print(res)\n",
        "\n",
        "print(\"All energies have been calculated\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5b6f1491-88bd-4fb1-a1fe-2e29ffa33f17",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-processing\" />\n",
        "\n",
        "## ステップ4：後処理\n",
        "\n",
        "シミュレーターと実際のハードウェアの両方について、各原子間距離について計算された基底状態エネルギーをプロットし、どこで最も低いエネルギーが得られるかを見ることができる。 これが自然界で見られる原子間距離のはずであり、実際に近い。 より滑らかな曲線は、他の解析やオプティマイザを試したり、各ジオメトリステップで複数回計算を実行し、いくつかのランダムな初期条件で平均化することで得られるかもしれない。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e81f3ead-27ac-415e-a9f2-64a51d4b7aa3",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/quantum-chem-with-vqe/geometry/extracted-outputs/e81f3ead-27ac-415e-a9f2-64a51d4b7aa3-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Here we can plot the results from this simulation.\n",
        "plt.plot(distances_sim, vqe_energies_sim, label=\"VQE Energy\")\n",
        "plt.xlabel(\"Atomic distance (Angstrom)\")\n",
        "plt.ylabel(\"Energy\")\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1efe0d76-3c6c-4369-9ace-fc890084e676",
      "metadata": {},
      "source": [
        "最適化のステップ数を単純に増やしても、シミュレータの結果が改善される可能性は高くない。\n",
        "\n",
        "サンプリングされた値の範囲がわずかに異なることを除けば、実際のハードウェアからの結果は同等である。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "de8f53cf-9547-4578-b6cb-20d2b5602ee0",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/quantum-chem-with-vqe/geometry/extracted-outputs/de8f53cf-9547-4578-b6cb-20d2b5602ee0-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(distances, vqe_energies, label=\"VQE Energy\")\n",
        "plt.xlabel(\"Atomic distance (Angstrom)\")\n",
        "plt.ylabel(\"Energy\")\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7e732366-425a-40e4-b6f2-d26f11f7d86b",
      "metadata": {},
      "source": [
        "H2、結合長は 0.74 オングストローム、全エネルギーは -1.17 ハートリーである。 実際のハードウェアの結果は、シミュレーターよりもこれらの値に近いことがわかる。 これは、どちらのケースでもノイズが存在した（あるいはシミュレートされた）ためと思われるが、実際のハードウェアの場合のみエラー緩和が行われた。\n",
        "\n",
        "<span id=\"closing\" />\n",
        "\n",
        "### クローズ中\n",
        "\n",
        "これで量子化学のためのVQE講座は終わりです。 量子コンピューティングで使用される基本的な情報理論に興味がある方は、ジョン・ワトラスの[量子情報の基礎に関する](/learning/courses/basics-of-quantum-information)コースをご覧ください。 VQEワークフローの簡単な例については、 [VQEによるハイゼンベルグ鎖の基底状態エネルギー推定チュートリアルを](/docs/tutorials/spin-chain-vqe)参照してください。 また、 [チュートリアルや](/docs/tutorials) [コースを](/learning)閲覧して、量子コンピューティングの最新技術に関する教材を探すこともできます。\n",
        "\n",
        "このコースの試験を受けることをお忘れなく。 80%以上のスコアを獲得すると、Credlyバッジが付与され、自動的にEメールで送信されます。 IBM クォンタム® ネットワークにご参加いただきありがとうございます！\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0ac6d249-d210-479d-a792-c8b4e94b8b88",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1.3.2\n",
            "0.35.0\n"
          ]
        }
      ],
      "source": [
        "import qiskit\n",
        "import qiskit_ibm_runtime\n",
        "\n",
        "print(qiskit.version.get_version_info())\n",
        "print(qiskit_ibm_runtime.version.get_version_info())"
      ]
    },
    {
      "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": 2
}