content/ex5/ex5-ja.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 演習5 - VQE(変分量子固有値ソルバー)\n",
    "***\n",
    "\n",
    "## 歴史的背景\n",
    "\n",
    "\n",
    "過去10年間で、量子コンピューターは急速に成熟し、量子的な手法で自然の法則をシミュレートできるコンピューティングシステムというファインマンの夢を実現し始めました。2014年の論文において、最初に、アルベルト・ペルッゾが **変分量子固有値ソルバー(VQE)** を発表しました。分子の基底状態エネルギー(最小エネルギー)をこれまでの手法より短い回路で見つけるアルゴリズムです。[1]　そして、2017年に、IBMの量子チームがVQEアルゴリズムを使って水素化リチウム分子の基底状態をシミュレートしました。[2]\n",
    "\n",
    "VQEのマジックは、問題の計算ワークロードのうちの一部を古典コンピューターにアウトソースすることです。アルゴリズムは、まず試行状態(ansatz、ベストな推測)と呼ばれるパラメーター化された量子回路から始め、古典オプティマイザーを使ってこの回路の最適なパラメーターを探します。VQEが古典アルゴリズムより優っている点は、量子回路が問題の厳密な波動関数を表現し保存できるということです。これは古典コンピューターでは指数関数的に難しい問題です。\n",
    "\n",
    "この演習５では、分子の基底状態と基底エネルギーを決定するために、変分量子固有値ソルバーを設定することで、ファインマンの夢をみなさんに実現してもらいます。この問題は、基底状態は、様々な分子の特性を計算するために使われるので興味深いことです。例えば、原子核における厳密な力は分子動力学シミュレーションで化学システムにおいて何が起こっているのか時間変化を伴って探究することができます。[3]\n",
    "\n",
    "\n",
    "### 参考文献\n",
    "\n",
    "1. Peruzzo, Alberto, et al. \"A variational eigenvalue solver on a photonic quantum processor.\" Nature communications 5.1 (2014): 1-7.\n",
    "2. Kandala, Abhinav, et al. \"Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.\" Nature 549.7671 (2017): 242-246.\n",
    "3. Sokolov, Igor O., et al. \"Microcanonical and finite-temperature ab initio molecular dynamics simulations on quantum computers.\" Physical Review Research 3.1 (2021): 013125.\n",
    "\n",
    "\n",
    "## はじめに\n",
    "\n",
    "VQEの実装において、特に試行状態(ansatz)の量子回路に注目して、どのようにシミュレーションを実装するか、その方法を選びます。\n",
    "ノイズのある量子コンピューターにおいてVQEを実行する際に重要な課題の一つは、忠実度(fidelity)をロスしないように、最もコンパクトな量子回路を基底状態として選ぶことです。この問題は、この課題に基づいて作られています。\n",
    "この課題は、実際には、精度を損わずに、２量子ビットゲート（例：CNOT）の数と変分パラメーターの数を削減することを意味します。\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "<b>目標</b> \n",
    "    \n",
    "与えられた問題における基底状態を正確に表現する最も短い試行状態(ansatz)を見つけてください。創造的に！\n",
    "    \n",
    "    \n",
    "<b>計画</b>   \n",
    "    \n",
    "はじめに、小さな分子を使ってVQEシミュレーションの構築の仕方を学びます。その後、より大きな分子の場合について学んだことを適用します。\n",
    "\n",
    "**1. チュートリアル -  H$_2$に対するVQE:** VQEに慣れるために、statevectorシミュレーターで実行して、試行状態(ansatz)と古典オプティマイザーのベストな組み合わせを選んでください。\n",
    "    \n",
    "    \n",
    "**2. チャレンジ - LiHに対するVQE:** 1のチュートリアルと同じような検討をしますが、statevectorシミュレーターのみに制限します。Qiskitに用意されている量子ビット数を削減するスキームを使って、このより大きな系に対して最適な回路を探してください。回路を最適化し、想像力を使って、パラメーター化された回路のベストなビルティングブロックを選ぶ方法を探してください。そして、Qiskitにすでにある基底状態のための試行回路よりもコンパクトで、最もコンパクトな試行状態回路を構築してください。\n",
    "    \n",
    "    \n",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "以下はVQEシミュレーションの理論の紹介です。VQEの実行の前に全てを理解する必要はありません。怖がらないで！    \n",
    "\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***\n",
    "\n",
    "## 理論\n",
    "\n",
    "\n",
    "下図に量子コンピューター上でVQEを使って分子シミュレーションを行う一般的なワークフローを示します。\n",
    "\n",
    "<img src=\"resources/workflow.png\" width=800 height= 1400/>\n",
    "\n",
    "量子-古典のハイブリッド手法のアイディアのコアは、 **CPU(古典プロセッシング・ユニット)** と **QPU(量子プロセッシング・ユニット)** にそれぞれベストな計算ができる部分をアウトソースすることです。CPUは、エネルギー計算のために測定する必要のある項目をリストすること、また回路のパラメーターを最適化することを担当します。QPUは、システムの量子状態を表現する量子回路を実装し、エネルギーを測定します。より詳細には以下のようになります：\n",
    "\n",
    "**CPU** は、電子のホッピングと相互作用(ハートリー・フォック計算による1体/ 2体積分)に関連したエネルギーを効率的に計算することができます。このエネルギーは、全エネルギーを表す演算子であるハミルトニアンとして表されます。[ハートリー・フォック (HF) 法](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method#:~:text=In%20computational%20physics%20and%20chemistry,system%20in%20a%20stationary%20state.) は、波動関数を1つのスレーター行列式によって表すことができると仮定することで、効率的に基底状態の波動関数の近似計算を行います。（例：4スピン軌道と4量子ビットのSTO-3G基底におけるH$_2$分子において、最もエネルギーの低いスピン軌道を電子が占有している場合、$|\\Psi_{HF} \\rangle = |0101 \\rangle$ です。）\n",
    "QPUが後でVQEで行うことは、欠落している電子相関に関連する他の状態を表すこともできる(回路とそのパラメーターに関連した)量子状態を見つけることです。(例：$|\\Psi \\rangle$における$\\sum_i c_i |i\\rangle$は、$c_{HF}|\\Psi_{HF} \\rangle + \\sum_i c_i |i\\rangle $に等しい。ここで$i$はビット列です。)\n",
    "\n",
    "HF計算の後、ハミルトニアンにおける演算子は、フェルミオン-量子ビット変換を使ってQPUにおける測定量にマップされます。（後述のハミルトニアンの章を参照してください。）\n",
    "量子ビット数の削減や試行回路(ansatz)を短くするためにシステムの特性を更に分析することができます：\n",
    "\n",
    "- Z2対称性と２量子ビット削減のためには、こちらをご覧ください：[Bravyi *et al*, 2017](https://arxiv.org/abs/1701.08213v1)\n",
    "- エンタングルメントを作るためには、こちらをご覧ください： [Eddins *et al.*, 2021](https://arxiv.org/abs/2104.10220v1)\n",
    "- 試行回路(ansatz)の適応のためには、こちらをご覧ください： [Tang *et al.*,2019](https://arxiv.org/abs/1911.10205). これらの論文のアイディアを活用して、量子回路を短くする方法を見つけることができるでしょう。\n",
    "\n",
    "**QPU** は、角度$\\vec\\theta$でパラメーター化された量子回路(下の試行状態の章を参照)を構築します。この回路は、様々な単一量子ビット回転とエンタングラー(2量子ビットゲートなど)を活用することで、基底状態の波動関数を表現します。量子超越性は、QPUが効率的に正確な波動関数を表現し、保持できることが元になっています。古典コンピューターでは、原子数が数個を超えるシステムの波動関数は扱いにくくなります。最終的には、QPUは選択された演算子(例：ハミルトニアンを表現する演算子)を測定します。\n",
    "\n",
    "以下で、VQEアルゴリズムの各コンポーネントにおける数学的な詳細について少し詳しくみていきます。[VQEのエピソード動画](https://www.youtube.com/watch?v=Z-A6G0WVI9w)も役に立つと思います\n",
    "\n",
    "\n",
    "### ハミルトニアン\n",
    "\n",
    "ここでは、与えられた系に対して測定する必要のある演算子をどうやって求めるかを説明します。\n",
    "これらは以下のように定義された分子のハミルトニアンに含まれています：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{H} &=\\sum_{r s} h_{r s} \\hat{a}_{r}^{\\dagger} \\hat{a}_{s} \\\\\n",
    "&+\\frac{1}{2} \\sum_{p q r s} g_{p q r s} \\hat{a}_{p}^{\\dagger} \\hat{a}_{q}^{\\dagger} \\hat{a}_{r} \\hat{a}_{s}+E_{N N}\n",
    "\\end{aligned}\n",
    "$$\n",
    "ここで\n",
    "$$\n",
    "h_{p q}=\\int \\phi_{p}^{*}(r)\\left(-\\frac{1}{2} \\nabla^{2}-\\sum_{I} \\frac{Z_{I}}{R_{I}-r}\\right) \\phi_{q}(r)\n",
    "$$\n",
    "$$\n",
    "g_{p q r s}=\\int \\frac{\\phi_{p}^{*}\\left(r_{1}\\right) \\phi_{q}^{*}\\left(r_{2}\\right) \\phi_{r}\\left(r_{2}\\right) \\phi_{s}\\left(r_{1}\\right)}{\\left|r_{1}-r_{2}\\right|} \n",
    "$$\n",
    "\n",
    "\n",
    "ここで$h_{r s}$と$g_{p q r s}$は、1体/2体の電子積分(ハートリー・フォック法で求める)であり、$E_{N N}$は、核反発エネルギーです。\n",
    "1体の電子積分は、電子の運動エネルギーと電子と原子核との相互作用を表します。\n",
    "2体の電子積分は、電子と電子の相互作用を表します。\n",
    "$\\hat{a}_{r}^{\\dagger}, \\hat{a}_{r}$の演算子は、スピン軌道$r$における電子の生成と消滅を表し、量子コンピューターでそれらを測定できるように演算子へマップされる必要があります。VQEは、電子のエネルギーを最小化するので、全体のエネルギーを計算するために、核反発ネルギー$E_{NN}$を取得して追加する必要があることに注意してください。\n",
    "\n",
    "よって、$ h_{r s}$と$g_{p q r s}$のテンソルのゼロでない要素すべてにおいて、以下のフェルミオンから量子ビットへの変換に従って、パウリ・ストリング(パウリ演算子のテンソル積)を構築することができます。例えば、軌道$r = 3$における、ジョルダン・ウィグナー変換で、以下のパウリ・ストリングを得ることができます：\n",
    "\n",
    "$$\n",
    "\\hat a_{3}^{\\dagger}= \\hat \\sigma_z \\otimes \\hat \\sigma_z \\otimes\\left(\\frac{ \\hat \\sigma_x-i \\hat \\sigma_y}{2}\\right) \\otimes 1 \\otimes \\cdots \\otimes 1\n",
    "$$\n",
    "\n",
    "ここで、$\\hat \\sigma_x, \\hat \\sigma_y, \\hat \\sigma_z$は、よく知られているパウリ演算子です。$\\hat \\sigma_z$演算子のテンソル積は、フェルミオンの反交換関係を満たすために配置されます。\n",
    "水分子の14スピン軌道と14量子ビットの間のジョルダン・ウィグナー変換は、以下のように与えられます：\n",
    "\n",
    "<img src=\"resources/mapping.png\" width=600 height= 1200/>\n",
    "\n",
    "\n",
    "従って、ハミルトニアンにおける1体/2体の励起(例えば、 $\\hat{a}_{r}^{\\dagger} \\hat{a}_{s}$, $\\hat{a}_{p}^{\\dagger} \\hat{a}_{q}^{\\dagger} \\hat{a}_{r} \\hat{a}_{s}$)を、対応するパウリ・ストリング(つまり、上の図の$\\hat{P}_i$)によって置き換えます。\n",
    "その結果、QPUによって測定可能な演算子のセットとなります。\n",
    "より詳しい情報は、[Seeley *et al.*, 2012](https://arxiv.org/abs/1208.5986v1)をご覧ください。\n",
    "\n",
    "### 試行状態(Ansatzes)\n",
    "\n",
    "化学の問題には主に２種類の試行状態(Ansatzes)を使うことができます。\n",
    "\n",
    "- **量子ユニタリー結合クラスター(q-UCC)試行状態** は、物理的な発想によるもので、量子回路に電子の励起を大まかにマップします。量子ユニタリー結合クラスターシングル・ダブル(q-UCCSD)試行状態 (Qiskitでは`UCCSD`)は、すべての可能な一電子励起と二電子励起まで考慮します。ペアダブルq-pUCCD (`PUCCD`) とシングレットのq-UCCD0 (`SUCCD`)は、これらの励起の一部のみを考慮し(したがってとても短い回路になります。)、解離プロファイルに良い結果をもたらすことが証明されています。例えば、q-pUCCD は、以下の図のように、一励起状態を持たず、第二励起状がペアになっています。\n",
    "\n",
    "- **経験的な試行状態(TwoLocal)** は、回路の深さ(Depth)を短くするために発明されましたが、この短い回路を使っても基底状態を表現できます。\n",
    "下の図にあるように、Rゲートがパラメーター化された1量子ビットの回転で、２量子ビットゲートによる$U_{CNOT}$ のエンタングラーがあります。特定の$D$回、（独立なパラメーターを使って）この同じブロックを繰り返した後に、基底状態に達するというアイディアです。\n",
    "\n",
    "より詳しい説明は、[Sokolov *et al.* (q-UCC ansatzes)](https://arxiv.org/abs/1911.10864v2) 、[Barkoutsos *et al.*](https://arxiv.org/pdf/1805.04340.pdf)をご覧ください。\n",
    "\n",
    "<img src=\"resources/ansatz.png\" width=700 height= 1200/>\n",
    "\n",
    "\n",
    "### VQE\n",
    "\n",
    "ハミルトニアン演算子$\\hat H$が、固有ベクトル$|\\psi_{min}\\rangle$に関連付けられた未知の最小固有値$E_{min}$を持つとして与えられた時、VQEは、$E_{\\theta}$を推定し、その最小値は、$E_{min}$です：\n",
    "\n",
    "\\begin{align*}\n",
    "    E_{min} \\le E_{\\theta} \\equiv \\langle \\psi(\\theta) |\\hat H|\\psi(\\theta) \\rangle\n",
    "\\end{align*}  \n",
    "\n",
    "ここで、$|\\psi(\\theta)\\rangle$は、$E_{\\theta}$に関連付けられた試行状態です。\n",
    "$U(\\theta)$で表された、パラメーター化された回路を任意の初期状態$|\\psi\\rangle$に対して適用することによって、アルゴリズムは、$U(\\theta)|\\psi\\rangle \\equiv |\\psi(\\theta)\\rangle$を$|\\psi_{min}\\rangle$に推定します。\n",
    "推定は、古典オプティマイザーによってパラメーター$\\theta$を変え、期待値$\\langle \\psi(\\theta) |\\hat H|\\psi(\\theta) \\rangle$を最小化するように、繰り返し最適化されます。\n",
    "\n",
    "VQEの応用として、分子動力学シミュレーションにおける可能性についてはこちら[Sokolov *et al.*, 2021](https://arxiv.org/abs/2008.08144v1)をご覧ください。また励起状態の計算については、こちら[Ollitrault *et al.*, 2019](https://arxiv.org/abs/1910.12890) にいくつか例があります。\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    " \n",
    "<b> より詳細な説明</b> \n",
    "    \n",
    "このアルゴリズムを実装した、Qiskit Natureのチュートリアルは[こちら](https://qiskit.org/documentation/nature/tutorials/01_electronic_structure.html)にありますが、十分ではないので、[githubレポジトリーの最初のページ](https://github.com/Qiskit/qiskit-nature)と[テストフォルダー](https://github.com/Qiskit/qiskit-nature/tree/main/test) をご覧になることをお勧めします。各コンポーネントに対するテストがあり、各機能の基本となるコードが置かれています。\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "# パート1: チュートリアル - H$_2$ 分子に対するVQE\n",
    "\n",
    "    \n",
    "このパートでは、PySCFドライバーとジョルダン・ウィグナー変換で、STO-3G基底を用いて、H$_2$分子をシミュレートします。\n",
    "以下でやり方を紹介するので、より難しい問題にも取り組めるようになるはずです。\n",
    "\n",
    "#### 1. ドライバー\n",
    " \n",
    "\n",
    "Qiskitに用意されている古典的な化学のコードへのインターフェースをドライバーと呼びます。例えば、`PSI4Driver`, `PyQuanteDriver`, `PySCFDriver`が用意されています。\n",
    "\n",
    "以下のセルで、ドライバーを実行させ(与えられた基底セットと分子構造においてハートリー・フォック法の計算をして)、量子アルゴリズムに適用する分子に関する必要なすべての情報を取得します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers import PySCFDriver\n",
    "\n",
    "molecule = \"H .0 .0 .0; H .0 .0 0.739\"\n",
    "driver = PySCFDriver(atom=molecule)\n",
    "qmolecule = driver.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "<b> 練習問題 1</b> \n",
    "    \n",
    "`qmolecule`の特性を調べて、以下の質問に答えましょう。\n",
    "    \n",
    "1. 分子の基本的特性を知る必要があります。与えられた系の電子は合計何個ですか？\n",
    "2. 分子軌道は何個あります？    \n",
    "3. スピン軌道は何個あります？  \n",
    "4. この分子をシミュレートするためにジョルダン・ウィグナー変換では何量子ビット必要ですか？\n",
    "5. 核反発エネルギーの値はいくつですか？\n",
    "\n",
    "解答は、このノートブックの最後にあります。\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# このラインの下にコードを書いてください \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# このラインの上にコードを書いてください "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 電子構造の問題\n",
    "\n",
    "次に、量子ビット（パウリ・ストリング）に変換する前に、フェルミ演算子のリストを生成する`ElectronicStructureProblem`を作ります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem\n",
    "problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "# 第二量子化演算子を作ります\n",
    "second_q_ops = problem.second_q_ops()\n",
    "\n",
    "# ハミルトニアン\n",
    "main_op = second_q_ops[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. 量子ビットへの変換\n",
    "シミュレーションで使用するマッピングを定義できます。\n",
    "別の変換も使うことができますが、シンプルな対応である`JordanWignerMapper`を使います：量子ビットを分子におけるスピン軌道として表現します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.mappers.second_quantization import ParityMapper, BravyiKitaevMapper, JordanWignerMapper\n",
    "from qiskit_nature.operators.second_quantization.qubit_converter import QubitConverter\n",
    "\n",
    "# 変換器と量子ビット変換を設定\n",
    "mapper_type = 'JordanWignerMapper'\n",
    "\n",
    "if mapper_type == 'ParityMapper':\n",
    "    mapper = ParityMapper()\n",
    "elif mapper_type == 'JordanWignerMapper':\n",
    "    mapper = JordanWignerMapper()\n",
    "elif mapper_type == 'BravyiKitaevMapper':\n",
    "    mapper = BravyiKitaevMapper()\n",
    "\n",
    "converter = QubitConverter(mapper=mapper, two_qubit_reduction=False)\n",
    "\n",
    "# フェルミオン演算子が量子ビット演算子に変換されます\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "qubit_op = converter.convert(main_op, num_particles=num_particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. 初期状態\n",
    "\n",
    "理論の節で説明したように、化学の問題における良い初期状態はHF状態(つまり$|\\Psi_{HF} \\rangle = |0101 \\rangle$)です。以下のように初期化できます："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.circuit.library import HartreeFock\n",
    "\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals\n",
    "init_state = HartreeFock(num_spin_orbitals, num_particles, converter)\n",
    "print(init_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5. 試行状態(Ansatz)\n",
    "\n",
    "最も重要な選択の一つが、基底状態を近似するために選ぶ量子回路です。\n",
    "回路として選択可能な量子回路ライブラリーの例を紹介します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.circuit.library import TwoLocal\n",
    "from qiskit_nature.circuit.library import UCCSD, PUCCD, SUCCD\n",
    "\n",
    "# ansatzを選ぶ\n",
    "ansatz_type = \"TwoLocal\"\n",
    "\n",
    "# q-UCC ansatzのためのパラメーター\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals\n",
    "\n",
    "# twolocalのための引数を設定\n",
    "if ansatz_type == \"TwoLocal\":\n",
    "    # 単一量子ビット回転は、すべての量子ビットに独立なパラメーターとして置かれます\n",
    "    rotation_blocks = ['ry', 'rz']\n",
    "    # エンタングルさせるゲート\n",
    "    entanglement_blocks = 'cx'\n",
    "    # 量子ビットをどのくらいエンタングルさせるか\n",
    "    entanglement = 'full'\n",
    "    # 独立なパラメーターを持った回転ブロックとエンタングルメントのブロックの繰り返し回数\n",
    "    repetitions = 3\n",
    "    # 最後の回転ブロック層をスキップする\n",
    "    skip_final_rotation_layer = True\n",
    "    ansatz = TwoLocal(qubit_op.num_qubits, rotation_blocks, entanglement_blocks, reps=repetitions, \n",
    "                      entanglement=entanglement, skip_final_rotation_layer=skip_final_rotation_layer)\n",
    "    # 初期状態を加える\n",
    "    ansatz.compose(init_state, front=True, inplace=True)\n",
    "elif ansatz_type == \"UCCSD\":\n",
    "    ansatz = UCCSD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"PUCCD\":\n",
    "    ansatz = PUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"SUCCD\":\n",
    "    ansatz = SUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"Custom\":    \n",
    "    # 回路の作り方の例    \n",
    "    from qiskit.circuit import Parameter, QuantumCircuit, QuantumRegister    \n",
    "    # 変分パラメーターの定義\n",
    "    theta = Parameter('a')\n",
    "    n = qubit_op.num_qubits\n",
    "    # 空の量子回路の作成\n",
    "    qc = QuantumCircuit(qubit_op.num_qubits)\n",
    "    qubit_label = 0\n",
    "    # アダマールゲートを置く\n",
    "    qc.h(qubit_label)\n",
    "    # CNOTを置く\n",
    "    for i in range(n-1):\n",
    "        qc.cx(i, i+1)\n",
    "    # 見た目のためのセパレーター\n",
    "    qc.barrier()\n",
    "    # 全量子ビットへrz回転\n",
    "    qc.rz(theta, range(n))\n",
    "    ansatz = qc\n",
    "    ansatz.compose(init_state, front=True, inplace=True)\n",
    "\n",
    "print(ansatz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6. バックエンド\n",
    "\n",
    "アルゴリズムを実行するシミュレーターまたはデバイスを特定するところです。\n",
    "このチャレンジでは`statevector_simulator`を使います。\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import Aer\n",
    "backend = Aer.get_backend('statevector_simulator')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7. オプティマイザー\n",
    "\n",
    "オプティマイザーは、試行回路のパラメーターの変化を導き、QPUで実行される測定の数が定義されるエネルギーの収束の評価のためにとても重要です。\n",
    "賢い選択をすることで、必要となるエネルギー評価の回数を劇的に激減できるでしょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA, SLSQP\n",
    "\n",
    "optimizer_type = 'COBYLA'\n",
    "\n",
    "# 各オプティマイザーのパラメーターを調整したいかもしれませんが\n",
    "# ここではデフォルト値を使います\n",
    "if optimizer_type == 'COBYLA':\n",
    "    optimizer = COBYLA(maxiter=1000)\n",
    "elif optimizer_type == 'L_BFGS_B':\n",
    "    optimizer = L_BFGS_B(maxfun=1000)\n",
    "elif optimizer_type == 'SPSA':\n",
    "    optimizer = SPSA(maxiter=500)\n",
    "elif optimizer_type == 'SLSQP':\n",
    "    optimizer = SLSQP(maxiter=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8. 厳密固有値ソルバー\n",
    "\n",
    "学習目的のために、問題をハミルトニアン行列の厳密な対角化によって解くことで、VQEがどこを目指しているかを知ることができます。\n",
    "もちろん、この行列の次元は分子軌道の数が増えるとともに指数関数的に増大するので、より大きな分子に対してこの方法を試すと、とても時間がかかることを確認できるでしょう。とても大きなシステムにおいては、その波動関数を保持するためのメモリーが足りなくなるでしょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.algorithms.ground_state_solvers.minimum_eigensolver_factories import NumPyMinimumEigensolverFactory\n",
    "from qiskit_nature.algorithms.ground_state_solvers import GroundStateEigensolver\n",
    "import numpy as np \n",
    "\n",
    "def exact_diagonalizer(problem, converter):\n",
    "    solver = NumPyMinimumEigensolverFactory()\n",
    "    calc = GroundStateEigensolver(converter, solver)\n",
    "    result = calc.solve(problem)\n",
    "    return result\n",
    "\n",
    "result_exact = exact_diagonalizer(problem, converter)\n",
    "exact_energy = np.real(result_exact.eigenenergies[0])\n",
    "print(\"Exact electronic energy\", exact_energy)\n",
    "print(result_exact)\n",
    "\n",
    "# よって目標とする電子のエネルギーは、-1.85336 Haです。\n",
    "# ご自分のVQEの結果を確認してください。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 9. VQEと試行状態のための初期パラメーター\n",
    "これで、VQEクラスをインポートし、アルゴリズムを実行できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms import VQE\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# リストのデータを出力して保存します\n",
    "def callback(eval_count, parameters, mean, std):\n",
    "    # 出力時に同じ行を上書きしてください\n",
    "    display(\"Evaluation: {}, Energy: {}, Std: {}\".format(eval_count, mean, std))\n",
    "    clear_output(wait=True)\n",
    "    counts.append(eval_count)\n",
    "    values.append(mean)\n",
    "    params.append(parameters)\n",
    "    deviation.append(std)\n",
    "\n",
    "counts = []\n",
    "values = []\n",
    "params = []\n",
    "deviation = []\n",
    "\n",
    "# ansatzの初期パラメーターを設定します\n",
    "# すべての参加者が似たような初期ポイントから始められるように\n",
    "# 固定の小さい変位を選びます\n",
    "try:\n",
    "    initial_point = [0.01] * len(ansatz.ordered_parameters)\n",
    "except:\n",
    "    initial_point = [0.01] * ansatz.num_parameters\n",
    "\n",
    "algorithm = VQE(ansatz,\n",
    "                optimizer=optimizer,\n",
    "                quantum_instance=backend,\n",
    "                callback=callback,\n",
    "                initial_point=initial_point)\n",
    "\n",
    "result = algorithm.compute_minimum_eigenvalue(qubit_op)\n",
    "\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 9. スコアリング関数\n",
    "\n",
    "あなたが選択した試行状態/オプティマイザーで、VQEシミュレーションがどれだけ良くできたかを判定する必要があります。\n",
    "そのために、以下のスコアリング関数を導入します：\n",
    "\n",
    "\n",
    "$$ score = N_{CNOT}$$\n",
    "\n",
    "\n",
    "ここで、$N_{CNOT}$はCNOTの数です。\n",
    "しかし、我々は、化学的精度$\\delta E_{chem} = 0.004$ Ha $= 4$ mHaに到達する必要があります。この値は、問題によっては到達が難しい場合があります。\n",
    "最も少ないCNOTの数でこの精度に到達した人がこのチャレンジの勝者です。スコアの値が小さい方が良いスコアです！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ディクショナリーに結果を保存します\n",
    "from qiskit.transpiler import PassManager\n",
    "from qiskit.transpiler.passes import Unroller\n",
    "\n",
    "# 作った回路をCNOTとUゲートにトランスパイルするUnroller\n",
    "pass_ = Unroller(['u', 'cx'])\n",
    "pm = PassManager(pass_)\n",
    "ansatz_tp = pm.run(ansatz)\n",
    "cnots = ansatz_tp.count_ops()['cx']\n",
    "score = cnots\n",
    "\n",
    "accuracy_threshold = 4.0 # mHaの単位\n",
    "energy = result.optimal_value\n",
    "\n",
    "if ansatz_type is \"TwoLocal\":\n",
    "    result_dict = {\n",
    "        'optimizer': optimizer.__class__.__name__,\n",
    "        'mapping': converter.mapper.__class__.__name__,\n",
    "        'ansatz': ansatz.__class__.__name__,\n",
    "        'rotation blocks': rotation_blocks,\n",
    "        'entanglement_blocks': entanglement_blocks,\n",
    "        'entanglement': entanglement,\n",
    "        'repetitions': repetitions,\n",
    "        'skip_final_rotation_layer': skip_final_rotation_layer,\n",
    "        'energy (Ha)': energy,\n",
    "        'error (mHa)': (energy-exact_energy)*1000,\n",
    "        'pass': (energy-exact_energy)*1000 <= accuracy_threshold,\n",
    "        '# of parameters': len(result.optimal_point),\n",
    "        'final parameters': result.optimal_point,\n",
    "        '# of evaluations': result.optimizer_evals,\n",
    "        'optimizer time': result.optimizer_time,\n",
    "        '# of qubits': int(qubit_op.num_qubits),\n",
    "        '# of CNOTs': cnots,\n",
    "        'score': score}\n",
    "else:\n",
    "    result_dict = {\n",
    "        'optimizer': optimizer.__class__.__name__,\n",
    "        'mapping': converter.mapper.__class__.__name__,\n",
    "        'ansatz': ansatz.__class__.__name__,\n",
    "        'rotation blocks': None,\n",
    "        'entanglement_blocks': None,\n",
    "        'entanglement': None,\n",
    "        'repetitions': None,\n",
    "        'skip_final_rotation_layer': None,\n",
    "        'energy (Ha)': energy,\n",
    "        'error (mHa)': (energy-exact_energy)*1000,\n",
    "        'pass': (energy-exact_energy)*1000 <= accuracy_threshold,\n",
    "        '# of parameters': len(result.optimal_point),\n",
    "        'final parameters': result.optimal_point,\n",
    "        '# of evaluations': result.optimizer_evals,\n",
    "        'optimizer time': result.optimizer_time,\n",
    "        '# of qubits': int(qubit_op.num_qubits),\n",
    "        '# of CNOTs': cnots,\n",
    "        'score': score}\n",
    "\n",
    "# 結果をプロットします\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.set_xlabel('Iterations')\n",
    "ax.set_ylabel('Energy')\n",
    "ax.grid()\n",
    "fig.text(0.7, 0.75, f'Energy: {result.optimal_value:.3f}\\nScore: {score:.0f}')\n",
    "plt.title(f\"{result_dict['optimizer']}-{result_dict['mapping']}\\n{result_dict['ansatz']}\")\n",
    "ax.plot(counts, values)\n",
    "ax.axhline(exact_energy, linestyle='--')\n",
    "fig_title = f\"\\\n",
    "{result_dict['optimizer']}-\\\n",
    "{result_dict['mapping']}-\\\n",
    "{result_dict['ansatz']}-\\\n",
    "Energy({result_dict['energy (Ha)']:.3f})-\\\n",
    "Score({result_dict['score']:.0f})\\\n",
    ".png\"\n",
    "fig.savefig(fig_title, dpi=300)\n",
    "\n",
    "# データを表示して保存します\n",
    "import pandas as pd\n",
    "import os.path\n",
    "filename = 'results_h2.csv'\n",
    "if os.path.isfile(filename):\n",
    "    result_df = pd.read_csv(filename)\n",
    "    result_df = result_df.append([result_dict])\n",
    "else:\n",
    "    result_df = pd.DataFrame.from_dict([result_dict])\n",
    "result_df.to_csv(filename)\n",
    "result_df[['optimizer','ansatz', '# of qubits', '# of parameters','rotation blocks', 'entanglement_blocks',\n",
    "    'entanglement', 'repetitions', 'error (mHa)', 'pass', 'score']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "<b>練習問題 2</b> \n",
    "    \n",
    "すべてのパラメーターを実験してから：\n",
    "\n",
    "1. ベストな(最も良いスコアの)経験的なansatz(試行回路)とオプティマイザーを見つけられましたか？(`TwoLocal` ansatzのパラメーターを修正することで)    \n",
    "2. ベストなq-UCC ansatzとオプティマイザーを見つけられましたか？(`UCCSD, PUCCD, SUCCD` ansatzの中から選択することで)   \n",
    "3. ansatzを定義するセルで、自分自身でゲートを置くことで`Custom` ansatzを修正して、`TwoLocal`の回路より良い回路を作れましたか？\n",
    "    \n",
    "各質問ごとに、`ansatz` オブジェクトを指定します。\n",
    "化学的精度$|E_{exact} - E_{VQE}| \\leq 0.004 $ Ha $= 4$ mHa に到達する必要があることを忘れないでください。\n",
    "\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# このラインの下にコードを書いてください \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# このラインの上にコードを書いてください "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "# パート2: チャレンジ -  LiH分子に対するVQE\n",
    "\n",
    "\n",
    "このパートでは、PySCFドライバーでSTG-3G基底を用いて、LiH分子をシミュレートします。\n",
    "\n",
    "</div>\n",
    "    \n",
    "<div class=\"alert alert-block alert-success\">\n",
    " \n",
    "<b>目標</b> \n",
    "    \n",
    "すべてのパラメーターを実験して、ベストな試行状態(ansatz)を見つけてください。好きなだけ創造力を発揮してください！\n",
    "\n",
    "すべての質問において、パート1と同じように`ansatz`オブジェクトを与えてください。最終スコアは、パート2の結果のみで判定されます。\n",
    "    \n",
    "</div>\n",
    "\n",
    "今回はシステムが大きくなったことに注意してください。このシステムに対して量子ビットが何個必要になるか、スピン軌道の数を得ることによって計算します。\n",
    "\n",
    "   \n",
    "\n",
    "### 問題サイズの削減\n",
    "\n",
    "シミュレーション実行時に量子ビット数を減らしたいかもしれません。\n",
    "- 化学に大きく貢献していないコア電子を凍結し、価電子のみを考慮することができます。Qiskitには、すでにこの機能が実装されています。よって、`qiskit_nature.transformers` において様々なtransformersを調べて、「フローズンコア近似」を実現するものを探してください。\n",
    "- 2量子ビットを削減するために、`two_qubit_reduction=True`とともに`ParityMapper`を使うことができます。\n",
    "- ハミルトニアンの対称性を調べることによって量子ビットの数を削減できる可能性があります。Qiskitの`Z2Symmetries`を使う方法を探してみてください。\n",
    "\n",
    "### 独自の試行状態(ansatz)\n",
    "\n",
    "[Grimsley *et al.*,2018](https://arxiv.org/abs/1812.11173v2), [H. L. Tang *et al.*,2019](https://arxiv.org/abs/1911.10205), [Rattew *et al.*,2019](https://arxiv.org/abs/1910.09694), [Tang *et al.*,2019](https://arxiv.org/abs/1911.10205)に提案されているアイディアを使ってみたいかもしれません。ベストな試行回路の生成する機械学習アルゴリズムも参考になるかもしれません。\n",
    "\n",
    "### シミュレーションの設定\n",
    "\n",
    "次に、ハートリー・フォック計算を実行させ、その後は、あなたにお任せします！\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "<b>注意事項</b> \n",
    "\n",
    "`driver`、`initial_point`、`initial_state` は与えられたものを使ってください。\n",
    "それ以外はQiskitに用意されているものを全て自由に使うことができます。    \n",
    "つまり、以下の初期点(すべてのパラメーターが0.01に設定されています)から始めてください：  \n",
    "    \n",
    "`initial_point = [0.01] * len(ansatz.ordered_parameters)`\n",
    "    または\n",
    "`initial_point = [0.01] * ansatz.num_parameters`\n",
    "\n",
    "また、初期状態はHF状態を使わなければいけません：\n",
    "    \n",
    "`init_state = HartreeFock(num_spin_orbitals, num_particles, converter)`\n",
    "    \n",
    "質問ごとに、`ansatz`オブジェクトを指定します。\n",
    "化学的な精度$|E_{exact} - E_{VQE}| \\leq 0.004 $ Ha$ = 4$ mHa に到達する必要があることを忘れないでください。\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers import PySCFDriver\n",
    "\n",
    "molecule = 'Li 0.0 0.0 0.0; H 0.0 0.0 1.5474'\n",
    "driver = PySCFDriver(atom=molecule)\n",
    "qmolecule = driver.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# このラインの下にコードを書いてください \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# このラインの上にコードを書いてください "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 以下のコードで答えを確認してください。\n",
    "from qc_grader import grade_ex5\n",
    "freeze_core = False # 核電子を凍結した場合はTrueにします。\n",
    "grade_ex5(ansatz,qubit_op,result,freeze_core)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 答えを提出します。何度でも再提出できます。\n",
    "from qc_grader import submit_ex5\n",
    "submit_ex5(ansatz,qubit_op,result,freeze_core)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# パート1の質問の解答"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "<b> 練習問題 1</b> \n",
    "    \n",
    "`qmolecule`の特性を調べて、以下の質問に答えましょう。\n",
    "    \n",
    "1. 分子の基本的特性を知る必要があります。与えられた系の電子は合計何個ですか？\n",
    "2. 分子軌道は何個あります？    \n",
    "3. スピン軌道は何個あります？  \n",
    "4. この分子をシミュレートするためにジョルダン・ウィグナー変換では何量子ビット必要ですか？\n",
    "5. 核反発エネルギーの値はいくつですか？\n",
    "\n",
    "    \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "<b>解答 </b> \n",
    "\n",
    "1. `n_el = qmolecule.num_alpha + qmolecule.num_beta`\n",
    "    \n",
    "2. `n_mo = qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "3. `n_so = 2 * qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "4. `n_q = 2* qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "5. `e_nn = qmolecule.nuclear_repulsion_energy`\n",
    "    \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "**Created by:** Igor Sokolov, Junye Huang, Rahul Pratap Singh\n",
    "\n",
    "**翻訳者:** 沼田祈史\n",
    "\n",
    "**Version:** 1.0.1"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}