tutorial/quantum_simulation/VQE_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Eigensolver\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "It is widely believed that one of the most promising applications of quantum computing in the near future is solving quantum chemistry problems [1-2]. **Variational Quantum Eigensolver** (VQE) is a strong proof to this possibility of studying quantum chemistry with **Noisy Intermediate-Scale Quantum** (NISQ) devices [1-4]. The core task is to solve the ground state of any molecular Hamiltonian $\\hat{H}$ by preparing a parametrized wave function ansatz $|\\Psi(\\boldsymbol\\theta)\\rangle$ on a quantum computer and adopt classical optimization methods (e.g. gradient descent) to adjust the parameters $\\boldsymbol\\theta$ to minimize the expectation value $\\langle \\Psi(\\boldsymbol\\theta)|\\hat{H}|\\Psi(\\boldsymbol\\theta)\\rangle$. This approach is based on the **Rayleigh-Ritz variational principle**. \n",
    "\n",
    "$$\n",
    "E_0 = \\min_{\\boldsymbol\\theta} \\langle \\Psi(\\boldsymbol\\theta)|\\hat{H}|\\Psi(\\boldsymbol\\theta)\\rangle.\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "where $E_0$ denotes the ground state energy. Numerically, it can be understood as finding the smallest eigenvalue $\\lambda_{\\min}$ of a **discretized** Hamiltonian $H$ (hermitian matrix) and its corresponding eigenvector $|\\Psi_0\\rangle$. How such a discretization can be done on a classical computer belongs to the art of quantum chemistry and is far beyond the scope of this tutorial. We will discuss this part with a few words in the background section. In general, such a Hamiltonian $H$ is expressed as a weighted sum of Pauli spin operators $\\{X,Y,Z\\}$ (native to quantum devices) such that this information can be processed on a quantum computer.\n",
    "\n",
    "$$\n",
    "H = \\sum_k c_k ~ \\bigg( \\bigotimes_{j=0}^{M-1} \\sigma_j^{(k)} \\bigg),\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "where $\\sigma_j^{(k)} \\in \\{I,X,Y,Z\\}$ and $M$ stands for qubit number. We refer this form of Hamiltonian as **Pauli strings**. For example, \n",
    "\n",
    "$$\n",
    "H= 0.12~Y_0 \\otimes I_1-0.04~X_0\\otimes Z_1.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "In the next section, we will provide a brief review on the electronic structure problem which essentially tells us how to calculate the Hamiltonian $H$. For those who are already familiar with this topic or only interested in how to implement VQE on Paddle Quantum, please skip this part and jump into the illustrative example of hydrogen molecule $H_2$.\n",
    "\n",
    "## Background: the electronic structure problem\n",
    "\n",
    "In this section, we focus on one of the fundamental problems in quantum chemistry --  **the electronic structure problem**. To be more specific, we are interested in the low lying energy eigenstates of any given molecule. These knowledge could help predict reaction rates and location of stable structures [5]. Suppose a molecule consists of $N_n$ nuclei and $N_e$ electrons, the first quantized (canonical quantization) Hamiltonian operator $\\hat{H}_{mol}$ describing the total energy of this molecular system can be written as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{H}_{\\text{mol}} & = -\\sum_{i}\\frac{\\nabla_{R_i}^2}{2M_i} - \\sum_{i} \\frac{\\nabla_{r_i}^2}{2} -\\sum_{i,j}\\frac{Z_i}{\\lvert R_i - r_j\\lvert} + \\sum_{i,j>i}\\frac{Z_iZ_j}{\\lvert R_i - R_j\\lvert} + \\sum_{i, j>i}\\frac{1}{\\lvert r_i - r_j\\lvert}, \n",
    "\\tag{4}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $R_i, M_i,$ and $Z_i$ denote the position, mass and atomic number (the number of protons) of the $i^{th}$ nucleus, and the positions of electrons are $r_i$. The first two sums describe the kinetic energy of nuclei and electrons, respectively. The third sum describes the attractive Coulomb interaction between the positively charged nuclei and the negatively charged electrons. The last two terms represent the repulsive nuclei-nuclei and electron-electron interactions.  Here, the molecular Hamiltonian $\\hat{H}_\\text{mol}$ is already in atomic units of energy, **Hartree**. $1$ Hartree is $[\\hbar^2/(m_ee^2a_0^2)] = 27.2$ eV or $630$ kcal/mol, where $m_e, e,$ and $a_0$ stand for the mass of electron, charge of electron, and Bohr radius. \n",
    "\n",
    "\n",
    "**Note:** The spin-orbit interaction and hyperfine interaction are not considered in this picture. They can be treated as perturbations if necessary. \n",
    "\n",
    "### Born-Oppenheimer approximation\n",
    "\n",
    "Since the nuclei are much heavier than electrons, the electrons will move much faster than the nuclei. It is reasonable to treat the positions of nuclei as fixed, $R_i =$ constants. This is known as the Born-Oppenheimer approximation by decoupling the behavior of nuclei and electrons in time scale. Consequently, the kinetic energy term of nuclei will disappear and the nuclei-nuclei repulsive interaction term can be viewed as an energy shift (independent of electron positions $r_i$). We could derive the electronic Hamiltonian $\\hat{H}_{\\text{electron}}$ as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{H}_{\\text{electron}} & =  - \\sum_{i} \\frac{\\nabla_{r_i}^2}{2} -\\sum_{i,j}\\frac{Z_i}{\\lvert R_i - r_j\\lvert} + \\sum_{i, j>i}\\frac{1}{\\lvert r_i - r_j\\lvert} \n",
    "\\tag{5},\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The energy levels of the electrons in the molecule can be found by solving the time independent Schrödinger equation\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\text{electron}} |\\Psi_n \\rangle = E_n |\\Psi_n \\rangle,\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "where $n$ stands for the energy level. Notice the electron repulsion terms scale as $N_e(N_e-1)/2$ which means for the Oxygen molecule $O_2$ carrying 16 electrons there will be 120 electron repulsion terms in total! In general, this problem cannot be solved analytically. As Dirac concluded in [Quantum mechanics of many-electron systems](https://royalsocietypublishing.org/doi/10.1098/rspa.1929.0094) [6],\n",
    "\n",
    "> *The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble.*                                                                                                                                                                  \n",
    ">\n",
    "> ​\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t-- Paul Dirac (1929)\n",
    "\n",
    "A straightforward numerical approach is discretizing the infinite-dimensional Hilbert space into equidistant grid points where linear algebra guides the whole calculation. Suppose each axis of space is discretized into $k$ points, the $N$-electron (drop the subscript e for simplicity) wave function can be written as [2]\n",
    "\n",
    "$$\n",
    "|\\Psi \\rangle = \\sum_{\\mathbf{x_1}, \\ldots, \\mathbf{x_N}} \\psi(\\mathbf{x_1}, \\ldots, \\mathbf{x_N}) \\mathcal{A}(|\\mathbf{x_1}, \\ldots, \\mathbf{x_N}\\rangle).\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "where coordinate $|\\mathbf{x_j}\\rangle = |r_j\\rangle |\\sigma_j\\rangle$ records the spatial location and spin of the $j^{th}$ electron, $|r_j\\rangle  = |x_j,y_j,z_j\\rangle$ for $j\\in \\{1,2,\\cdots,N\\}$, $x_j,y_j,z_j \\in \\{0,1,\\cdots,k-1\\}$ and  $\\sigma_j \\in \\{\\downarrow,\\uparrow\\}$ for spin down or up.  There will be $k^{3N}\\times 2^{N}$ complex amplitudes in total. Here, $\\mathcal{A}$ denotes anti-symmetrization and a consequence of the Pauli exclusion principle (electrons are fermion), and $\\psi(\\mathbf{x_1}, \\mathbf{x_2}, \\ldots, \\mathbf{x_N})=\\langle\\mathbf{x_1}, \\mathbf{x_2}, \\ldots, \\mathbf{x_N}|\\Psi\\rangle$. One can see that storing such a wave function already requires **exponentially growing memory** with respect to the number of electrons $N$. This would make classical simulation methods based on this naive numerical approach intractable for systems size larger than few tens of electrons. Now, the question becomes can we prepare such a wave function $|\\Psi\\rangle$ directly on a quantum computer and measure the expectation value $E_0$? In the next section, let's take the simplest molecular system -- hydrogen molecule $H_2$ as a concrete example.\n",
    "\n",
    "\n",
    "\n",
    "**Note:** A detailed review on quantum chemistry and existing classical computational methods are far beyond the scope of this tutorial, we refer the enthusiastic readers to the standard textbooks *'Molecular Electronic-Structure Theory'* [5] by Helgaker and *'Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory'* [7] by Szabo & Ostlund. To bridge to knowledge gap between quantum chemistry and quantum computing, please check the following review papers [Quantum computational chemistry](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.92.015003) [2] and [Quantum Chemistry in the Age of Quantum Computing](https://pubs.acs.org/doi/10.1021/acs.chemrev.8b00803) [1].\n",
    "\n",
    "**Note:** For energy calculation, it is desired to reach the **chemical accuracy** of $1.6\\times10^{-3}$ Hartree or 1 kcal/mol . "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ground state of the hydrogen molecule $H_2$\n",
    "\n",
    "### Building electronic Hamiltonian\n",
    "\n",
    "First of all, let us import the necessary libraries and packages.`qchem` module in Paddle Quantum is developed basing on `psi4` and `openfermion`, so you need to install these two packages before executing the following codes. We strongly encourage you to read the tutorial [Building molecular Hamiltonian](./BuildingMolecule_EN.ipynb) first, which introduces how to utilize our quantum chemistry toolkit `qchem`.\n",
    "\n",
    "**Note: As to the environment setting, please refer to [README.md](https://github.com/PaddlePaddle/Quantum/blob/master/README.md).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:44.970178Z",
     "start_time": "2021-04-30T09:14:40.895128Z"
    }
   },
   "outputs": [],
   "source": [
    "import paddle\n",
    "import paddle_quantum.qchem as qchem\n",
    "from paddle_quantum.utils import Hamiltonian\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy\n",
    "from numpy import pi as PI\n",
    "from numpy import savez, zeros\n",
    "\n",
    "# Eliminate noisy python warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To analyze specific molecules, we need several key information such as **geometry**, **basis set** (such as STO-3G), **multiplicity**, and **charge** to model the molecule and achieve more information about the molecule, such as one-body integrations, two-body integrations, and others. Next, through our built-in quantum chemistry toolkit, we could set a `Hamiltonian` class to carry molecule Hamiltonian information, which may simplify the following calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:44.982005Z",
     "start_time": "2021-04-30T09:14:44.975892Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FCI energy for H2_sto-3g_singlet (2 electrons) is -1.137283834485513.\n",
      "\n",
      "The generated h2 Hamiltonian is \n",
      " -0.09706626861762556 I\n",
      "-0.04530261550868938 X0, X1, Y2, Y3\n",
      "0.04530261550868938 X0, Y1, Y2, X3\n",
      "0.04530261550868938 Y0, X1, X2, Y3\n",
      "-0.04530261550868938 Y0, Y1, X2, X3\n",
      "0.1714128263940239 Z0\n",
      "0.16868898168693286 Z0, Z1\n",
      "0.12062523481381837 Z0, Z2\n",
      "0.16592785032250773 Z0, Z3\n",
      "0.17141282639402394 Z1\n",
      "0.16592785032250773 Z1, Z2\n",
      "0.12062523481381837 Z1, Z3\n",
      "-0.2234315367466399 Z2\n",
      "0.17441287610651626 Z2, Z3\n",
      "-0.2234315367466399 Z3\n"
     ]
    }
   ],
   "source": [
    "geo = qchem.geometry(structure=[['H', [-0., 0., 0.0]], ['H', [-0., 0., 0.74]]])\n",
    "# geo = qchem.geometry(file='h2.xyz')\n",
    "\n",
    "# Save molecule information in to variable molecule, including one-body integrations, one-body integrations, molecular and Hamiltonian\n",
    "molecule = qchem.get_molecular_data(\n",
    "    geometry=geo,\n",
    "    basis='sto-3g',\n",
    "    charge=0,\n",
    "    multiplicity=1,\n",
    "    method=\"fci\",\n",
    "    if_save=True,\n",
    "    if_print=True\n",
    ")\n",
    "# Recall Hamiltonian\n",
    "molecular_hamiltonian = qchem.spin_hamiltonian(molecule=molecule,\n",
    "                                               filename=None, \n",
    "                                               multiplicity=1, \n",
    "                                               mapping_method = 'jordan_wigner',)\n",
    "# Print results\n",
    "print(\"\\nThe generated h2 Hamiltonian is \\n\", molecular_hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** This Hamiltonian is generated with an interatomic distance of $d = 74$ pm. \n",
    "\n",
    "In addition to inputting molecular geometry directly, inputting molecular geometry file (.xyz file) is also allowed. For more information about quantum chemistry toolkit, please refer to the tutorial [Building molecular Hamiltonian](./BuildingMolecule_EN.ipynb). If you need to test the geometric configuration of more molecules, please check out this [database](http://smart.sns.it/molecules/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building QNN and trial wave function\n",
    "\n",
    "To implement VQE, we firstly need to design a quantum neural network QNN to prepare the wave function ansatz $|\\Psi(\\boldsymbol\\theta)\\rangle$. Here, we provide a 4-qubit quantum circuit template with a depth of $D$ blocks. The dotted frame in the figure below denotes a single block:\n",
    "\n",
    "\n",
    "![Utheta.jpg](https://release-data.cdn.bcebos.com/PIC%2FUtheta.jpg)\n",
    "\n",
    "Next, we use the `UAnsatz` class and the built-in `real_entangled_layer(theta, D)` circuit template in Paddle Quantum to realize this QNN.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:50.083041Z",
     "start_time": "2021-04-30T09:14:50.062255Z"
    }
   },
   "outputs": [],
   "source": [
    "def U_theta(theta, Hamiltonian, N, D):\n",
    "    \"\"\"\n",
    "    Quantum Neural Network\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialize the quantum neural network according to the number of qubits N\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # Built-in {R_y + CNOT} circuit template\n",
    "    cir.real_entangled_layer(theta[:D], D)\n",
    "    \n",
    "    # Lay R_y gates in the last row\n",
    "    for i in range(N):\n",
    "        cir.ry(theta=theta[D][i][0], which_qubit=i)\n",
    "        \n",
    "    # The quantum neural network acts on the default initial state |0000>\n",
    "    fin_state = cir.run_state_vector()\n",
    "    \n",
    "    # Calculate the expectated value of the given Hamiltonian\n",
    "    expectation_val = cir.expecval(Hamiltonian)\n",
    "\n",
    "    return expectation_val, cir, fin_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the loss function and model\n",
    "\n",
    "Now that we have the target Hamiltonian and QNN, we will further define the training model and loss function. By applying the QNN $U(\\theta)$ on the initial state $|0..0\\rangle$, we get the output state $|\\psi(\\boldsymbol\\theta)\\rangle $. Then, the loss function to be minimized is the expectation value, \n",
    "\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol\\theta}  \\mathcal{L}(\\boldsymbol \\theta) = \\min_{\\boldsymbol\\theta} \\langle \\Psi(\\boldsymbol\\theta)|H |\\Psi(\\boldsymbol\\theta)\\rangle\n",
    "= \\min_{\\boldsymbol\\theta} \\sum_k c_k~\\langle \\Psi(\\boldsymbol\\theta)| \\bigotimes_j \\sigma_j^{(k)}|\\Psi(\\boldsymbol\\theta)\\rangle.\n",
    "\\tag{8}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:50.183996Z",
     "start_time": "2021-04-30T09:14:50.167892Z"
    }
   },
   "outputs": [],
   "source": [
    "class StateNet(paddle.nn.Layer):\n",
    "    \"\"\"\n",
    "    Construct the model net\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, shape, dtype=\"float64\"):\n",
    "        super(StateNet, self).__init__()\n",
    "        \n",
    "        # Initialize the list of theta parameters, filling the initial values with a uniform distribution of [0, 2* PI]  \n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "    # Define loss function and forward propagation mechanism\n",
    "    def forward(self, N, D):\n",
    "        \n",
    "        # Calculate loss function/expected value\n",
    "        loss, cir, fin_state = U_theta(self.theta, molecular_hamiltonian.pauli_str, N, D)\n",
    "\n",
    "        return loss, cir, fin_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hyper-parameters\n",
    "\n",
    "Before training the QNN, we also need to set some training hyper-parameters, mainly the learning rate (LR), the number of iterations (ITR), and the depth (D) of repeated blocks. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:50.222465Z",
     "start_time": "2021-04-30T09:14:50.187093Z"
    }
   },
   "outputs": [],
   "source": [
    "ITR = 80  # Set the number of optimization iterations\n",
    "LR = 0.4   # Set the learning rate\n",
    "D = 2      # Set the depth of the repetitive calculation module in QNN\n",
    "N = molecular_hamiltonian.n_qubits # Set number of qubits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training\n",
    "\n",
    "After all the training model parameters are set, we convert the data into Tensor in the Paddle, and then train the quantum neural network. The results of the training process is stored in the summary_data file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:15:52.165788Z",
     "start_time": "2021-04-30T09:15:29.625076Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 20 loss: -1.0510\n",
      "iter: 20 Ground state energy: -1.0510 Ha\n",
      "iter: 40 loss: -1.1298\n",
      "iter: 40 Ground state energy: -1.1298 Ha\n",
      "iter: 60 loss: -1.1361\n",
      "iter: 60 Ground state energy: -1.1361 Ha\n",
      "iter: 80 loss: -1.1371\n",
      "iter: 80 Ground state energy: -1.1371 Ha\n",
      "\n",
      "Circuit after training:\n",
      "--Ry(4.710)----*--------------x----Ry(1.569)----*--------------x----Ry(-0.01)--\n",
      "               |              |                 |              |               \n",
      "--Ry(1.507)----x----*---------|----Ry(1.566)----x----*---------|----Ry(4.392)--\n",
      "                    |         |                      |         |               \n",
      "--Ry(5.984)---------x----*----|----Ry(4.476)---------x----*----|----Ry(1.692)--\n",
      "                         |    |                           |    |               \n",
      "--Ry(0.132)--------------x----*----Ry(7.773)--------------x----*----Ry(3.194)--\n",
      "                                                                               \n"
     ]
    }
   ],
   "source": [
    "# Determine the parameter dimensions of the network \n",
    "net = StateNet(shape=[D + 1, N, 1])\n",
    "\n",
    "# In general, we use the Adam optimizer to obtain relatively good convergence,\n",
    "# Of course you can change it to SGD or RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# Record optimization results\n",
    "summary_iter, summary_loss = [], []\n",
    "\n",
    "# Optimize iterations\n",
    "for itr in range(1, ITR + 1):\n",
    "\n",
    "    # Forward propagation calculates the loss function\n",
    "    loss, cir, fin_state = net(N, D)\n",
    "\n",
    "    # Back propagation minimizes the loss function\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # Update optimization results\n",
    "    summary_loss.append(loss.numpy())\n",
    "    summary_iter.append(itr)\n",
    "\n",
    "    # Print results\n",
    "    if itr % 20 == 0:\n",
    "        print(\"iter:\", itr, \"loss:\", \"%.4f\" % loss.numpy())\n",
    "        print(\"iter:\", itr, \"Ground state energy:\", \"%.4f Ha\" \n",
    "                                            % loss.numpy())\n",
    "    if itr == ITR:\n",
    "        print(\"\\nCircuit after training:\") \n",
    "        print(cir)\n",
    "\n",
    "# Save the training results in the Output folder\n",
    "os.makedirs(\"output\", exist_ok=True)\n",
    "savez(\"./output/summary_data\", iter = summary_iter, \n",
    "                               energy=summary_loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmarking\n",
    "We have now completed the training of the quantum neural network, and the estimated value of the ground state energy obtained is $E_0 \\approx -1.137 $ Hartree, we compare it with the theoretical value $E_0 = -1.1371$ to benchmark our model. The estimation obtained with VQE agree with a full configuration-interaction (FCI) calculation within chemical accuracy $\\varepsilon = 1.6 \\times 10^{-3}$ Hartree.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:15:18.096944Z",
     "start_time": "2021-04-30T09:15:17.481250Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyKUlEQVR4nO3deXxU1fn48c9DQgKEhF1AEBIURAJJIGETQdkUV1ChaLHiglYrKrVa8av9ldbWtda61X2tWLQoAoqCgCiuGJB9VQgSZQ87hkB4fn+cmziEySSQ5c4kz/v1mtfce+fMvU9mJvPMPeeec0RVMcYYY4pTw+8AjDHGhDdLFMYYY0KyRGGMMSYkSxTGGGNCskRhjDEmpGi/AyhvjRs31sTERL/DMMaYiDJ//vxtqtok2GNVLlEkJiaSmZnpdxjGGBNRRGR9cY9Z1ZMxxpiQLFEYY4wJyRKFMcaYkKpcG4Ux1dXBgwfJzs4mNzfX71BMGKtVqxYtW7akZs2apX6OJQpjqojs7Gzi4+NJTExERPwOx4QhVWX79u1kZ2eTlJRU6udZ1ZMxVURubi6NGjWyJGGKJSI0atTomM86LVEYU4VYkjAlOZ7PiCWKAvv2wRtvwOrVfkdijDFhxRJFgcOH4b//hRUr/I7EGGPCiiWKAnXrggjs2eN3JMYYE1YsURQQgfh42L3b70iMiXj9+vXj0KFDIcv8/PPPnHnmmeTn5wOQn5/PrbfeSnJyMp06dWLt2rXk5eXRp0+fI/Z11llnkZWVBcCzzz7LjTfeeMR+O3bsyAqvZiCwbHnHsmPHDi6++OJSvR6RzhJFoPh4O6MwpoyWLVtGo0aNiI4OffX9Sy+9xCWXXEJUVBQA999/P23atGHZsmXccsst/Pvf/yYmJob+/fvz5ptvBt3HkiVL6NKlS+F6bm4uWVlZtGvX7phiPp5YGjRoQE5ODtu3bz+mY0Ui60cRKCHBzihM1fD887B2bfnus00buO66EotNnjyZIUOGANC5c2c++OADnnzySU455RSSkpJ4+umnmTBhAuPHj+eNN94AYN++fUyaNIn58+cDkJSUxPvvvw/AkCFDuOuuuxgxYsRRx1q8eDFXX3114fqSJUto165d4Rd+oIqI5fzzz2fq1KlcddVVpXkFI5YlikDx8bB1q99RGBPRpk2bxnvvvcehQ4fIycmhWbNmLFq0iKFDh/Lpp5+SmppKXl4ea9eupWBKgJkzZ7JhwwbS0tIAyMnJYcCAAYCrSvrmm2+CHmvZsmVccsklhZd87t27lwsuuOCochUVy+DBg7nzzjstUVQrCQnl/yvMGD+U4pd/Rdi/fz95eXnUr1+fpUuX0r59ewCWL19Ohw4deOKJJ7jkkkvYtm0b9evXL3zewoUL+etf/8oNN9wAwKhRo0hJSQEgKiqKmJgY9uzZQ3x8fOFzNmzYQJMmTVi5cmXhttGjRwftcbxy5coKieXUU09l1apV5fDKhTdrowhkjdnGlEmdOnUQEfbu3cuqVas49dRTycnJoW7dusTExJCZmUnXrl2pXbv2Eb2Dd+zYQZ06dQD363/GjBlceOGFhY8fOHCAWrVqHXGsJUuWkJycfMS25cuXF36pB6qoWNavX39MQ2FEKksUgRISIC8PDhzwOxJjItY555zDhx9+SExMDCtXriQzM5PU1FRef/11EhMTOeGEE2jQoAH5+fmFX9Dt2rXjq6++AuDRRx/l/PPPL/wC3r59O40bNz5qELvFixfToUOHI7YtW7aMTp06HRVTRcUyefJkBg8eXF4vXdiyRBGo4LTWrnwy5rgNHjyYd999l0GDBtG+fXtGjBjBnDlzyMzM5LXXXissd/bZZ/PZZ58BcPnll7NgwQJOOeUUFi9ezD//+c/Cch9//DHnn3/+UcdZsmTJEYkiJycHVaVZs2ZHla2oWKZOnVotEgWqWqVu6enpetw+/1z1ggtUv//++PdhjE+WL1/udwiFOnXqpAcPHlRV1auuukpnzJhxVJn58+frFVdcUeK+Lr74Yl21alXh+plnnqnr1q0rVRxFy5ZnLDk5Odq7d+9SxRFugn1WgEwt5nvVzigCJSS4ezujMKZMFi9eXNiPYvHixUHbDbp06ULfvn0LO7kFk5eXx5AhQ465X0SouMorlgYNGvDpp5+WS1zhzq56ClSQKKxB25hyU9AfIZhrrrkm5HNjYmK48sorj9h21VVXHXGVUihFy5Z3LNWFJYpA1kZhTNg7lj4LVb1/Q2XxtepJRAaJyCoR+U5ExgZ5PFZE3vQe/1pEEis0oIJEYWcUxhhTyLdEISJRwFPAuUAH4HIR6VCk2LXADlU9BXgUeLBCg4qOhjp17IzCGGMC+HlG0Q34TlXXqmoeMAEoep3ZYOBVb3ki0F8qegov63RnjDFH8DNRtAA2BKxne9uCllHVQ8AuoFHRHYnI9SKSKSKZW8s6VlNCgp1RGGNMgCpxeayqPqeqGaqa0aRJk7LtzM4ojDHmCH4mih+BkwLWW3rbgpYRkWigHlCxg7/bnBTGHLfNmzfz61//mjZt2pCenk7Pnj2ZNGlSpcaQlZVFx44dS11+zpw5fPHFF+VWriryM1F8A7QVkSQRiQEuA6YUKTMFGOktDwVmez0IK45VPRlzXFSVIUOG0KdPH9auXcv8+fOZMGEC2dnZR5Utafa7yhSJiaKyXz/fEoXX5jAamA6sAN5S1WUi8lcRucgr9iLQSES+A24DjrqEttwlJMC+fRBGH2RjIsHs2bOJiYkpHJ4boHXr1tx8880AvPLKK1x00UX069eP/v37k5OTw5AhQ0hJSaFHjx4sXrwYgHHjxvGPf/yjcB8dO3YkKyuLrKwsTjvtNK677jqSk5M5++yz+fnnnwHXkS41NZXU1FSeeuqpYmN8/PHH6dChAykpKVx22WVkZWXxzDPP8Oijj5KWlsbcuXOZOnUq3bt3p3PnzgwYMIDNmzcHLbd161YuvfRSunbtSteuXfn888+POl5+fj533HEHXbt2JSUlhWeffRZwSeess85i6NChhWNQFfwGnj9/PmeeeSbp6emcc845bNy4EXDTuo4ZM4aMjAwee+wxvvnmG1JSUkhLS+OOO+4oPIvq06cPCxcuLIzhjDPOYNGiRcf8fh6huLE9IvVWprGeVFXfe8+N97RjR9n2Y0wlKzp+z9ixqjNnuuWDB9367NluPTfXrX/6qVvfu9etf/65W9+1y61//bVbz8kp+fiPPfaYjhkzptjHX375ZW3RooVu375dVVVHjx6t48aNU1XVWbNmaWpqqqqq/vnPf9aHH3648HnJycm6bt06XbdunUZFRem3336rqqrDhg3T//znP6rqxpb65JNPVFX19ttv1+Tk5KAxNG/eXHNzc1VVdYf3P170eDk5OXr48GFVVX3++ef1tttuC1ru8ssv17lz56qq6vr167V9+/ZHHe/ZZ5/Ve++9V1VVc3NzNT09XdeuXasff/yxJiQk6IYNGzQ/P1979Oihc+fO1by8PO3Zs6du2bJFVVUnTJigV199taq6catuvPHGI16XL774QlVV77zzzsK/+ZVXXtFbb71VVVVXrVqlwb4Tj3WsJ+uZXVRg7+xSDhNgjDnaTTfdxGeffUZMTEzhrHADBw6kYcOGAHz22We8/fbbAPTr14/t27ezu4QLSZKSkgpnnktPTycrK4udO3eyc+dO+vTpA8BvfvMbPvjgg6DPT0lJYcSIEQwZMqRwutaisrOzGT58OBs3biQvL6/Y+SZmzpzJ8uXLC9d3797N3r17qVu3buG2GTNmsHjxYiZOnAjArl27WLNmDTExMXTr1o2WLVsCkJaWRlZWVuGETwMHDgTcGUnz5s0L9zd8+HAAdu7cyZ49e+jZsycAv/71r3nvvfcAGDZsGPfeey8PP/wwL730Urn0TrdEUZSN92SqiPvv/2U5OvrI9djYI9fj4o5cT0g4cr1Bg5KPl5ycXPjFD/DUU0+xbds2MjIyAo4TV+J+oqOjOXz4cOF64KRCsbGxhctRUVGFVU/Fufrqq/n222858cQTmTZtGu+//z6ffvopU6dO5e9//ztLliw56jk333wzt912GxdddBFz5sxh3LhxQfd9+PBhvvrqq6MmVAqkqjzxxBOcc845R2yfM2fOUX/LoUOHUFWSk5P58ssvg+6vNK9fnTp1GDhwIJMnT+att94KOb5VaVWJy2PLlY0ga8xx6devH7m5uTz99NOF2/bv319s+d69ezN+/HjAfXE2btyYhIQEEhMTWbBgAQALFixg3bp1IY9bv3596tevXzifRME+AV5++WUWLlzItGnTOHz4MBs2bKBv3748+OCD7Nq1i7179xIfH8+egP/3Xbt20aKF69L16quvFm4vWu7ss8/miSeeKFwPbBcocM455/D0009z8OBBAFavXs2+ffuK/VtOPfVUtm7dWpgoDh48yLJly4L+zfHx8Xz99dcATJgw4YjHR40axS233ELXrl1pUJosXwJLFEXZeE/GHBcR4d133+WTTz4hKSmJbt26MXLkSB58MPjIO+PGjWP+/PmkpKQwduzYwi/lSy+9lJycHJKTk3nyySdLNcT4yy+/zE033URaWlpho3BR+fn5XHHFFXTq1InOnTtzyy23UL9+fS688EImTZpU2Eg9btw4hg0bRnp6Oo0bNy58ftFyjz/+OJmZmaSkpNChQweeeeaZo445atQoOnToQJcuXejYsSO//e1vQ16xFBMTw8SJE7nzzjtJTU0lLS2t2CutXnzxRa677jrS0tLYt28f9erVK3wsPT2dhIQErr766hJfu9KQ4l7USJWRkaGZmZnHv4PcXBg2DEaOhKFDyy8wYyrYihUrOO200/wOw1SSwPaQBx54gI0bN/LYY48B8NNPP3HWWWexcuVKatQ4+nwg2GdFROarasZRhbEziqPFxkLNmlb1ZIwJa++//z5paWl07NiRuXPncs899wDw2muv0b17d/7+978HTRLHwxqzixKx3tnGmLA3fPjwwqugAl155ZXlPsGSnVEEk5BgbRQmIlW1qmRT/o7nM2KJIhhLFCYC1apVi+3bt1uyMMVSVbZv3x7ykt5grOopmPh4WL/e7yiMOSYtW7YkOzubMg+1b6q0WrVqFXb0Ky1LFMFYG4WJQDVr1iy2F7ExZWFVT8EUjCBrp/DGGGOJIqiEBDh8GEL0KjXGmOrCEkUw1jvbGGMKWaIIxhKFMcYUskQRjA0MaIwxhSxRBGOJwhhjClmiCMaqnowxppAlimDi4tyYT5YojDHGEkVQNjCgMcYUskRRHEsUxhgDWKIong0MaIwxgE+JQkQaishHIrLGuz9qUlcRSRORL0VkmYgsFpGjB16vSJYojDEG8O+MYiwwS1XbArO89aL2A1eqajIwCPiXiNSvtAit6skYYwD/EsVg4FVv+VVgSNECqrpaVdd4yz8BW4AmlRWgJQpjjHH8ShRNVXWjt7wJaBqqsIh0A2KA7ys6sEIJCZCXBwcOVNohjTEmHFXYfBQiMhNoFuShuwNXVFVFpNjxvEWkOfAfYKSqHi6mzPXA9QCtWrU67piPENg7Oza2fPZpjDERqMIShaoOKO4xEdksIs1VdaOXCLYUUy4BeB+4W1W/CnGs54DnADIyMspnEonA3tmNG5fLLo0xJhL5VfU0BRjpLY8EJhctICIxwCTgNVWdWImxOQVnFDt3VvqhjTEmnPiVKB4ABorIGmCAt46IZIjIC16ZXwF9gKtEZKF3S6u0CAvmlN2wodIOaYwx4ciXObNVdTvQP8j2TGCUt/w68Holh/aLevWgfn1Yt863EIwxJhxYz+xQEhNh/Xq/ozDGGF9ZogglMRF++MHNn22MMdWUJYpQEhNdX4qNG0ssaowxVZUlilASE929tVMYY6oxSxShnHSSm5vC2imMMdWYJYpQYmLcZbJ2RmGMqcYsUZSkdWvIyvI7CmOM8Y0lipIkJcHmzfDzz35HYowxvrBEUZLWrd29tVMYY6opSxQlSUpy91b9ZIyppixRlKRJE6hd2xKFMabaskRREhEbysMYU61ZoiiNxER3iayWz1QXxhgTSSxRlEZiIuzbB9u3+x2JMcZUOksUpVEwlIe1UxhjqiFLFKVRcImsJQpjTDVkiaI04uLc1U+WKIwx1ZAlitJKTLREYYyplkJOhSoitYALgN7AicDPwFLgfVVdVvHhhZGkJFiwAA4dgmhfZpA1xhhfFHtGISJ/AT4HegJfA88CbwGHgAdE5CMRSamUKMPBySdDfj6sWOF3JMYYU6lC/TSep6p/Luaxf4rICUCrCogpPHXp4oYd//xz6NTJ72iMMabSFHtGoarvh3qiqm5R1czyDylM1aoFGRnwxRfW8c4YU62U2JgtIk1E5B8iMk1EZhfcKiO4sNOrF+zYAcuX+x2JMcZUmtJc9TQeWAEkAX8BsoBvynJQEWnotXGs8e4bhCibICLZIvJkWY5ZLrp2hZo1XfWTMcZUE6VJFI1U9UXgoKp+oqrXAP3KeNyxwCxVbQvM8taLcy/waRmPVz5q14b0dKt+MsZUK6VJFAe9+40icr6IdAYalvG4g4FXveVXgSHBColIOtAUmFHG45WfXr3cmE8rV/odiTHGVIrSJIq/iUg94A/A7cALwO/LeNymqrrRW96ESwZHEJEawCPeMUMSketFJFNEMrdu3VrG0ErQtavrR2HVT8aYaqLEnmOq+p63uAvoW9odi8hMoFmQh+4usn8VkWD1OL8DpqlqtoiUFONzwHMAGRkZFVsnFBfnLpX9/HO49lo3X4UxxlRhxSYKEXkCKPZLV1VvCbVjVR0QYt+bRaS5qm4UkebAliDFegK9ReR3QF0gRkT2qmqo9ozK0asXzJsHq1fDqaf6HY0xxlSoUGcUgX0k/gIU1/nueEwBRgIPePeTixZQ1REFyyJyFZARFkkCoHv3X6qfLFEYY6q4YhOFqhY0NiMiYwLXy8EDwFsici2wHviVd5wM4AZVHVWOxyp/cXGQlgaffQZXXQU1bGxFY0zVVdpvuHKt91fV7araX1XbquoAVc3xtmcGSxKq+oqqji7PGMps4EDYuhU++cTvSIwxpkLZT+Hj1bMntGkDb7zhRpQ1xpgqKtTosXtEZLeI7AZSCpYLtldijOFJBK64AjZtgtnVc0QTY0z1EGpQwHhVTfBu0QHL8aqaUJlBhq2MDNeY/d//wsGDJZc3xpgIFOqMom5JTy5NmSqt4Kxi2zaYPt3vaIwxpkKEaqOYLCKPiEgfEYkr2CgibUTkWhGZDgyq+BDDXGoqdOwIb70FBw74HY0xxpS7UFVP/XED9v0WWCYiu0RkO/A6rsf1SFWdWDlhhrGCs4odO2DaNL+jMcaYchdyCA9VnQbYt19JkpNdv4qpU2HIEBvWwxhTpdjlseWlVy/XryI72+9IjDGmXFmiKC9durj7BQv8jcMYY8qZJYrycsIJ0KIFfPut35EYY0y5Ks2c2Y+ISHJlBBPx0tNhyRLIy/M7EmOMKTelOaNYATwnIl+LyA3eJEYmmC5dXJJYtszvSIwxptyUmChU9QVV7QVcCSQCi0XkDREp9SRG1UbHjlCzprVTGGOqlFK1UYhIFNDeu20DFgG3iciECowt8sTGuktl58/3OxJjjCk3pWmjeBRYBZwH3Keq6ar6oKpeCHSu6AAjTpcusGGDG9bDGGOqgNKcUSwGUlX1t6o6r8hj3SogpshWcJmsXf1kjKkiQvbM9iwCTpUjexvvAtar6q4KiSqStWoFjRq56qeBA/2Oxhhjyqw0ieLfQBfcmYUAHYFlQD0RuVFVZ1RgfJFHBDp3hi+/hPx8iIryOyJjjCmT0lQ9/QR0VtUMVU3HtUusBQYCD1VkcBErPR327YM1a/yOxBhjyqw0iaKdqhZ2DFDV5UB7VV1bcWFFuNRUd2ZhVz8ZY6qA0lQ9LReRp4GCS2GHe9tiAZvWLZj4eDj5ZOt4Z4ypEkpzRjES+A4Y493WAlfhkoR1uivOKafAunWg6nckxhhTJiHPKLyOdtNUtS/wSJAieyskqqogKQk+/NANPX7CCX5HY4wxxy3kGYWq5gOHy3t8JxFpKCIficga775BMeVaicgMEVkhIstFJLE846hQJ5/s7tet8zcOY4wpo9JUPe0FlojIiyLyeMGtjMcdC8xS1ba46VbHFlPuNeBhVT0N17lvSxmPW3lat3YN2t9/73ckxhhTJqVpzH7Hu5WnwcBZ3vKrwBzgzsACItIBiFbVjwBUNbKquWrVghNPhLV2cZgxJrKVmChU9VURqQ20UtVV5XTcpqq60VveBDQNUqYdsFNE3gGSgJnAWK867Agicj1wPUCrVq3KKcRycPLJsHKl31EYY0yZlGZQwAuBhcCH3nqaiEwpxfNmisjSILfBgeVUVYFglwZFA72B24GuQBvc1VZHUdXnvA6BGU2aNCkptMqTlARbtsCePX5HYowxx600VU/jcO0DcwBUdaGItCnpSao6oLjHRGSziDRX1Y0i0pzgbQ/ZwMKCjn0i8i7QA3ixFDGHh8AG7ZQUf2MxxpjjVJrG7INBBv87XMbjTsH1z8C7nxykzDdAfREpOEXoBywv43ErVxsvn1o7hTEmgpUmUSwTkV8DUSLSVkSeAL4o43EfAAaKyBpggLeOiGSIyAtQeGnu7cAsEVmCG5Dw+TIet3LVqwcNG1qiMMZEtNJUPd0M3A0cAP4LTAfuLctBVXU70D/I9kxgVMD6R0Bk19m0aWN9KYwxEa00Vz3txyWKuys+nCqoTRs3h3ZeHsTE+B2NMcYcsxIThYi0w1UBJQaWV9V+FRdWFdKmDRw+DD/84MZ/MsaYCFOaqqf/Ac8ALwBH9WEwJQhs0LZEYYyJQKVJFIdU9ekKj6SqatYMate2Bm1jTMQqzVVPU0XkdyLS3BvMr6GINKzwyKoKEdfxzhKFMSZCleaMoqC/wx0B2xTXU9qUxsknw0cfubkpRPyOxhhjjklprnpKqoxAqrSkJMjNhY0b3UCBxhgTQYqtehKRPwYsDyvy2H0VGVSVUzCUh1U/GWMiUKg2issClu8q8tigCoil6jrpJIiOtpFkjTERKVSikGKWg62bUGrWhPR0mDsX8u0KY2NMZAmVKLSY5WDrpiT9+kFODixa5HckxhhzTEIlilQR2S0ie4AUb7lgvVMlxVd1dO0K8fEwa1bocqqwfDnMmOGWjTHGZ8Ve9aSqUZUZSJVXsyb07g0zZ8K+fRAXd+Tj27bB7Nnu8Y3e5H95eXDBBZUfqzHGBChNhztTXvr3d1/+XxQZpT0zE0aNgv/8Bxo3ht//3rVpvPwy/PijP7EaY4zHEkVlatsWWrY8svopJwcefRRatYIXXoD77nPtGbfcArGx7jFrADfG+MgSRWUScUlg2TLYtMm1QTz6KBw4AH/8IzRt+kvZhg3hxhth1Sp4+23/YjbGVHuWKCpb374uYcyeDe+8AwsXwvXXuzONonr3hj594I03rLOeMcY3ligqW+PGkJoK06a5NolevWDgwOLL33CDm1L1H/9ww4AYY0wls0Thh379YNcuV700enTogQLj413jdnY2PPGEXTJrjKl0pRk91pS300+HJUtg0CCoW7fk8mlp8JvfwGuvuQbxIUMqOkJjjClkicIPsbHuqqZjMXQorFnjLplt0wZSUiomNmOMKcKqniKFCIwZ44Ypf+gh10HPGGMqgSWKSFKnDtx9t+u0N2YMPP44fPmlNXIbYyqUL1VP3lSqbwKJQBbwK1XdEaTcQ8D5uIT2EXCrajVvzW3ZEv76V5gyxfXw/ugjNzxIjx5w6aW/zH1hjDHlxK82irHALFV9QETGeut3BhYQkdOBXkBBZfxnwJnAnEqMMzy1b+9uhw65AQS/+sr19p47Fzp3du0ZnTrZtKvGmHLhV9XTYOBVb/lVYEiQMgrUAmKAWKAmsLkygosY0dGuUfv66+Gll2DkSFi3zlVPTZrkd3TGmCrCr0TRVFW9IVLZBDQtWkBVvwQ+BjZ6t+mquiLYzkTkehHJFJHMrVu3VlTM4S0uzp1JvPiiG1Dwf/+D/fv9jsoYUwVUWKIQkZkisjTIbXBgOa/N4ah2BxE5BTgNaAm0APqJSO9gx1LV51Q1Q1UzmjRpUgF/TQSJiYERI2DvXtf72xhjyqjC2ihUdUBxj4nIZhFprqobRaQ5sCVIsYuBr1R1r/ecD4CewNwKCbgqadvWnVVMmuTms6hVy++IjDERzK+qpynASG95JDA5SJkfgDNFJFpEauIasoNWPZkghg+H3bvhgw/8jsQYE+H8ShQPAANFZA0wwFtHRDJE5AWvzETge2AJsAhYpKpT/Qg2Ip12mmvofucd1+/CGGOOky+JQlW3q2p/VW2rqgNUNcfbnqmqo7zlfFX9raqepqodVPU2P2KNaJddBjt3wvTpfkdijIlg1jO7KuvYETp0cBMfHTzodzTGmAhliaIqE3FnFdu3w2ef+R2NMSZCWaKo6tLS3LwXX33ldyTGmAhliaKqE4Hu3WHBAmvUNsYcF0sU1UH37m6E2UWL/I7EGBOBLFFUBykpULs2fP2135EYYyKQJYrqoGZN11N73jybc9sYc8wsUVQX3bvDjh2werXfkRhjIowliuoiPR1q1LDqJ2PMMbNEUV3Ex7sOeHaZrDHmGFmiqE66d4cNG2DjxpLLGmOMxxJFddKjh7u3swpjzDGwRFGdnHACJCZaO4Ux5phYoqhuevSA5cth1y6/IzHGRAhLFNVNr16uL8Xs2X5HYoyJEJYoqpvERDf0+LRp1vnOGFMqliiqowsugE2bYP58vyMxxkQASxTVUc+e0KABvP++35EYYyKAJYrqKDoaBg1yZxTWp8IYUwJLFNXVoEFuSI8PPvA7EmNMmLNEUV01bOiqoGbMgAMH/I7GGBPGLFFUZ+efD/v2waef+h2JMSaMWaKozpKT3eWyU6fC4cN+R2OMCVO+JAoRGSYiy0TksIhkhCg3SERWich3IjK2MmOsFkTgkktg3Tq4/36bU9sYE5RfZxRLgUuAYus8RCQKeAo4F+gAXC4iHSonvGqkb1+4/no3/tPdd8Pu3X5HZIwJM74kClVdoaqrSijWDfhOVdeqah4wARhc8dFVQxdeCGPHwvffwx//6DrjGWOMJ5zbKFoAGwLWs71tRxGR60UkU0Qyt27dWinBVTmnnw5/+5sbLPD222HNGr8jMsaEiQpLFCIyU0SWBrmV+1mBqj6nqhmqmtGkSZPy3n310aEDPPwwxMa6MwwbjtwYQwUmClUdoKodg9wml3IXPwInBay39LaZitSyJTzyiLsa6u9/d1dEGWOqtXCuevoGaCsiSSISA1wGTPE5puqhfn247z7o1g2eew5eecVGmjWmGvPr8tiLRSQb6Am8LyLTve0nisg0AFU9BIwGpgMrgLdUdZkf8VZLsbHwf/8H550Hb78Nb73ld0TGGJ9E+3FQVZ0ETAqy/SfgvID1acC0SgzNBKpRA264AXJz4fXXIT7eJQ5jTLXiS6IwEUQEbrnFDfXxzDMQFwdnnul3VMaYShTObRQmXERFuf4VHTvCo4/Chx/CoUN+R2WMqSSWKEzpxMTAPfdAu3bw1FOuN/eUKa5ayhhTpYlWsatZMjIyNDMz0+8wqi5VyMx0DdzLlkHdutC0KeTnu1uNGpCWBv36QVKSq7oyxoQ9EZmvqkHH3rM2CnNsRKBrV3dbuRKmTYO9e12CiIpyZxjvvw+TJ0Pr1jBggBsiJCrK78iNMcfJEoU5fu3bu1tRe/bAZ5/B7Nnw4ouQlQW33mpnF8ZEKEsUpvzFx8O557rbhAkwfry7WmrUKEsWxkQgSxSmYg0f7qqmJk927RmXX+53RMaYY2SJwlQsEbj2WtcP4403oE4dGGyjxRsTSSxRmIonAqNHw/798MILsGMHXHmlawAPZds2N0Vr48YllzXGVBhLFKZyREXBHXfA88+7S2vXr3fzXsTFHVkuLw+++gqmT4fFi9226Gg44QRo3hx69YKzzoKaNSv9TzCmurJ+FKbyTZvmRqVt3hyuu85Nv7pxI/z4IyxY4K6aatoUBg6EBg3cY5s2uaunsrOhYUO46CIYNOjoRFOSnBw3KdPq1W6fO3e6M5ydO10fkdhYd6tVy8Vw0knuMt9WrdwQ7NH228pUTaH6UViiMP5YsgTuv98lBXDVU40bu8ttzz4bUlOPvkJKFRYtcmckCxe6L/OMDOjSBdLTXQIJlJsL330Hq1a52+rVsH27e6xGDZcIGjRwt/r13bbcXDhwAH7+GX76ySWpw4fdc2JioE0baNsWTjkFTjwRmjWDevXsai4T8SxRmPCUk+Pm6W7WzH1px8SU/rlr17qOfZmZbj/gzlDAfdHn5bkG9ILPd/PmbviRglubNqU73sGD7kwnK8slndWrXcx5eb+UiY2FJk1cwkhIcPcNGrhjNmvm7i2ZmDBnicJUXaquvSMz032RR0W5L+6YGNefoyAx1KtXfsfMz/+lOqzgtm2bm2989+5f7gP/t+LiXPXVSSe5+wYN3BlR7dou1v37XfVXwfNzc91ZTW6uG4AxLu7IW5067rl16rgzoQMHfrnl5h55A/e6FNxq1jzyVpDAgt3XqOGq2wpuBfuoUcPdij6n4LGCskWfKxI6YRb08C/Yf8ExAp+jGnwirZL2bUKyITyOwV13uVEn+vd3/59/+pOrCenb1/0PjhvnpmTo3dv9YP3b39wIFaef7v7H778fLr7YTQ63Ywc89BAMHepqRrZtc7OMDh/uhkPatAkeewxGjHADs/74Izz5pLsg6LTT3PffM8/ANde42o61a11b8HXXuR/Ea9bASy+5KSNat4YVK+C119wFRi1awNKlrq/brbe6H7YLF8Kbb8If/uBqeebPh4kT3cCwDRrAvHkwaZJ7DRIS4Isv3Eyo99zjvpvmznXNC+PGue/ijz+GGTPg3nvd98CsWTBzpnsNwLVHz53rXiNwz503zz0f3JiCixa51xjcsVeudMcHF9vatS4+cH33fvzRxQ/ub9u6VRgzJhESE3n1VVeTNXq0e/yll+BAJtzY1a0//7y7v+46d//00+7vuOYat/7kky63jBzp1v/1L3eiMGKEW3/kEfe6XnZZFLRsyUNvtKRNGxg6yj1+//3Qvpt7/8nL497/+5nU5lu4qN1KyM5m3Pi2dItewHnxLwBwz8or6N1wGeec8K377K24kgFNFtL/hKUcio3jTytHcHbiavq22sCBPXmMm3cu5zX4it6NlrPvUCx/WzOcC5vO4/SGK9l9sDb3fzeMi5t9SbcGa9iRF8dDa4cyNHEp6Q2z2La/Do+suoDhzT8lLWEtm3Lr89i6ixjR4kM6JvzAjz835MmsC7iy5WxOi89m/f4mPLP+XK456SPa1t3I2n1Nef6Hc7iu1XTaxG1mzd7mvLRhIDe0/oDWdbayYk9LXsvux+jE92hRO4elu1sx/sezuDVpCs1q7WThriTe/Kk3fzj5XRrH7Gb+zpOZuLEXfzz5bRrE7GPejrZM2tSTu075Hwk1f+aLnPZM3dyNe9q+SVz0AebmJDNtSwbj2o4nNuoQH2/rxIytnbn31NeJrnGYWdtSmLk1jfs7joeoKKZv6czcnA78LXUi1KjBtJ/SmLetDeM6vQ2qTMnuwqKdrflT8jvus5fdlZV7WnBXp/egRg0m/tCNtXua8MfU6SDChO+78uP+Bvyh0wxQZfz33dmaG8+Y5JkAvLrmdPYcrMXoDrPdZ2/1GRzIj+bG0+a4z96qPu6zd+qn7rO34ixiow5xTbvPQIQnl/clvuYBRrb9wn32lg6gSa09jDjFzVv/yJKzaRG3k8tO/gZEeGjRObSJ38rQpPnus7fwXNrX38TF/Xb98g9UjixRGFMRYmKgXgy0rQcXtnXbNgNd+8MZv3UZ7b46kNIN0i90mfjfJ8GgC+D8WpAv8Cfg7POhL3AAGAeccxGk74dtP8PDMXBGV0jdBwdrwwuN4YL+0Ksm7I+FR2vCsKGQDmwDHgF+dQl0yofsQ/C4wCUD4bR8+EnghVowfAC0PQQ/CLxaG4afCa3zYC0wvg4M7QItDsC6aHinPgxNhuZ5sDYW3msIw06FRrmwOgY+bASXtICE/bCqNsxpCkMSICEPvqsHnzeHS+Oh7kFYXR++bAaX1oHYA7C8PmQ2hyFAzEFY0QgWtYRLfgWxwPKmsLgZ/OrXEKWwtCksPgEuvdT9wlvUDFY3hT59XBvT0lawvimccYZ7LxYnwo+N3RV0IrCgNWysBz16uPJRSbA13v2CU4UDLWFXXejQwT3/wEmwt7Z7HODnlpBb85fy+1rAoahfHt9zorsvWN91IkTnQ3KyW89pDrF5v+x/W1OIS/hlfUtTqBcH7b02vZ+aQMOa0H6fi/+nJnBCNLSuXw4f3qNZ1ZMxxpiQVU/Wi8kYY0xIliiMMcaEZInCGGNMSJYojDHGhGSJwhhjTEiWKIwxxoRkicIYY0xIliiMMcaEVOU63InIVmD9MTylMa7fargJ17ggfGML17ggfGML17jAYjseZYmrtao2CfZAlUsUx0pEMovrjeincI0Lwje2cI0Lwje2cI0LLLbjUVFxWdWTMcaYkCxRGGOMCckSBTzndwDFCNe4IHxjC9e4IHxjC9e4wGI7HhUSV7VvozDGGBOanVEYY4wJyRKFMcaYkKptohCRQSKySkS+E5GxPsfykohsEZGlAdsaishHIrLGu2/gQ1wnicjHIrJcRJaJyK1hFFstEZknIou82P7ibU8Ska+99/VNEYmp7Ni8OKJE5FsReS/M4soSkSUislBEMr1t4fB+1heRiSKyUkRWiEjPMInrVO+1KrjtFpExYRLb773P/lIR+a/3P1Ehn7NqmShEJAp4CjgX6ABcLiIdfAzpFWBQkW1jgVmq2haY5a1XtkPAH1S1A9ADuMl7ncIhtgNAP1VNBdKAQSLSA3gQeFRVTwF2ANf6EBvArcCKgPVwiQugr6qmBVxvHw7v52PAh6raHkjFvXa+x6Wqq7zXKg03qex+YJLfsYlIC+AWIENVOwJRwGVU1OdMVavdDegJTA9Yvwu4y+eYEoGlAeurgObecnNgVRi8bpOBgeEWG1AHWAB0x/VKjQ72PldiPC1xXx79gPcACYe4vGNnAY2LbPP1/QTqAevwLq4Jl7iCxHk28Hk4xAa0ADYADYFo73N2TkV9zqrlGQW/vMgFsr1t4aSpqm70ljcBTf0MRkQSgc7A14RJbF71zkJgC/AR8D2wU1UPeUX8el//BfwROOytNwqTuAAUmCEi80Xkem+b3+9nErAVeNmrrntBROLCIK6iLgP+6y37Gpuq/gj8A/gB2AjsAuZTQZ+z6pooIoq6nwe+XccsInWBt4Exqro78DE/Y1PVfHVVAi2BbkB7P+IIJCIXAFtUdb7fsRTjDFXtgqt2vUlE+gQ+6NP7GQ10AZ5W1c7APopU5YTB/0AMcBHwv6KP+RGb1yYyGJdkTwTiOLr6utxU10TxI3BSwHpLb1s42SwizQG8+y1+BCEiNXFJYryqvhNOsRVQ1Z3Ax7hT7foiEu095Mf72gu4SESygAm46qfHwiAuoPCXKKq6BVfX3g3/389sIFtVv/bWJ+ISh99xBToXWKCqm711v2MbAKxT1a2qehB4B/fZq5DPWXVNFN8Abb0rBGJwp5RTfI6pqCnASG95JK59oFKJiAAvAitU9Z9hFlsTEanvLdfGtZ2swCWMoX7Fpqp3qWpLVU3Efa5mq+oIv+MCEJE4EYkvWMbVuS/F5/dTVTcBG0TkVG9Tf2C533EVcTm/VDuB/7H9APQQkTre/2nBa1YxnzM/G4f8vAHnAatx9dp3+xzLf3H1jAdxv66uxdVrzwLWADOBhj7EdQbulHoxsNC7nRcmsaUA33qxLQX+n7e9DTAP+A5XTRDr4/t6FvBeuMTlxbDIuy0r+NyHyfuZBmR67+e7QINwiMuLLQ7YDtQL2OZ7bMBfgJXe5/8/QGxFfc5sCA9jjDEhVdeqJ2OMMaVkicIYY0xIliiMMcaEZInCGGNMSJYojDHGhGSJwkQUEVEReSRg/XYRGVdO+35FRIaWXLLMxxnmjZD6cZHtJ4rIRG85TUTOK8dj1heR3wU7ljElsURhIs0B4BIRaex3IIECesOWxrXAdaraN3Cjqv6kqgWJKg3XZ6W8YqgPFCaKIscyJiRLFCbSHMLNC/z7og8UPSMQkb3e/Vki8omITBaRtSLygIiMEDefxRIROTlgNwNEJFNEVnvjNhUMPviwiHwjIotF5LcB+50rIlNwvWKLxnO5t/+lIvKgt+3/4ToyvigiDxcpn+iVjQH+Cgz35kAY7vWqfsmL+VsRGew95yoRmSIis4FZIlJXRGaJyALv2IO93T8AnOzt7+GCY3n7qCUiL3vlvxWRvgH7fkdEPhQ378JDx/xumSrhWH4FGRMungIWH+MXVypwGpADrAVeUNVu4iZjuhkY45VLxI1/dDLwsYicAlwJ7FLVriISC3wuIjO88l2Ajqq6LvBgInIibm6AdNy8ADNEZIiq/lVE+gG3q2pmsEBVNc9LKBmqOtrb33244UCu8YYumSciMwNiSFHVHO+s4mJV3e2ddX3lJbKxXpxp3v4SAw55kzusdhKR9l6s7bzH0nCjBh8AVonIE6oaOPKyqQbsjMJEHHUj2L6Gm7iltL5R1Y2qegA3bEvBF/0SXHIo8JaqHlbVNbiE0h43JtKV4oY0/xo3fENbr/y8oknC0xWYo27QtkPAeKBPkHKldTYw1othDlALaOU99pGq5njLAtwnIotxQ0u0oOQhsM8AXgdQ1ZXAeqAgUcxS1V2qmos7a2pdhr/BRCg7ozCR6l+4yYpeDth2CO/Hj4jUAAKngTwQsHw4YP0wR/4fFB3TRnFfvjer6vTAB0TkLNyQ2JVBgEtVdVWRGLoXiWEE0ARIV9WD4kaxrVWG4wa+bvnYd0a1ZGcUJiJ5v6Df4sipHrNwVT3g5g6oeRy7HiYiNbx2iza4mcymAzeKG3IdEWnnjb4ayjzgTBFpLG7q3cuBT44hjj1AfMD6dOBmb6RQRKRzMc+rh5sP46DX1lBwBlB0f4Hm4hIMXpVTK9zfbQxgicJEtkeAwKufnsd9OS/CzU1xPL/2f8B9yX8A3OBVubyAq3ZZ4DUAP0sJv6zVzX42Fjfs8yJgvqoey5DPHwMdChqzgXtxiW+xiCzz1oMZD2SIyBJc28pKL57tuLaVpUUb0YF/AzW857wJXOVV0RkDYKPHGmOMCc3OKIwxxoRkicIYY0xIliiMMcaEZInCGGNMSJYojDHGhGSJwhhjTEiWKIwxxoT0/wFZaZQgViVXjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "result = numpy.load('./output/summary_data.npz')\n",
    "\n",
    "eig_val, eig_state = numpy.linalg.eig(\n",
    "                     molecular_hamiltonian.construct_h_matrix())\n",
    "min_eig_H = numpy.min(eig_val.real)\n",
    "min_loss = numpy.ones([len(result['iter'])]) * min_eig_H\n",
    "\n",
    "plt.figure(1)\n",
    "func1, = plt.plot(result['iter'], result['energy'], \n",
    "                  alpha=0.7, marker='', linestyle=\"-\", color='r')\n",
    "func_min, = plt.plot(result['iter'], min_loss, \n",
    "                  alpha=0.7, marker='', linestyle=\":\", color='b')\n",
    "plt.xlabel('Number of iteration')\n",
    "plt.ylabel('Energy (Ha)')\n",
    "\n",
    "plt.legend(handles=[\n",
    "    func1,\n",
    "    func_min\n",
    "],\n",
    "    labels=[\n",
    "        r'$\\left\\langle {\\psi \\left( {\\theta } \\right)} '\n",
    "        r'\\right|H\\left| {\\psi \\left( {\\theta } \\right)} \\right\\rangle $',\n",
    "        'Ground-state energy',\n",
    "    ], loc='best')\n",
    "\n",
    "#plt.savefig(\"vqe.png\", bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determining the interatomic distance\n",
    "\n",
    "Recall the above calculation is done with an interatomic distance $d = 74$ pm between two hydrogen atoms. Another interesting aspect we can try with VQE is determining the true interatomic distance by modifying the `h2.xyz` file. The results are summarize in figure below,\n",
    "\n",
    "![vqe-fig-dist](figures/vqe-fig-distance.png)\n",
    "\n",
    "The lowest value is found around $d = 74$ pm (1 pm = $1\\times 10^{-12}$m), which is consistent with the [experimental data](https://cccbdb.nist.gov/exp2x.asp?casno=1333740&charge=0) $d_{exp} (H_2) = 74.14$ pm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## References\n",
    "\n",
    "[1] Cao, Yudong, et al. Quantum Chemistry in the Age of Quantum Computing. [Chemical reviews 119.19 (2019): 10856-10915.](https://pubs.acs.org/doi/10.1021/acs.chemrev.8b00803)\n",
    "\n",
    "[2] McArdle, Sam, et al. Quantum computational chemistry. [Reviews of Modern Physics 92.1 (2020): 015003.](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.92.015003)\n",
    "\n",
    "\n",
    "[3] Peruzzo, A. et al. A variational eigenvalue solver on a photonic quantum processor. [Nat. Commun. 5, 4213 (2014).](https://www.nature.com/articles/ncomms5213)\n",
    "\n",
    "[4] Moll, Nikolaj, et al. Quantum optimization using variational algorithms on near-term quantum devices. [Quantum Science and Technology 3.3 (2018): 030503.](https://iopscience.iop.org/article/10.1088/2058-9565/aab822)\n",
    "\n",
    "[5] Helgaker, Trygve, Poul Jorgensen, and Jeppe Olsen. Molecular electronic-structure theory. John Wiley & Sons, 2014.\n",
    "\n",
    "[6] Dirac, Paul Adrien Maurice. Quantum mechanics of many-electron systems. [Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 123.792 (1929): 714-733.](https://royalsocietypublishing.org/doi/10.1098/rspa.1929.0094)\n",
    "\n",
    "[7] Szabo, Attila, and Neil S. Ostlund. Modern quantum chemistry: introduction to advanced electronic structure theory. Courier Corporation, 2012."
   ]
  }
 ],
 "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "426.667px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}