{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Eigensolver (VQE)\n",
    "In this section, the variational quantum eigensolver (VQE) is run on the simulator using Qulacs to find the ground state of the quantum chemical Hamiltonian obtained using OpenFermion/PySCF.\n",
    "\n",
    "This notebook is translated from https://dojo.qulacs.org/ja/latest/notebooks/6.2_qulacs_VQE.html\n",
    "\n",
    "Necessary packages:\n",
    "\n",
    "* qulacs\n",
    "* openfermion\n",
    "* openfermion-pyscf\n",
    "* pyscf\n",
    "* scipy\n",
    "* numpy\n",
    "\n",
    "## Install and import necessary packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please execute if various libraries are not installed\n",
    "## When use Google Colaboratory, please ignore 'You must restart the runtime in order to use newly installed versions'.\n",
    "## Crash when restarting runtime.\n",
    "!pip install qulacs pyscf openfermion openfermionpyscf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qulacs\n",
    "from openfermion.transforms import get_fermion_operator, jordan_wigner\n",
    "from openfermion.linalg import get_sparse_operator\n",
    "from openfermion.chem import MolecularData\n",
    "from openfermionpyscf import run_pyscf\n",
    "from scipy.optimize import minimize\n",
    "from pyscf import fci\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Hamiltonian\n",
    "Create Hamiltonian by PySCF in the same procedure as described in the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis = \"sto-3g\"\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "distance  = 0.977\n",
    "geometry = [[\"H\", [0,0,0]],[\"H\", [0,0,distance]]]\n",
    "description  = \"tmp\"\n",
    "molecule = MolecularData(geometry, basis, multiplicity, charge, description)\n",
    "molecule = run_pyscf(molecule,run_scf=1,run_fci=1)\n",
    "n_qubit = molecule.n_qubits\n",
    "n_electron = molecule.n_electrons\n",
    "fermionic_hamiltonian = get_fermion_operator(molecule.get_molecular_hamiltonian())\n",
    "jw_hamiltonian = jordan_wigner(fermionic_hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert Hamiltonian to qulacs Hamiltonian\n",
    "In Qulacs, observables like Hamiltonians are handled by the `Observable` class. There is a function `create_observable_from_openfermion_text` that converts OpenFermion Hamiltonian to Qulacs `Observable`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qulacs import Observable\n",
    "from qulacs.observable import create_observable_from_openfermion_text\n",
    "qulacs_hamiltonian = create_observable_from_openfermion_text(str(jw_hamiltonian))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct ansatz\n",
    "Construct a quantum circuit on Qulacs. Here, the quantum circuit is modeled after the experiments with superconducting qubits (A. Kandala et. al. , “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets“, Nature **549**, 242–246)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qulacs import QuantumState, QuantumCircuit\n",
    "from qulacs.gate import CZ, RY, RZ, merge\n",
    "\n",
    "depth = n_qubit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def he_ansatz_circuit(n_qubit, depth, theta_list):\n",
    "    \"\"\"he_ansatz_circuit\n",
    "    Returns hardware efficient ansatz circuit.\n",
    "\n",
    "    Args:\n",
    "        n_qubit (:class:`int`):\n",
    "            the number of qubit used (equivalent to the number of fermionic modes)\n",
    "        depth (:class:`int`):\n",
    "            depth of the circuit.\n",
    "        theta_list (:class:`numpy.ndarray`):\n",
    "            rotation angles.\n",
    "    Returns:\n",
    "        :class:`qulacs.QuantumCircuit`\n",
    "    \"\"\"\n",
    "    circuit = QuantumCircuit(n_qubit)\n",
    "    for d in range(depth):\n",
    "        for i in range(n_qubit):\n",
    "            circuit.add_gate(merge(RY(i, theta_list[2*i+2*n_qubit*d]), RZ(i, theta_list[2*i+1+2*n_qubit*d])))\n",
    "        for i in range(n_qubit//2):\n",
    "            circuit.add_gate(CZ(2*i, 2*i+1))\n",
    "        for i in range(n_qubit//2-1):\n",
    "            circuit.add_gate(CZ(2*i+1, 2*i+2))\n",
    "    for i in range(n_qubit):\n",
    "        circuit.add_gate(merge(RY(i, theta_list[2*i+2*n_qubit*depth]), RZ(i, theta_list[2*i+1+2*n_qubit*depth])))\n",
    "\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define VQE cost function\n",
    "As explained in Section 5-1, VQE obtains an approximate ground state by minimizing the expected value \n",
    "\n",
    "\\begin{equation}\n",
    "\\left<H(\\theta)\\right>=\\left<\\psi(\\theta)|H|\\psi(\\theta)\\right>\n",
    "\\end{equation}\n",
    "\n",
    "of Hamiltonian for state $\\left|\\psi(\\theta)\\right>=U(\\theta)\\left|0\\right>$ output from the quantum circuit $U(\\theta)$ with parameters. The following defines a function that returns the expectation value of this Hamiltonian.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(theta_list):\n",
    "    state = QuantumState(n_qubit) #Prepare |00000>\n",
    "    circuit = he_ansatz_circuit(n_qubit, depth, theta_list) #Construct quantum circuit\n",
    "    circuit.update_quantum_state(state) #Operate quantum circuit on state\n",
    "    return qulacs_hamiltonian.get_expectation_value(state) #Calculate expectation value of Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run VQE\n",
    "Now everthing is prepared, run VQE. For optimization, the BFGS method implemented in scipy is applied, and initial parameters are randomly selected. This process should end in tens of seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cost_history = []\n",
    "init_theta_list = np.random.random(2*n_qubit*(depth+1))*1e-1\n",
    "cost_history.append(cost(init_theta_list))\n",
    "method = \"BFGS\"\n",
    "options = {\"disp\": True, \"maxiter\": 50, \"gtol\": 1e-6}\n",
    "opt = minimize(cost, init_theta_list,\n",
    "               method=method,\n",
    "               callback=lambda x: cost_history.append(cost(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results can be plotted to see that they converge to the correct solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEdCAYAAACsS3i2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1dnA8d8TIBnCjoCyL6IoUhowCPIqbri8WhUVpRZQXIoVqbtQFRUV9wVcKyiLLFaKlVdpFZcqxdZaNlERAS37HkB2AoQ87x/nTpjMTEJuMpOZZJ7v53M/k5xz751nLiFPzrnnniOqijHGGJMM0hIdgDHGGBNkSckYY0zSsKRkjDEmaVhSMsYYkzQsKRljjEkalpSMMcYkDUtKxhhjkoYlJWOMMUnDV1ISkeYiMk5E1orIARE52ytv6JV3iU+YxhhjUkGJk5KItAbmAVcA3wNVgnWqmgNkAzfGOkBjjDGpo6qPfR8D8oEOwD5gc1j9B8DFMYrLGGNMCvLTfdcTeFVV1wDRJsxbBTSLSVTGGGNSkp+WUm1gQzH16T7PVyk0aNBAW7VqlegwjDGmQpk/f/4WVW0YXu4niawBTiqmvhvwk9/AKrpWrVoxb968RIdhjDEVioisilbup/vuXeB6EekQUqbeya8ArgT+XOoIjTHGpDw/SekxYC3wH2AyLiH9QUT+jUtG3wDPxTxCY4wxKaPESUlVdwKnAm/ghn8LcC7QDngVOEtVc+MRZDQikiYid4jIEhHJFZE1IvKciNQo4fFaxLY73rEbY4yJztfABC8x3QbcJiINcYkpRxOzfO1I4FZgOq6FdqL3fScR6amq+SU4xxfAmLCygzGN0hhjTImVerSc98BsQojIScDvgXdV9YqQ8hXAi8CvgbdKcKrlqjo5PlEaY4zxq8RJSUR6lGQ/VZ1d+nBK7GpcK21UWPnrwJNAP0qWlBCRdCBdVa3bzhhjEsxPS2kW0R+aDVflyLuUWRfc7BJzQgtVNVdEFnr1JdEbl8CqiEgOMBUYpqo7YhmsMcaYkvGTlK4r4vhjgQHASmB02UMqkSbAFlXdH6VuHdBdRNJV9UAx55gDTMM9W1UbuBAYDJwhIt3j3nKaNAn27IHf/S6ub2PK344dO9iyZQsHDhT342dM5VKlShVq1apF/fr1ycjIKPV5SpyUVPXNoupE5BlgQamj8C8TiJaQAHJD9inyt4Kqdg0rmigi3+KGvt/mvUYlIgOBgQAtWrQoYchhpk6FDRssKVUyubm5bNq0iWbNmlG9enVEJNEhGRN3qsrBgwfZuXMnq1evpkWLFqVOTDFZT0lVf8YNFR8Si/OVwF6gqE8cCNnHr2dwieyi4nZS1TGqmq2q2Q0bRsySUTKBAOwvKq+aiionJ4eGDRuSmZlpCcmkDBEhPT2dBg0aUK9ePbZt21bqc8Vykb+fgTYxPF9x1gMNRCRaYmqK69rz3XeiqgeD5y5jfEeWkQG55fZYlyknubm51KxZM9FhGJMwtWvXZteuXaU+PiZJSUQCQH9gYyzOVwJzcbGfEiWOLNy6T755xzcDNpU1wCMKBCwpVUJ5eXlUrZpy8xIbU6BatWocOnSo1Mf7GRI+roiq+riZHhoC95Q6En+mAvcBt+MegA36Le5e0pRggYgcC1RT1SUhZUep6tYo530Ud01mxCPoQqz7rtKybjuTysr68+/nT7oBRZRvA5YBd6hqiZ4NKitV/U5EXgEGi8i7uAUGgzM6/IPCzyj9HWiJe64paJiIdAM+B1YDNXGj787Cze33Utw/hHXfGWNMBD+j72J5/ykWbscNQx+IG5iwBZdMHizBFEOzgPbAtcBRwCHgR+B+4PlymcPPuu+MMSZChe38VtVDuDnvip2ZXFVbRSl7D3gvPpGVUCAAeXlw6BBUKY/njY0xJvklW+sndQTH8Nt9JVMBXXnllYgICxcuLHIfVaV169bUrVuXffv2FZS/9957XHTRRTRq1Ij09HSaNGlC7969+eKLL6Ke58wzz0REitxGjBgR889nEqfIlpKIfFaK86mqnlOGeFJHwHucKjcXMjMTG4sxPt1www288847jB8/nhdeeCHqPp9//jkrV67kpptuonr16hw6dIjrrruOSZMm0b59e2677TYaN27MqlWrmDhxIj169GDYsGE8+uijEefKyMjgjTfeiPo+WVlZMf1sJrGK675rQ8nmujOlEUxK1lIyFdB5551H8+bNmTJlCs888wzp6ekR+4wfPx5wCQxg+PDhTJo0iWuuuYaxY8cWGjo/dOhQLr30UkaMGEHbtm259tprC52ratWq9OvXL46fyCSLIrvvVLWVqrb2u5Vn8BVasPvOBjuYCigtLY0BAwawdetW3n///Yj6nTt38pe//IUOHTrQpUsXNm/ezLPPPkuLFi0YPXp0xLNcmZmZTJkyhRo1anDffffZvIEpzO4pJUpo950xFdB1112HiBS0iEK9/fbb7Nu3r6CV9Le//Y3c3Fz69+9PIPizH6ZRo0ZceumlrF+/nq+++iqifsuWLVG3vLy82H4wk1AVdvRdhWfdd6nl9tuhmEEBCZGVBaPClyQrudatW3PWWWfx0UcfsWHDBho3blxQN378eNLT0wu63BYtWgRA586diz3nySefzFtvvcW3335Ljx6Hl3Dbs2cPRc0zOXfuXLKzs0v9OUxy8ZWURKQecAPQFahHZEvLBjqUlHXfmUrghhtu4LPPPmPixIkMHToUgCVLlvDVV1/Ru3dvGjRw00ju3LkTgDp16hR7vtq1awNEzJ0WCASYMSP6RCvt2rUr02cwycXPNEMtgX/h1jLagVuDaBuHk9MWYE8cYqycrPsutZShRZLMLr/8curWrcv48eMLktK4cW5Gsuuvv75gv2Cy2bGj+PUzg8nr6KOPLlRepUoVevbsGbO4TfLyc09pBFAXOAc4DjdtTx9ccnoC2AWcHusAKy3rvjOVQCAQ4De/+Q1Lly7lyy+/5NChQ0yaNIlmzZpx/vnnF+zXoUMHABYsKH7ZtWB927Zt4xe0SWp+ktI5wOuq+jmHh4qLqu5V1fuB74CnYh1gpWXdd6aSCA5mGD9+PB9++CEbN27k2muvJS3t8K+Xiy66iEAgwOTJk9lfxB9iOTk5vPfee7Rs2ZLTTjutXGI3ycdPUjoKWOR9fdB7rR5S/wlwbiyCSgnWfWcqic6dO5OVlcXUqVN55ZVXEJFCXXfgRtbdfffdrFq1iptvvjliaYN9+/bRv39/du/ezSOPPFIooZnU4megQw5umQpwXXW5QKuQ+nQKJylTHOu+M5XIDTfcwO9//3tmzpzJmWeeSZs2ket9Dh8+nJUrVzJ+/Hjmzp3L1VdfXWhGh1WrVjF8+HCuueaaiGPz8vKYPHly1Pdu06YN3bt3j/lnMonhJyl9D/wS3BA7EZkDDBKR93EtroHAkmKON6Gs+85UIn379uWee+4hNzc3opUUVKVKFSZNmsQVV1zBmDFjGDlyJNu2bSM/Px8RYebMmZx33nlRj92/fz/9+/cv8r0tKVUefpLSe8BdIlJdVfcBjwAfASu8egUuj3F8lZd135lKpF69eoUmXS1Or1696NWrV8H3U6dOpW/fvjz88MN07949Yjn5WbNmxTJUk+RK3HGrqq+q6rFeQkJVP8OtOPsC8DzQQ1Uj5xsx0Vn3nTEA9OnTh3HjxvHVV19xySWXlDi5mcqpTDM6qOo8YF6MYkkt1n1nTIFrrrkm6r0kk3pK3FISkedFpGM8g0kpVau6xf0sKRljTAE/4y5vB74WkYUicruINIpXUCkjELDuO2OMCeEnKZ0APImb1eF5YK2IzBCR3iISuZiKObKMDGspGWNMCD8DHZap6v2q2go3u8MUoAfwZ2CjiLwqIqfGJ8xKKhCwpGSMMSFK9di0qn6uqtcBxwD9gTnAb4EvYhhb5Wfdd8YYU0iZ5vLwhoevAzbgZniQWASVMqz7zhhjCinVkHAROR64BugLtAAOAR8Cb8YutBRg3XfGGFOIn/WU6gFX45JRF1yraCEwCnhLVXPiEmFlZt13xhhTiJ+W0kZv/0240Xdvquqi4g8xxbLuO2OMKcRPUpqO6577SFXz4xRPagkEYOvWREdhjDFJw8+Q8F+r6oeWkGLIuu+MSSmzZs1CRJgwYUJczn/mmWfSqlWruJy7vNhKWolk3Xemggr+ci1qq1q1TNNqxtSsWbMYPnw427dv93XcjBkzOPfcc2nWrBkZGRk0btyY7t27M2TIELZs2RKnaI9s1KhRcUtqySB5fnJSkY2+MxXc1VdfzYUXXhhRnkwrx86aNYuHH36YAQMGULdu3RIdM3ToUJ5++mk6duzIoEGDOProo1m/fj3fffcdr732GldddRUNGjSIc+TRjRo1ilatWjFgwICIuo8//hhVLf+gYsiSUiJZ952p4Dp37ky/fv0SHUZMbd68mWeffZYuXbrwr3/9i2rVqhWq3717d4IiO7L09Io/41vy/DmTiqz7zqSAIUOGICJMmjSpUPm3335L9erVOeuss8jPd7eq169fz1133UVWVhb16tUjEAjQvn17nnrqKQ4dOhRx7gMHDvD000+TlZVFZmYmderUITs7m5dffhmAAQMG8PDDDwPQunXrgu7F4cOHFxnv8uXLyc/Pp0ePHhEJCaBmzZoRCxFu2bKFW265hebNm5Oenk7z5s255ZZb2FqCgUwTJkxARKIuZhh+j0hEWLVqFf/4xz8KdZeuXLky6v5Bs2fP5txzz6VOnTpUr16dzp07M3bs2CLfb/369Vx99dXUq1ePzMxMzj//fJYtW3bEzxIL1lJKJOu+MxXc3r17o95fSU9Pp3bt2gA89thjzJ49m0GDBtGtWzeOO+449u7dS58+fahRowaTJ08u6O779ttveffdd7nssss49thjOXjwIDNnzuQPf/gDy5cvZ/To0QXvceDAAc4//3xmzZrFeeedR79+/QgEAnz33Xe8++67DB48mJtuuomdO3cyffp0Ro4cWdDl1rFj0avwtGnTBoC//vWv3HnnnTRp0qTYa7Bjxw66d+/OTz/9xPXXX0/nzp35+uuv+eMf/8hnn33GnDlzqFWrlr8LW4RJkyZxxx130KBBA+6///6C8oYNGxZ5zIwZM7jssss45phjuOuuu6hVqxZvv/02N954I8uXL+exxx4rtP+ePXvo0aMH3bp14/HHH2fFihW88MILXHrppSxatIgqVarE5LMUSVVtK8N28skna6k9+KAqqObnl/4cJqksXry4yLozzjgjYnvllVdUVXXPnj1R68ePH6+qqjk5OVHr3377bVVVXb16ddT6999/X1VVlyxZElFXFp9//rkCRW4XXXRRof2XL1+uderU0c6dO+v+/fv1+uuvV6AgvqC9e/dqfpT/D/369dO0tDRdv359QdlTTz2lgN57770R+x86dKjg64ceekgBXbFiRYk/3+DBgxXQ9PR0Pf300/Wee+7RadOm6bZt2yL2ve+++xQo+LcMevnllxXQYcOGFZQFr1vw31VVdfz48Qro559/HnHuM844Q1u2bFmorGXLlkX++4Xvn5eXpy1atNA6derounXrCsr379+v3bt317S0NF22bFmh4wF96qmnCp336aefVkBnzpwZ9X3DFff/IAiYp1F+p/ruvhOR40TkQhHpLyLXhG9lT5MpJLj67IEDiY3DmFIaOHAgn3zyScQW/td369atGTNmDAsWLODss89m3Lhx3HrrrVx88cWF9qtevToibgrNAwcOsG3bNrZs2cL5559Pfn4+8+YdXuh6ypQp1KtXjwcffDAirrIOtHjxxReZOHEi3bt3Z86cOTzzzDNceeWVNG7cmKFDhxbqSpw+fToNGzZk4MCBhc5x00030bBhQ6ZPn16mWMpi/vz5rF69muuvv75Qiy89PZ0hQ4aQn5/Pe++9V+iYtLQ0br311kJlZ599NgA//vhj3GP2M83Q0biHZ88NFkXZTYGJMYgrNQQC7jU393CCMpVWtHsGQZmZmcXWN2jQoNj65s2bF1vfrl27YutL67jjjqNnz54l2veqq67i/fffZ8qUKXTo0IGnn346Yp+8vDyefPJJJk6cyE8//RQxkuznn38u+PrHH38kKyuLQPD/UQyJCP3796d///4cOHCAb7/9lo8//phRo0bx9NNPU7duXe69914AVqxYQXZ2dsQw+KpVq3L88cezYMGCmMdXUitWrADgpJNOiqgLli1fvrxQeZMmTSKu6VFHHQVQontkZeXnntLLuIT0R+AzwKYiKKvgP7yNwDMpYPv27fzzn/8E3ICGzZs307x580L73Hnnnbz00kv06dOH+++/n0aNGlGtWjUWLFjA0KFDCwZElKf09HSys7PJzs7miiuu4MQTT2Ts2LEFSamsgi3DaPLy8mLyHn4Ud88o/I+EePCTlM4FXlPVwfEKJuUEW0c22MGkgBtuuIG1a9fy0ksvcc8999CvXz8+++yzQr8EJ02aRI8ePXj77bcLHfvTTz9FnO/4449nyZIl7N+/n4xiehqK+6XvV7t27ahXrx7r1q0rKGvTpg1Lly4lLy+vUGspLy+PZcuWFQycKEr9+vUB2LZtW0TdihUrIkYA+vk8wff+/vvvI+oWL15caJ9k4afjNQ34Jl6BpKTQ7jtjKrHXXnuNd999l2HDhjF48GCeffZZZs+ezYgRIwrtV6VKlYi/xvfs2cPIkSMjztm3b19+/vnniHNA4b/og8O3o/3Sj2bjxo0sXLgwat0XX3zBtm3baN++fUFZr169yMnJ4Y033ii07+uvv05OTg6XXXZZse93/PHHA/Dpp58WKv/Tn/7E+vXrI/avWbNmiT9L586dadGiBePHj2fjxo0F5QcPHuSZZ55BRLj00ktLdK7y4qel9AXwy3gFkpKs+85UcAsWLGDy5MlR63r16kXNmjVZtGgRd955Jz169OCBBx4A4JZbbuGTTz7h0Ucf5ZxzzuG0004DoHfv3owePZo+ffrQs2dPNm3axLhx4wruaYS67bbbmDFjBiNGjGDu3Lmcd955BAIBvv/+e5YuXVrwS75bt26Am6Whb9++BAIBOnToQIcOHaLGvXbtWrp06ULXrl0555xzaNOmDfv37+ebb75hypQpVKtWjccff7xg/yFDhjBt2jRuueUWFixYQKdOnfj6668ZO3Ys7dq1Y8iQIcVew3bt2tGzZ09Gjx6NqpKVlcXChQuZPn06bdu25eDBg4X279atG2PHjuWBBx7gxBNPJC0tjYsvvpgaNWpEnLtKlSq8/PLLXHbZZXTp0oWBAwdSq1Ytpk6dyldffcV9993HcccdV2x85S7akLxoG9AOWA9cUdJjUmEr05DwGTPckPA5c0p/DpNUSjIUtjI40pBwQH/88Ufdu3evnnTSSVq/fn1ds2ZNoXNs3bpVmzVrpi1atCgYar1nzx69++67tUWLFpqRkaFt27bVJ554Qj/99NOIodSqqvv27dMRI0Zo+/btNSMjQ+vUqaPZ2dkRw7Ofeuopbd26tVatWlUBfeihh4r8bLt27dJXXnlFe/XqpW3atNEaNWpoenq6tmzZUvv27asLFiyIOGbz5s168803a9OmTbVq1aratGlTHTRokObk5ES9buGfY8OGDdq7d2+tVauW1qhRQy+44AJdvHhx1CHhmzZt0ssvv1zr1aunIlJouHu0/VVVZ82apT179tRatWppRkaGZmVl6RtvvBGxX1HHr1ix4ojXLVRZhoSLlvDGlYh8BjQDjvWS03LcirNhOU7PKUOOrHCys7M1dJiqL59+CueeC7Nnw+mnxzYwkxA//PADJ554YqLDMCahSvL/QETmq2p2eLmf7rs2uL+AVnvft/BxrInGuu+MMaaQEiclVW0VxzhSk42+M8aYQmxC1kSy0XfGGFOI7wlZRaQ20BPXnQfu3tInqrorloGlBOu+M8aYQnwlJRG5EXgOqMnhaYYU2C0id6pq5FzopmjWfWeMMYX4mfvuEmAMrmX0ABB8RPgk4PfAGBHZrKozYh5lZWXdd8YYU4ifltIQ4Aegq6qGLr34dxEZD3wFDAUsKZWUdd9VSqoa06ltjKlISvqYUVH8DHT4JTAhLCEFg9iFm0G8XGd8EJE0EblDRJaISK6IrBGR50Qk8tHmOBxfZtZ9V+lUrVo1IZNoGpMsDh48WKaFAP0kpSP96Rf/6WMjjQSeBxbjuhCnAbcCM0SkJJ+trMeXTXq6e7WkVGkEAgF27474u82YlLFz584yrbTrp/vuG2CAiLyqqntCK0SkJjCAcpywVUSC97LeVdUrQspXAC8CvwbeitfxMSHiuvCs+67SaNiwIatXryYjI6PQgnXGVGaqysGDB9m5cyc///wzLVqUfm4FP0npGeBdYIGIvIhrXcDhgQ5tgctLHYl/V+Nab6PCyl8HngT6UXxSKevxsZGRYS2lSiQQCHD00UezceNG9tsfGyaFVKlShVq1atGiRYtilxI5Ej8zOvyfiAwGngJe4nB3nQB7gMGq+l5Rx8dBFyAfmBMWZ66ILPTq43l8bAQClpQqmTp16lCnTp1Eh2FMheTrOSVVfVVE3sIt+NfaKw4+PLsj1sEdQRNgi6pG+3N0HdBdRNJV9UCsjxeRgcBAoEzNVMC674wxJoTvGR1UdTtuQECiZQJF/TbPDdmnqKRU6uNVdQzumS2ys7PLNsDDWkrGGFOgIs99txcoquMyELJPvI6PDbunZIwxBYpsKXnrJylwvqrmed8fSXmup7QeaC8iGVG64JriuuaKaiXF4vjYsO47Y4wpUFz3XRvcQAAJ+T4RzyIVZS5wHnAKbql2AEQkAGQBs+N8fGxY950xxhQosvtOVVupahtVPRjyfesjbeUXOlNxSfL2sPLf4u4FTQkWiMixInJCaY+PK+u+M8aYAr4HOiQLVf1ORF4BBovIu8AHwIm4GRn+QeFnjP4OtCRkVgqfx8dPIABbtpTLWxljTLIr8UAHETkkIr8ppr6PiByKTVgldjtwN+4B3ldwszC8BPxKVfPL4fiys+47Y4wp4KeldKT5Usp9PhVVPYRb3+m5I+zXqizHx5V13xljTIFYDglvAdjqs37Z6DtjjClQbEtJRC4FLg0pGigiPaPsWh+3RPo/YxhbarDuO2OMKXCk7rss3Ozf4Eaq9fC2cLuBL4HBMYssVVj3nTHGFCi2+05VH1bVNFVNw90z6hf8PmyrrarnqepP5RN2JWLdd8YYU8DPQIfWQE68AklZgQAcOgR5eVC1wo7QN8aYmCjxQAdVXaWq8Z8LLtXYkujGGFPA15/mIlIPuAHoCtQjMqmV59x3lUPAm/t1/36oWTOxsRhjTIKVOCmJSEvgX7h1iHYAtYFtHE5OW3CL/Rk/gknJWkrGGOPrOaURQF3gHOA43MCHPrjk9ATuGaXTYx1gpWfdd8YYU8BPUjoHeF1VPydkKXRV3auq9wPf4ZZKN36Edt8ZY0yK85OUjgIWeV8f9F6rh9R/glsm3fhh3XfGGFPAT1LKwc3cAK6rLhdoFVKfTuEkZUrCuu+MMaaAn6T0PfBLcEPsgDnAIBFpISKtgIHAklgHWOlZ950xxhTwMyT8PeAuEamuqvuAR4CPgBVevQKXxzi+ys+674wxpkCJk5Kqvgq8GvL9ZyJyKvAb4BAwXVW/jH2IlZx13xljTIEyzWujqvOAeTGKJTVZ950xxhTws/LschG5pJj6X4nI8tiElUKs+84YYwr4GejQCihuHpwaQMsyRZOKrPvOGGMKxHLl2aMBm7DVL+u+M8aYAkdaebYHcGZI0eUi0jbKrvWBXwMLYxdairDuO2OMKXCkgQ5nAQ95XweHfBc17Psn4I4YxZU6rPvOGGMKHCkpjQIm4CZfXQ7cjnteKZQCu1V1W8yjSwVpaVCtmnXfGWMMR0hKqroDt0wFInIWsFhVbfXZWAsErKVkjDH4G+jwHdC4qEoR6egtAmj8ysiwpGSMMfhLSk/juvKKMh63rpLxKxCw7jtjjMFfUjoLmFFM/ftAz7KFk6Ks+84YYwB/SakJsLqY+rXePsYv674zxhjAX1LaQ/EzNrQErA+qNKz7zhhjAH9J6T/AtSJSK7zCK7sGt8aS8cu674wxBvCXlJ4FmgFfikhvEWnrbb2BL726Z+IRZKVn3XfGGAP4W0/pcxEZBLwATA2rPggMVtVPYxlcyggEYOfOREdhjDEJ52s9JVUdLSJ/Ba4CgnPgLQPeUdV1sQ4uZVj3nTHGAKVY5M9LPiPjEEvqsu47Y4wBSpGURKQGcCpuqYpPVXVTzKNKNTb6zhhjAJ/rKYnIzcA64GNgInCSV95IRHJF5LexDzEFWPedMcYA/pZDvwJ4BfgcuBE3czgAqroZmAn0inWAKcG674wxBvDXUroH+FxVLyNy+QqAeUCHmESVaqz7zhhjAH9J6RfA9GLqNwCNyhZOigoE4MAByM9PdCTGGJNQfpLSoSPs3wQ3FZHxK7j6rLWWjDEpzk9S+gY4P1qFiKQBVwJzYxFUygkE3KslJWNMivOTlF4G/ldEHgXqB48XkXbANNxIvBdjHF9qCCYlG+xgjElxfqYZmioivwDuB+71imfiRuEJMFxVP4x9iCkg2H1nSckYk+L8TjM0TETeBfoCJ+CS0Y/AJFWdF4f4UoN13xljDFC6aYYWAAviEEvqsu47Y4wBSpGUAEQkk8ML/q1S1b2xCykFWfedMcYA/qcZai8iHwDbgUXetl1EPhCRk+IRYEqw7jtjjAF8tJREpBMwC6gJfAIs9qpOAs4D/kdEzlDVhbEOstKz7jtjjAH8dd89A+QDXbz7SgVEpDPwmbfPubELL0VY950xxgD+uu+6AS+HJyQoGPzwCm5Ji3IjIteIyNcisk9ENonIGyLS0MfxK0VEi9gaxDP2Qqz7zhhjAH8tpVxgYzH164F9ZQun5ETkDuB54B/AbUAz4E7gVBE5RVVLOuXREuCxKOW7YhJoSVj3nTHGAP6S0gfAJbgWUTSXAOXy8KzXihmBm9boHFU95JXPBd7HJanHS3i6Tao6OS6BlpR13xljDOCv++5O4CgRmSYiXUSklredIiLv4KYeuiM+YUboBWQCLwUTEoCqzgCWA/38nExEqopI7diG6GnEOwIAAB3cSURBVEPduiACG4triBpjTOXnp6W0GVCgM3B5WF1wwb/NIhJarqpaqmehjqCL9/rvKHVfAVeLSE1V3V2Cc3UF9gLVRGQHbq2oe1V1fWxCLYGaNeGEE2DOnHJ7S2OMSUZ+EsZEXFJKBk2813VR6tbhkmQTYNkRzvM98AbwA1ANOBO3qu453n2pqIlJRAYCAwFatGjhN/boTjkFPvgAVF2ryRhjUpCfCVkHxPrNRaQucLuPQ15U1W24rjuAaMPVgjdmMqPUFaKqF4UVvS0is4EpwMPAb4s4bgwwBiA7Ozs2ibprV3jzTVi1Clq1iskpjTGmoolp15qIVFXVPB+H1AUe8rH/ZGAbrrsNIIPIEX/eUDZKNfWRqr4lIo8B4Qkrvk45xb3+5z+WlIwxKavEAx1EZII3511R9W2BL/28uaquVFXxsf3kHRrsVmsa5bRNcd2MZbkntBIov+eUADp2dEPD7b6SMSaF+Rl91x+YLyK/DK8Qkf7AfOD4WAV2BMEVbqM9rNsNWFrCQQ5FaQtsKsPx/lWrBp07u5aSMcakKD9J6ULcsO+vRORWABGpISITgQm4h1A7xzzC6N7DddsNFpEqwUIRuRhog7snREh5CxE5QUSqhZTVJwoRuQX3IO6MeARerFNOgQUL4ODBcn9rY4xJBiVOSqr6EdAR+CcwUkQ+BBbiFvx7DvgfVV0elygjY8kBHgBOAT4VkYEi8jDwJ1xyHBV2yETcCLvQ7r5rROQ7EXlGRG4RkdtEZDpu2ff/4u9eV2x07Qr79sGiReX+1sYYkwz8rjy7SUTOA2YD5+Pu3dyqqkXN8hA3qvqciGzFPbD7IrAT+DPwhxJ23c0Fzgb6AA1xw8hXAE8BT6rq9rgEXpzQwQ6dOpX72xtjTKKJaslHNHtdXuOBX+FaTB2BKsAtqjopLhEmuezsbJ03L0YrwatCo0Zw8cUwblxszmmMMUlIROaranZ4uZ/Rdz1w3XUX4FojZ+DuIS0GJojIRBGpEauAU5KIay3ZYAdjTIryM9DhM+AgcLqqPgOgqiuA/8Gto9QXiFjWwvjUtSv88APs3JnoSIwxptz5SUrvAJ1UtdCDNKp6SFX/gGtB1YplcCnplFNcN16sugSNMaYC8TP67teqWuSf76r6CRDxDJPxKXSwgzHGpBg/LaUCIpIhIk1FJD203Buqbcqifn1o29ZmdjDGpCRfSUlEOovIZ7hVWVcDp3nljUTk7yLSMw4xpp6uXV1LKT8/0ZEYY0y58jP6Lgv4AjgW9zBqAVXdDFQHro1pdKnqggtgwwa4+253f8kYY1KEn4dnH8FNctoJNxP39WH1fweuilFcqa1vX9d9N3Ik1KsHDzyQ6IiMMaZc+ElKpwNPqOpuEcmIUr+aw4vvmbIQgVGjYMcOePBBqFMHbr010VEZY0zc+UlKAWBHMfW1yxiLCZWWBmPHuueVbrsNGjSA3/wm0VEZY0xc+Rno8F/g5GLqz8bN7mBipWpV+NOf4LTTYNAg2Lgx0REZY0xc+UlKbwH9w0bYKYCI3IV7eDYl57+Lq0AA3njDzR5+992JjsYYY+LKT1J6FvgK+Ag3S7jilrBYBzwNfAK8GvMIDbRrB0OHwpQp8NlniY7GGGPixs+MDgeAc4G7cQvs5eJWmt0CDAF+par2YE283HsvtGkDN98M+/cnOhpjjIkLXw/Pqmqeqo5U1WxVraGqmar6S1V9TlXz4hWkAapXh1degWXL4JlnEh2NMcbERammGTIJcsEFcOWV8Nhj1o1njKmULClVNKNGQfPmcM45ritv165ER2SMMTFjSamiadIEFi6EO++E0aOhQweYMcPmyTPGVAqWlCqizEx47jn417/c15dc4kboPfssbNmS6OiMMabURG3CzzLJzs7WeYlckG//fvjLX+C11+CLLyAjAzp1gmOPdVubNtCokZsRomFDOOYY9+yTMcYkkIjMV9XsiHJLSmWT8KQU6vvvYdw417333//C6tXRZxlv2tQlq7ZtYfhwaNGi3EM1xqS2opJSiee+E5FhwFhV3RDTyEzsnHSS69YL2r8f1qyBnBzXrZeTA+vWuYS1fDmMH++S07BhiYvZGGNC+F264iERmQmMBWao6qH4hGViIiPDtYbato1e37ata1UZY0yS8DPQoSsuGZ0O/AVYKyJPisjxcYnMxF9WliUlY0xS8TPN0FxV/R3QGLgOWIabXugHEZktIv1FpHqc4jTx0KmT68rbUdyKJMYYU358DwlX1X2qOlFVzwDa4SZjPRaYAGwQkVe9pdNNssvy/pm+/TaxcRhjjKeszymtAOYDPwAC1AR+C8wXkb+JSOMynt/EU6dO7vXrrxMbhzHGeEqVlETkJBF5HlgPTAVOAEYAbYDmwGPAWcC4GMVp4qFxY/fskt1XMsYkCT9DwmsCVwM3AF2AfGAmMAb4W9iyFQ+KyG7goRjGamJNxLWWrKVkjEkSfoaEbwICwFrc8PCxqrq2mP1XATbwIdllZcHIkXDgAKSnJzoaY0yK89N99wlwCdBaVR8+QkJCVaeqqs2tl+w6dYKDB2Hx4kRHYowxvoaE91LV8G46U9EFR+DZfSVjTBKwlkyqO+44N9O43VcyxiSBEiclEckXkUNH2HaLyGIRGSUiTeIZuImRKlWgY0drKRljkoKfltJE4Dvc80hLgPe8balX9i3wIZAH3Ap8LSJtYhqtiY9OnVxSsoUCjTEJ5jcptQYuVNWTVPVyb2sP/Mqre0VVOwIXA3Vxo/RMssvKgp07YeXKREdijElxfpLSCGC0qs4Mr1DVD4DXgSe87/8GjAfOiUWQJs5sZgdjTJLwk5SycNMKFWU50DHk+6+B+qUJypSzDh3cvSW7r2SMSTA/SWk7xbd8egI7Q76vA9j00xVB9epwwgnWUjLGJJyfGR3eBm4TkdeAkcBPgALHAXcAvYAXQvY/C7AnMiuKTp1g5kx44QWXpKpXh3PPhWOOSXRkxpgU4icp3Y9bqmIgbibw4FCtNNzou4+8fRCRAK777ouYRWri69xzYfJkuP32w2XHHgvz50OdOomLyxiTUkRV/R0gciGHR9sBrMQtjf5BbEOrGLKzs3XevHmJDiM29u07vH3zDVxyCfTqBdOmuclbjTEmRkRkvqpmh5eXqKUkIlWApsBuL/mkZAKq9ILddgBNm8ITT8CQIfDyy/D73yc2NmNMSijpQIdquNF1N8QxFpNs7roLLr7Yvc6dm+hojDEpoERJSVVzgS3AnviGY5JKWhpMmABNmsBVV8H27YmOyBhTyfkZEv4B7l6SSSX168Of/uRmexg9OtHRGGMqOT9JaQjQWETeFJFfeCPsTCo49VQ4+2x49VXIy0t0NMaYSsxPUtqMm7GhP7AQ2BNllvBy+40lIjeJyBQRWeK9t79hhIfPc6GIfCkie0Rkm4hME5HWRz4yxQweDKtXw1//muhIjDGVmJ/nlCbiHpZNFvcCR+Geh6oBNPN7AhG5HHgH+Aa4BzcLxe3Av0QkW1XXxy7cCu7ii6F5czcSr1evREdjjKmkSpyUVHVAHOMojTOB1aqaLyJ/xWdSEpFqwEvAGuB0Vd3tlX8IzAeG4x4UNgBVq8LNN8N998EPP8CJJyY6ImNMJVRhV55V1ZVlXJr9DKAJ8EYwIXnnXQjMAvp4icsE3XgjZGS41pIxxsSBr6QkIlVE5BoRmSwin4hIJ6+8nlfeND5hxkUX7/XfUeq+AmoDx5dfOBVAw4bw61/Dm2/CDptr1xgTe36WQ88E/gFMAC4FzgbqedU7gSeBm2McXzwFl2tfF6UuWFaRkmz5GDwY9uyBiRMTHYkxphLyM9BhOJANXAZ8CWwKVqjqIRF5FzgfGFbSE4pIXdzAgpJ6UVW3+di/OJne6/4odblh+xQiIgPx7je1aNEiRuFUENnZ0K0bvPQSDBrk1mEyxpgY8ZOUrgTGqOp7InJUlPqfgD4+378u8JCP/ScDsUpKe73XjCh1gbB9ClHVMcAYcBOyxiieiuPuu6F3b5g6FX7zm0RHY4ypRPzcU2qCGzpdlL1ALT9v7g1WEB/bT37OfwTB4d7RuuiCZdG69sxll8EvfgGPPAKHDiU6GmNMJeInKW2l+HssJ3H4F31FEJxh9NQodd1w98mWlV84FUhaGjz4ICxd6lpLxhgTI36S0t+B67wBD4V4MyBcD8yMVWCxJCKNReSEsNj/AWwAbhSRmiH7/hL3DNQ0VT1YvpFWIJdfDh06wKOPWmvJGBMzfpLSw7jRdnNxo+wUuEBEngAW4AYMPBHzCIsgIheLyDARGQa09cqGedvgsN2fAH4ATgkWeAnnNqA58IWIDBKRPwAfAzn4u9eVeoKtpSVL4M9/TnQ0xphKwtfKsyJyMjAO+EVY1SKgv6oWd88ppkRkAnBtEdWrVLVVlH3PUtVZYef5FW7EYEdcYv07MFRV/1uSOCrVyrN+5edDx46upbRokY3EM8aUWFErz/peDt07WQfgRECAH1X167KHWDGldFICt1T6VVe5wQ+nngrt28Nxx0HdulCzplvJ1pZSN8aEKdNy6OFUdRGudWRS3RVXwHXXwYcfwvTpkfVpaVCrFtSp47batd1URVWrQrVqbktPd2UZGRAIQGbm4aXZg+UZGe77GjXcVrOmO1/9+lCvnjuPMabCK1VS8gYMHIVrKRWiqqvLGpSpQNLSYNw49/W2bW6y1uXLYdeuwtuOHYe3/fvdrBB5eXDggPs+uOXmwr597ms/atVyK+Q2bQrNmkGbNm7Y+i9+4b62rkVjKoQSd9+JSBpuob/fA8cUtZ+qptT//pTvvouX/HyXoEIT1r59Lpnt2QO7d7sEt20b/Pwz5OTA+vWwdi2sW+degz/bmZnQubObiaJbNzjlFJe4rFvRmISJRffdk8DdwPfAX3DPLRkTH2lpLplkRp3p6cj27IHFi+G77+Dbb+E//4EXX4Rnn3X19eu7QRodO7rWVWam6xbMzHRdjLVru9ZXcKtZ03UhWiIzJq78JKV+wExVvTBewaSaM888M6LsqquuYtCgQezdu5cLL4y81AMGDGDAgAFs2bKF3r17R9TffPPN9OnThzVr1tC/f/+I+rvuuouLL76YpUuXctNNN0XUDxs2jJ49e7Jw4UJuvz1yWsLHH3+c7t278+WXX3LfffdF1I8aNYqsrCw+/fRTRowYEVE/evRo2rVrx4wZM3juueci6idNmkTz5s2ZOnUqf/zjHyPq33nnHRo0aMCECROYMGFCRP0HH3xAZmYmr775Jn8OHaqekQFduzLrmWdg7lyenTyZv86fD7Nnu1YZUB340Nv9UdwwzFBHAX/JzIRq1bj3wAH+nZfnkpSXqJpVrcrkRo0AuH3rVhYeOFDo+OOrVWOMVz8wJ4dlofUiZGVkMKphQwD6bdzI2rCl508NBHiiYUMQ4Yp169ga9nzYOTVq8ECDBgD875o17MvPL5REf1WzJncf5WYIO3PVqohrd1WtWgyqX5+9+flcuGZNRP2AOnUYULcuW/Ly6L0ucrKTm+vVo0/t2qw5eJD+6yOfo7+rfn0urlWLpfv3c9PGjRH1wxo0oGeNGizMzeX2TZsi6h9v2JDumZl8uXcv9+XkRNSPOvposgIBPt2zhxFbtkTUjz7mGNplZDBj1y6e2xY5W9mkJk1oXq0aU3fu5I8//xxR/07TpjSoWpUJ27czIcos+R80b05mWhqvbtvGn3ftiqif1bIlAM9u3cpfd+8uVFddhA+9eTQfzcnh73sLz3B2VJUq/KWZWzLu3s2b+fe+fYXqm1WtyuSmbm6D2zduZGFY9/fx6emMadwYgIEbNhT+2QOyAgFGHeM6wPqtW8fag4Uf0Ty1enWeOPpoAK5Ys4athw4xq1UrmDHDdY/HkJ+kVA94L6bvbkx5SkuDrl3dlpt7eGn3/Hw3rL1aNXj1VXcP7NVXYd48Vx7c0tPdCrwHD7pktnHj4S5CVdfKOu009/3cua5rMVTt2u69VV3LLfwXV926cPLJ7usvv3TdlaHd6w0auJYduJZg+H23Ro3gpJPc19u3Rz7UfMwx0K6d+zrKL22aNIG2bd29vii/tGnWDFq1cu+7fXtkffPmbtu7N/rSJi1buvcI3mcM16oVHH20O3fYL23A/fJr0MDFvjfKtJRt27pruGmTu3bhjj/etXrXr49+z7JdO/dvuGaNu9cZ7sQT3R83K1e6axSufXs3gOenn1z3cbgOHdzr0qWwYUPhuipVDtcvXgybNxeuT08/XP/dd7A1rKOqevXD9Xl5kf8+NWsers/NPXx9gz9fdeu6+MHVhV/fo446XB+8L9y+vbseMebnntIc4ANVHR7zKCowu6dkjDH+FXVPye+MDr8TkeaxC8sYY4w5zE/33cnAKmCxiEwHVgDhk56pqj4aq+CMMcakFr+L/AX1K2Ifxd0jNsYYY3zzk5Raxy0KY4wxBh9JSVUjx5AaY4wxMeRnoEOxRCRTRGI7YN0YY0xKKTYpicgBEfl1yPe1ROR9EQlfugLgMuDHWAdojDEmdRyppVQ1bJ904FdAw7hFZIwxJmWVapZwc9j8+fO3iEhp77c1AKI8Wm9C2DUqnl2f4tn1ObJEXaOW0QotKZWRqpa61Sgi86I90WwOs2tUPLs+xbPrc2TJdo1iNtDBGGOMKStLSsYYY5JGSbrvLhSR4KJ+mbhZG64Ukayw/U6OaWSpYUyiA6gA7BoVz65P8ez6HFlSXaNiZwkXkXyf59NUW3nWGGNM7ByppXRWuURhjDHG4GM9JWOMMSbebKBDORORNBG5Q0SWiEiuiKwRkedEpEaiYytPInK8iDwiIl+JSI6I7BKRhSJyf7RrISLtROT/RORnEdkjIl+IyNmJiD1RvKm8louIisjLUepT7hqJSH0ReVZEfvL+P+WIyOcicnrYfl1F5FPv52yniMyMcl+80hGRmiJyn4h85332LSLypYgMEBEJ2zcprpE9p1T+RgK3AtOB54ATve87iUhPVfV7H6+iuh64BXgfmAIcxHUXjwCuEpFuqroPQESOBb4E8oCngR3Ab4GPROR/VfXTBMSfCI9QxGwqqXiNRKQlMAuoCYwFlgF1gI5A05D9unn7rQMe9IoHA1+ISHdV/a78oi4/IpIGfAh0B94EXsINVrsaGI/73TPU2zd5rpGq2lZOG3ASkA/8Jaz897hRjb9JdIzleC2ygTpRykd412JwSNmfcQtKZoWU1cQtOrkUrxu6Mm9AZ1zCudO7Pi+H1afcNQK+ANYAjY+w3xxgJ9A0pKypV/Zxoj9HHK/Pqd7Pysiw8nRgObA9Ga+Rdd+Vr6sBAUaFlb8O7KXoxRMrHVWdp6o7olRN9V47AHhdeZcAs1R1Ycjxu4E3gOOBLnEON6FEpAruZ2Qm8G6U+pS7RiLSAzgNeFpVN4hINRHJjLJfW9xnn6aq64Ll3tfTgJ4hj7xUNrW91/Whhap6ADet0B5IvmtkSal8dcG1lOaEFqpqLrCQSvaLo5Saea+bvNeOQAbw7yj7fuW9VvbrdgdwAq47JZpUvEYXeq+rRWQGsA/YIyLLRCT0j7vg5y7q2giV9xnLOcB2YIiIXCkiLUTkBBF5AveZh3v7JdU1sqRUvpoAW1R1f5S6dUADEUkv55iShtcieADXTfWWV9zEe10X5ZBgWdModZWCiLQGHgYeUdWVReyWiteonff6OlAfuBZ3n/IAMElErvPqU/HaAKCqP+Na0Ntw3burgB9w93KvUNXXvV2T6hrZQIfylQlES0gAuSH7HCifcJLOKFw/+H2qutQrC3bJRLtuuWH7VEav4fr/ny9mn1S8RrW8113AWV6XFCLyf7jr9biIvElqXptQu4FFuAFFX+IS+C3AWyJyqap+QpJdI0tK5Wsv0KiIukDIPilHRB7FdU+NUdUnQqqC1yMjymGV+pp53VDnAj1U9WAxu6biNdrnvf4pmJDAtQ5E5H3gGlxrKhWvDQDeYqxfAneo6msh5X/CJarXvVGbSXWNrPuufK3HddFF+8dviuvaS7lWkogMB4bhhqn+Lqw6eJM2WvdBsCxat0OF5v2MPA98AGwUkbbeDengGjR1vLK6pOY1Wuu9boxSt8F7rUdqXpugO3BJZVpooaruBf6G+1lqRZJdI0tK5Wsu7pqfElooIgEgC5iXiKASyUtID+Geo7hRvbGoIb7DdSucGuXwbt5rZbxu1XHPJF0E/BiyzfLq+3nf30hqXqPgYKFmUeqCZZtx/+eg6GujwPzYhpY0ggkl2nykVUNek+saJXosfSptwC8o/jmlfomOsZyvx4Pe554IpBWz3zTcMzi/DCkLPoOzjMr5DE41oHeU7Wbvmn3ofX98Kl4jXCtoJ67FVDOkvDHuPsrSkLK53r5NQsqaeGWfJvqzxPEajfR+VoaElQdb19uAKsl2jWzuu3ImIi/h7p1Mx3XNBGd0+BdwtqbIjA4icgvwMrAaN+Iu/HNvUncTNvgcxRzcrA8jcf9RfotL8hep6kflFXeiiUgrYAXwiqoODilPuWskIgOB0cD3wDjcQ6E34xLTr1T1Y2+/7sDnuAT2knf474Gjgf9R1W/KOfRy4c14sQCXwKfgfsfUx/1ctAJuUdVXvX2T5xolOpun2oZrSt+Fe8p+P66v9nlC/tpLhQ2YgPsrrqhtVtj+JwLv4Z672Av8E+iZ6M+RgOvWiigzOqTqNQIuxz1Lswc3Eu9j3C/R8P1OBf6Oa0XtAj4COic6/nK4PsfiusbX4v5g2QnMBi5P1mtkLSVjjDFJwwY6GGOMSRqWlIwxxiQNS0rGGGOShiUlY4wxScOSkjHGmKRhSckYY0zSsKRkjDEmaVhSMsYAICJnioiKyIBEx2JSlyUlY2Ik5Jf63d73dUVkuIicmeDQCohIlhdTq0THYkw0tp6SMfFTFzcDOhye3TvRsnAxzQJWhtXNxs1OXtzaTcbElbWUjKmgRKTWkfcqOVXNV9VcVT0Uy/Ma44clJWPiwOuyW+F9+5DXracisjJsvz4i8k8R2SUie0XkPyLSO8r5VEQmiMg53v67gRleXRMReU5EForIzyKSKyKLRWSoiFQJOcdw3EKKAJ+HxDQhGHO0e0oiUkNEnhCR/4rIfhHZKCITvVmoC33m4PEicp2IfO/tv0pEhpT+appUYt13xsTHD7iVP0filil51yvfHdxBREYA9wMzObx8x2XANBEZrKqvhJ0zG7gCeB0383NQR9xs2dOB/+LWYroAeBJoA9zk7fcublmHgcDjXox4x0QlItVws0X/D/AO8BxwHG6JiPNEJFtV14Yd9jvckgdjcTOW9wOeEpG1qvpWUe9lDGBLV9hmW6w24EzcshJ3e9+38r4fHmXfzl7d41Hq/g+3xECtkLLgkh4RS1Hg7gNFLOIHTMIt/Nc4pGyAd54zi4l/QEjZb72yp8P2vcgrnxTl+PVAnZDyTCAH+Hei/41sS/7Nuu+MSYy+uF/gb4pIg9ANeB+oReTy1N+o6qfhJ1LVfaqqACKSLiL1vfN8hOuizy5DnJfhWnBPhL3n34CFwKUiEv57ZLyq7gjZdy9uzaPjyhCHSRHWfWdMYpwICLCkmH2ODvt+WbSdRKQq8AfgGqCtd95Q9UoZI0BrYL2q/hyl7nvcaL4GwOaQ8uVR9t0KHFWGOEyKsKRkTGIIrqX0v7gutmi+D/t+bxH7PY9bunoq8BguQRzEdRE+RfkPaLLRe6bULCkZEz/FLev8I24wwmpV/aGY/UqiPzBbVX8dWigibX3GFM1y4AIRqauq28Pq2uPufW3xeU5jimT3lIyJn+BIu/pR6iZ5r4+HDtsOEpHwrrviHCKsy05EauBG//mJKZr/w/2e+EPY+f8X6AS8r6r5PmI1pljWUjImTlR1q4j8BPxaRP4LbAL2qOoMVZ3rPTc0HFgoItNwo9YaAycDFwLpJXyrd4CbRGQq8CnuXtT1uPs44ebiBi7cLyL1gD3AClX9TxHnngBcCwz1piaajbtvNcj7PPeVMEZjSsSSkjHx1Rf3rNLjuKHRq/AeelXVh0VkHnArcDtQA3c/aJFXVlJ3AruAq4BLgTXAGFwCKjRaT1VXi8j1wFDgj7hnmt4EoiYlVT0oIucDw4A+uOehtgPTgGGqusZHnMYckXgjSY0xxpiEs3tKxhhjkoYlJWOMMUnDkpIxxpikYUnJGGNM0rCkZIwxJmlYUjLGGJM0LCkZY4xJGpaUjDHGJA1LSsYYY5KGJSVjjDFJ4/8BmIx1g4YbhioAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"font.size\"] = 18\n",
    "plt.plot(cost_history, color=\"red\", label=\"VQE\")\n",
    "plt.plot(range(len(cost_history)), [molecule.fci_energy]*len(cost_history), linestyle=\"dashed\", color=\"black\", label=\"Exact Solution\")\n",
    "plt.xlabel(\"Iteration\")\n",
    "plt.ylabel(\"Energy expectation value\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are interested, you can calculate the ground state by varying the `distance` between hydrogen atoms to find the interatomic distance at which the hydrogen molecule is most stable. (It should be about 0.74 angstroms, depending on the performance of the ansatz.)"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}