challenge-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f9369076",
   "metadata": {
    "tags": []
   },
   "source": [
    "## IBM Quantum Challenge Fall 2021\n",
    "# Challenge 2: Calculate bandgap of OLED molecules\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "We recommend that you switch to **light** workspace theme under the Account menu in the upper right corner for optimal experience."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d6d42da",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Organic Light Emitting Diodes or OLEDs have become increasingly popular in recent years as the basis for fabrication of thin, flexible TV and mobile phone displays that emit light upon application of an electric current. Recent studies ([**Gao et al., 2021**](https://www.nature.com/articles/s41524-021-00540-6)) have been looking at electronic transitions of high energy states in phenylsulfonyl-carbazole (PSPCz) molecules, which could be useful thermally activated delayed fluorescence (TADF) emitters for OLED technology. TADF emitters could potentially produce OLEDs that perform with 100 percent internal quantum efficiency (IQE), i.e the fraction of the charge carriers in a circuit or system that emit absorbed photons, compared with conventional fluorophores currently used to make OLEDs whose quantum efficiencies are limited to 25 percent. That large boost in efficiency means manufacturers could produce OLEDs for use in devices requiring low-power consumption, such as cell phones, which could in turn lead to future developments where virtually any surface can be converted into a cheap and energy-efficient lighting source covering vast areas of homes, offices, museums and more!\n",
    "\n",
    "<center><img src=\"resources/JSR_img6_1920w.jpg\" width=\"600\"></center>\n",
    "\n",
    "### Why quantum?\n",
    "\n",
    "Quantum computers could be invaluable tools for studying the electronic structure and dynamical properties of complex molecules and materials as it makes more sense to model quantum mechanical systems on a quantum device than on a classical computer. A recent joint research project by IBM Quantum and partners was successful in developing methods to improve accuracy for the calculation of excited TADF states for efficient OLEDs, making it the world's first research case of applying quantum computers to the calculation of excited states of commercial materials (see paper linked above for reference). \n",
    "    \n",
    "With this background information, we are interested in describing quantum computations of the “excited states,” or high energy states, of industrial chemical compounds that could potentially be used in the fabrication of efficient OLED devices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddc39149",
   "metadata": {},
   "source": [
    "## Challenge\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Goal**\n",
    "\n",
    "The goal of this challenge is to use quantum algorithms to reliably predict the excited states energies of these TADF materials. Along the way, this challenge introduces state-of-the-art hybrid classical-quantum embedded chemistry modelling allowing the splitting of the work-load between classical approximations and more accurate quantum calculations. \n",
    "    \n",
    "1. **Challenge 2a & 2b**: Understanding the atomic orbitals (AO), molecular orbitals (MO) and how to reduce the number of orbitals using active space transformation.\n",
    "2. **Challenge 2c & 2d**: Calculating ground state energy of PSPCz molecule using NumPy and Variational Quantum Eigensolver (VQE).\n",
    "3. **Challenge 2e**: Calculating excited state energy of PSPCz module using quantum Equation-of-Motion (QEOM) algorithm.\n",
    "4. **Challenge 2f**: Running VQE on the cloud (simulator or real quantum system) using Qiskit Runtime.\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "Before you begin, we recommend watching the [**Qiskit Nature Demo Session with Max Rossmannek**](https://youtu.be/UtMVoGXlz04?t=38) and check out the corresponding [**demo notebook**](https://github.com/qiskit-community/qiskit-application-modules-demo-sessions/tree/main/qiskit-nature) to learn how to define electronic structure calculations.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f81686b",
   "metadata": {},
   "source": [
    "### 1. Driver\n",
    "\n",
    "The interfaces to the classical chemistry codes that are available in Qiskit are called drivers. We have for example `PSI4Driver`, `PyQuanteDriver`, `PySCFDriver` are available.\n",
    "\n",
    "By running a driver (Hartree-Fock calculation for a given basis set and molecular geometry), in the cell below, we obtain all the necessary information about our molecule to apply then a quantum algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e15447fe",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/pyscf/lib/misc.py:47: H5pyDeprecationWarning: Using default_file_mode other than 'r' is deprecated. Pass the mode to h5py.File() instead.\n",
      "  h5py.get_config().default_file_mode = 'a'\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.drivers import Molecule\n",
    "from qiskit_nature.drivers.second_quantization import ElectronicStructureDriverType, ElectronicStructureMoleculeDriver\n",
    "\n",
    "# PSPCz molecule\n",
    "geometry = [['C', [  -0.2316640,    1.1348450,    0.6956120]],\n",
    "            ['C', [  -0.8886300,    0.3253780,   -0.2344140]],\n",
    "            ['C', [  -0.1842470,   -0.1935670,   -1.3239330]],\n",
    "            ['C', [   1.1662930,    0.0801450,   -1.4737160]],\n",
    "            ['C', [   1.8089230,    0.8832220,   -0.5383540]],\n",
    "            ['C', [   1.1155860,    1.4218050,    0.5392780]],\n",
    "            ['S', [   3.5450920,    1.2449890,   -0.7349240]],\n",
    "            ['O', [   3.8606900,    1.0881590,   -2.1541690]],\n",
    "            ['C', [   4.3889120,   -0.0620730,    0.1436780]],\n",
    "            ['O', [   3.8088290,    2.4916780,   -0.0174650]],\n",
    "            ['C', [   4.6830900,    0.1064460,    1.4918230]],\n",
    "            ['C', [   5.3364470,   -0.9144080,    2.1705280]],\n",
    "            ['C', [   5.6895490,   -2.0818670,    1.5007820]],\n",
    "            ['C', [   5.4000540,   -2.2323130,    0.1481350]],\n",
    "            ['C', [   4.7467230,   -1.2180160,   -0.5404770]],\n",
    "            ['N', [  -2.2589180,    0.0399120,   -0.0793330]],\n",
    "            ['C', [  -2.8394600,   -1.2343990,   -0.1494160]],\n",
    "            ['C', [  -4.2635450,   -1.0769890,    0.0660760]],\n",
    "            ['C', [  -4.5212550,    0.2638010,    0.2662190]],\n",
    "            ['C', [  -3.2669630,    0.9823890,    0.1722720]],\n",
    "            ['C', [  -2.2678900,   -2.4598950,   -0.3287380]],\n",
    "            ['C', [  -3.1299420,   -3.6058560,   -0.3236210]],\n",
    "            ['C', [  -4.5179520,   -3.4797390,   -0.1395160]],\n",
    "            ['C', [  -5.1056310,   -2.2512990,    0.0536940]],\n",
    "            ['C', [  -5.7352450,    1.0074800,    0.5140960]],\n",
    "            ['C', [  -5.6563790,    2.3761270,    0.6274610]],\n",
    "            ['C', [  -4.4287740,    3.0501460,    0.5083650]],\n",
    "            ['C', [  -3.2040560,    2.3409470,    0.2746950]],\n",
    "            ['H', [  -0.7813570,    1.5286610,    1.5426490]],\n",
    "            ['H', [  -0.7079140,   -0.7911480,   -2.0611600]],\n",
    "            ['H', [   1.7161320,   -0.2933710,   -2.3302930]],\n",
    "            ['H', [   1.6308220,    2.0660550,    1.2427990]],\n",
    "            ['H', [   4.4214900,    1.0345500,    1.9875450]],\n",
    "            ['H', [   5.5773000,   -0.7951290,    3.2218590]],\n",
    "            ['H', [   6.2017810,   -2.8762260,    2.0345740]],\n",
    "            ['H', [   5.6906680,   -3.1381740,   -0.3739110]],\n",
    "            ['H', [   4.5337010,   -1.3031330,   -1.6001680]],\n",
    "            ['H', [  -1.1998460,   -2.5827750,   -0.4596910]],\n",
    "            ['H', [  -2.6937370,   -4.5881470,   -0.4657540]],\n",
    "            ['H', [  -5.1332290,   -4.3740010,   -0.1501080]],\n",
    "            ['H', [  -6.1752900,   -2.1516170,    0.1987120]],\n",
    "            ['H', [  -6.6812260,    0.4853900,    0.6017680]],\n",
    "            ['H', [  -6.5574610,    2.9529350,    0.8109620]],\n",
    "            ['H', [  -4.3980410,    4.1305040,    0.5929440]],\n",
    "            ['H', [  -2.2726630,    2.8838620,    0.1712760]]]\n",
    "\n",
    "molecule = Molecule(geometry=geometry, charge=0, multiplicity=1)\n",
    "\n",
    "driver = ElectronicStructureMoleculeDriver(molecule=molecule, \n",
    "                                           basis='631g*', \n",
    "                                           driver_type=ElectronicStructureDriverType.PYSCF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05144bc5",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    \n",
    "**Challenge 2a**\n",
    "    \n",
    "Question: Find out these numbers for the PSPCz molecule.\n",
    "    \n",
    "1. What is the number of C, H, N, O, S atoms?\n",
    "1. What is the total number of atoms?\n",
    "1. What is the total number of atomic orbitals (AO)?\n",
    "1. What is the total number of molecular orbitals (MO)?\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "**How to count atomic orbitals?**\n",
    "\n",
    "The number depends on the basis. The number below is specific to `631g*` basis which we will use for this challenge. \n",
    "\n",
    "- C: 1s, 2s2p, 3s3p3d  = 1+4+9 = 14\n",
    "\n",
    "- H: 1s, 2s  = 1+1 = 2\n",
    "\n",
    "- N: 1s, 2s2p, 3s3p3d  = 1+4+9 = 14\n",
    "\n",
    "- O: 1s, 2s2p, 3s3p3d  = 1+4+9 = 14\n",
    "\n",
    "- S: 1s, 2s2p, 3s3p3d, 4s4p = 1+4+9+4 = 18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a32b91c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'C': 24, 'H': 17, 'N': 1, 'O': 2, 'S': 1, 'atoms': 45, 'AOs': 430, 'MOs': 430}\n"
     ]
    }
   ],
   "source": [
    "num_ao = {\n",
    "    'C': 14,\n",
    "    'H': 2,\n",
    "    'N': 14,\n",
    "    'O': 14,\n",
    "    'S': 18,\n",
    "}\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "num_C_atom = 24\n",
    "num_H_atom = 17\n",
    "num_N_atom = 1\n",
    "num_O_atom = 2\n",
    "num_S_atom = 1\n",
    "num_atoms_total = 45\n",
    "num_AO_total = 430\n",
    "num_MO_total = 430\n",
    "\n",
    "##############################\n",
    "\n",
    "answer_ex2a ={\n",
    "    'C': num_C_atom,\n",
    "    'H': num_H_atom,\n",
    "    'N': num_N_atom,\n",
    "    'O': num_O_atom,\n",
    "    'S': num_S_atom,\n",
    "    'atoms': num_atoms_total,\n",
    "    'AOs': num_AO_total,\n",
    "    'MOs': num_MO_total\n",
    "}\n",
    "\n",
    "print(answer_ex2a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d85ff593",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitting your answer for 2a. Please wait...\n",
      "Congratulations 🎉! Your answer is correct and has been submitted.\n"
     ]
    }
   ],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2a\n",
    "grade_ex2a(answer_ex2a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98047d20",
   "metadata": {},
   "source": [
    "As you found out yourself in the exercise above, PSPCz is a large molecule, consisting of many atoms and many atomic orbitals. Direct calculation of a large molecule is out of reach for current quantum systems. However, since we are only interested in the bandgap, calculating the energy of Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) is sufficient. Here we applied a technique called active space transformation to reduce the number of molecular orbitals to only 2 (HOMO and LUMO):\n",
    "\n",
    "$$E_g = E_{LUMO} - E_{HOMO}$$\n",
    "\n",
    "\n",
    "<center><img src=\"resources/Molecule_HOMO-LUMO_diagram.svg\" width=\"600\"></center>\n",
    "\n",
    "\n",
    "Each circle here represents an electron in an orbital; when light or energy of a high enough frequency is absorbed by an electron in the HOMO, it jumps to the LUMO.\n",
    "\n",
    "For PSPCz molecules, we limit this excited state to just the first singlet and triplet states. In a singlet state, all electrons in a system are spin paired, giving them only one possible orientation in space. A singlet or triplet excited state can form by exciting one of the two electrons to a higher energy level. The excited electron retains the same spin orientation in a singlet excited state, whereas in a triplet excited state, the excited electron has the same spin orientation as the ground state electron. \n",
    "\n",
    "<center><img src=\"resources/spin.jpg\" width=\"300\"><figcaption>Spin in the ground and excited states</figcaption></center>\n",
    "\n",
    "One set of electron spins is unpaired in a triplet state, meaning there are three possible orientations in space with respect to the axis. LUMO (a-c) and HOMO (e-f) orbitals of the triplet state optimized structures of PSPCz (a, d) and its variants 2F-PSPCz (b, e) and 4F-PSPCz (c, f) respectively would then look something like this.\n",
    "\n",
    "\n",
    "<center><img src=\"resources/oled_paper_fig2.jpg\" width=\"600\"></center>\n",
    "\n",
    "<center><img src=\"resources/oled_paper_fig1.jpg\" width=\"600\"></center>\n",
    "\n",
    "\n",
    "By using the active space transformer method, we will manage to exclude non-core electronic states by restricting calculations to the singlet and triplet, i.e. the smallest possible active space and manage to compute this energy with a small number of qubits while keeping a high-quality description of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b4c60b72",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers.second_quantization import HDF5Driver\n",
    "\n",
    "driver_reduced = HDF5Driver(\"resources/PSPCz_reduced.hdf5\")\n",
    "properties = driver_reduced.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a63297f8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ElectronicEnergy\n",
      "\t(AO) 1-Body Terms:\n",
      "\t\tAlpha\n",
      "\t\t<(430, 430) matrix with 184900 non-zero entries>\n",
      "\t\t[0, 0] = -11.481107571585675\n",
      "\t\t[0, 1] = -2.6982522446048134\n",
      "\t\t[0, 2] = -2.237143188610541\n",
      "\t\t[0, 3] = 0.0017433998087159669\n",
      "\t\t[0, 4] = 0.0007741436199762753\n",
      "\t\t... skipping 184895 entries\n",
      "\t\tBeta\n",
      "\t\t<(430, 430) matrix with 184900 non-zero entries>\n",
      "\t\t[0, 0] = -11.481107571585675\n",
      "\t\t[0, 1] = -2.6982522446048134\n",
      "\t\t[0, 2] = -2.237143188610541\n",
      "\t\t[0, 3] = 0.0017433998087159669\n",
      "\t\t[0, 4] = 0.0007741436199762753\n",
      "\t\t... skipping 184895 entries\n",
      "\t(MO) 1-Body Terms:\n",
      "\t\tAlpha\n",
      "\t\t<(2, 2) matrix with 4 non-zero entries>\n",
      "\t\t[0, 0] = -0.4968112637934733\n",
      "\t\t[0, 1] = 0.00027750088691888997\n",
      "\t\t[1, 0] = 0.00027750088691825913\n",
      "\t\t[1, 1] = -0.1843594001763901\n",
      "\t\tBeta\n",
      "\t\t<(2, 2) matrix with 4 non-zero entries>\n",
      "\t\t[0, 0] = -0.4968112637934733\n",
      "\t\t[0, 1] = 0.00027750088691888997\n",
      "\t\t[1, 0] = 0.00027750088691825913\n",
      "\t\t[1, 1] = -0.1843594001763901\n",
      "\t(MO) 2-Body Terms:\n",
      "\t\tAlpha-Alpha\n",
      "\t\t<(2, 2, 2, 2) matrix with 16 non-zero entries>\n",
      "\t\t[0, 0, 0, 0] = 0.22795982746869856\n",
      "\t\t[0, 0, 0, 1] = -0.00027753808830176344\n",
      "\t\t[0, 0, 1, 0] = -0.00027753808830176615\n",
      "\t\t[0, 0, 1, 1] = 0.13689436105642472\n",
      "\t\t[0, 1, 0, 0] = -0.0002775380883017597\n",
      "\t\t... skipping 11 entries\n",
      "\t\tBeta-Alpha\n",
      "\t\t<(2, 2, 2, 2) matrix with 16 non-zero entries>\n",
      "\t\t[0, 0, 0, 0] = 0.22795982746869856\n",
      "\t\t[0, 0, 0, 1] = -0.00027753808830176344\n",
      "\t\t[0, 0, 1, 0] = -0.00027753808830176615\n",
      "\t\t[0, 0, 1, 1] = 0.13689436105642472\n",
      "\t\t[0, 1, 0, 0] = -0.0002775380883017597\n",
      "\t\t... skipping 11 entries\n",
      "\t\tBeta-Beta\n",
      "\t\t<(2, 2, 2, 2) matrix with 16 non-zero entries>\n",
      "\t\t[0, 0, 0, 0] = 0.22795982746869856\n",
      "\t\t[0, 0, 0, 1] = -0.00027753808830176344\n",
      "\t\t[0, 0, 1, 0] = -0.00027753808830176615\n",
      "\t\t[0, 0, 1, 1] = 0.13689436105642472\n",
      "\t\t[0, 1, 0, 0] = -0.0002775380883017597\n",
      "\t\t... skipping 11 entries\n",
      "\t\tAlpha-Beta\n",
      "\t\t<(2, 2, 2, 2) matrix with 16 non-zero entries>\n",
      "\t\t[0, 0, 0, 0] = 0.22795982746869856\n",
      "\t\t[0, 0, 0, 1] = -0.00027753808830176344\n",
      "\t\t[0, 0, 1, 0] = -0.00027753808830176615\n",
      "\t\t[0, 0, 1, 1] = 0.13689436105642472\n",
      "\t\t[0, 1, 0, 0] = -0.0002775380883017597\n",
      "\t\t... skipping 11 entries\n",
      "\tEnergy Shifts:\n",
      "\t\tActiveSpaceTransformer = -4042.866322560092\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.properties.second_quantization.electronic import ElectronicEnergy\n",
    "\n",
    "electronic_energy = properties.get_property(ElectronicEnergy)\n",
    "print(electronic_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "292393c5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-18T02:17:23.308949Z",
     "iopub.status.busy": "2021-09-18T02:17:23.308711Z",
     "iopub.status.idle": "2021-09-18T02:17:23.321441Z",
     "shell.execute_reply": "2021-09-18T02:17:23.320065Z",
     "shell.execute_reply.started": "2021-09-18T02:17:23.308925Z"
    },
    "tags": []
   },
   "source": [
    "You can see that `(AO) 1-Body Terms` contains a (430 x 430) matrix which describes the original molecule with 430 atomic orbitals which translate to 430 molecular orbitals (?). After `ActiveSpaceTransformation` (pre-calculated), the number of molecular orbitals `(MO) 1-Body Terms` is reduced to a (2x2) matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b9f88d1",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    \n",
    "**Challenge 2b**\n",
    "    \n",
    "Question: Use property framework to find out the answer for the questions below. \n",
    "    \n",
    "1. What is the number of electrons in the system after active space transformation?\n",
    "1. What is the number of molecular orbitals (MO)?\n",
    "1. What is the number of spin orbitals (SO)?\n",
    "1. How many qubits would you need to simulate this molecule with Jordan-Wigner mapping?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5f9f2668",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'electrons': 2, 'MOs': 2, 'SOs': 4, 'qubits': 4}\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.properties.second_quantization.electronic import ParticleNumber\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "particle_number = properties.get_property(ParticleNumber)\n",
    "num_electron = 2 \n",
    "num_MO = 2 \n",
    "num_SO = 4 \n",
    "num_qubits = 4 \n",
    " \n",
    "\n",
    "##############################\n",
    "\n",
    "answer_ex2b = {\n",
    "    'electrons': num_electron,\n",
    "    'MOs': num_MO,\n",
    "    'SOs': num_SO,\n",
    "    'qubits': num_qubits\n",
    "}\n",
    "\n",
    "print(answer_ex2b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "58cd92be",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitting your answer for 2b. Please wait...\n",
      "Congratulations 🎉! Your answer is correct and has been submitted.\n"
     ]
    }
   ],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2b\n",
    "grade_ex2b(answer_ex2b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93e47c9",
   "metadata": {},
   "source": [
    "### 2. Electronic structure problem\n",
    "\n",
    "You can then create an ElectronicStructureProblem that can produce the list of fermionic operators before mapping them to qubits (Pauli strings). This is the first step in defining your molecular system in its ground state. You can read more about solving for the ground state in [**this tutorial**](https://qiskit.org/documentation/nature/tutorials/03_ground_state_solvers.html).\n",
    "\n",
    "<center><img src=\"resources/H2_gs.png\" width=\"300\"></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6c6f3b58",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fermionic Operator\n",
      "register length=4, number terms=26\n",
      "  (0.01572205126528473+0j) * ( +_0 -_1 +_2 -_3 )\n",
      "+ (-0.01572205126528473+0j) * ( +_0 -_1 -_2 +_3 )\n",
      "+ (0.00027750088691888997+0j) * ( +_0 -_1 )\n",
      "+ (0.0003149147870892302+0j) * ( +_0 -_1 +_3 -_3 )\n",
      "+  ...\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.problems.second_quantization import ElectronicStructureProblem\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "es_problem = ElectronicStructureProblem(driver_reduced)\n",
    "\n",
    "##############################\n",
    "second_q_op = es_problem.second_q_ops()\n",
    "print(second_q_op[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e17e2c51",
   "metadata": {},
   "source": [
    "### 3. QubitConverter\n",
    "\n",
    "Allows to define the mapping that you will use in the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "71327193",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.45781773131305903 * IIII\n",
      "- 0.009666607989543467 * ZIII\n",
      "+ 0.12689900731767084 * IZII\n",
      "+ 0.030293077447785 * ZZII\n",
      "- 0.009666607989543479 * IIZI\n",
      "+ 0.03732964036584735 * ZIZI\n",
      "+ 0.034223590264106186 * IZZI\n",
      "+ 0.12689900731767084 * IIIZ\n",
      "+ 0.034223590264106186 * ZIIZ\n",
      "+ 0.05698995686717464 * IZIZ\n",
      "+ 0.030293077447785 * IIZZ\n",
      "+ 0.00014809461815615455 * XXII\n",
      "+ 0.00014809461815615455 * YYII\n",
      "- 7.872869677230731e-05 * XXZI\n",
      "- 7.872869677230731e-05 * YYZI\n",
      "+ 6.938452207544002e-05 * XXIZ\n",
      "+ 6.938452207544002e-05 * YYIZ\n",
      "+ 0.00014809461815615455 * IIXX\n",
      "- 7.872869677230731e-05 * ZIXX\n",
      "+ 6.938452207544002e-05 * IZXX\n",
      "+ 0.00014809461815615455 * IIYY\n",
      "- 7.872869677230731e-05 * ZIYY\n",
      "+ 6.938452207544002e-05 * IZYY\n",
      "+ 0.003930512816321183 * XXXX\n",
      "+ 0.003930512816321183 * YYXX\n",
      "+ 0.003930512816321183 * XXYY\n",
      "+ 0.003930512816321183 * YYYY\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.converters.second_quantization import QubitConverter\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper, ParityMapper, BravyiKitaevMapper\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "qubit_converter = QubitConverter(JordanWignerMapper())\n",
    "\n",
    "##############################\n",
    "\n",
    "qubit_op = qubit_converter.convert(second_q_op[0])\n",
    "print(qubit_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a0c1f33",
   "metadata": {},
   "source": [
    "### 4. Initial state\n",
    "\n",
    "A good initial state in chemistry is the HartreeFock state. We can initialize it as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f5b00153",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAH0AAADWCAYAAAD8W40IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJqklEQVR4nO3df0wU+R3G8WdmF0SDCWdIMMpBXXQNoiBqleMSWVtD+IP2SBXRtEQpiXtCm0hqrU2qCcWSVqnSXppom5zchRxNRG3uDJJIFPAXTaiA1SbdgIuEhHoRiBHPoDLbP8wu4mqZpbvMDJ/nlfjHLbnd7+WdGZdFn1N8Pp8PJIpq9AFo9jG6QIwuEKMLxOgCMbpAjC4QowvE6AIxukCMLhCjC8ToAjG6QIwuEKMLxOgCMbpAjC4QowvE6AIxukCMLhCjC8ToAjG6QHajD2B2P/v3v9Dz5Ikhr52xcCF+v3JV2J+X0afR8+QJ2kdHjD5GWPH2LhCjC8ToAjG6QIwuEKMLxOgCMbpAjC6QqaNrmoaamhqsWLECMTExyMjIQFtbG1auXIm9e/cafby38r14gRcf/wQTp/8y5fGJC3/Dix/thm9szKCTTTJ19NLSUlRVVcHtduPSpUvYsWMHdu3ahfv372P9+vVGH++tlKgo2A/9HNrFJmhd3QAAn9cL7dPPYDt4AEpsrLEHhIk/e29oaEBdXR1aW1uRk5MDANiyZQtu376N8+fPY926dQaf8N2UbyVD/fFuTNSchPLJSbz87XGoH30Pavoao48GwMRXenV1NfLy8gLB/ZYvX46oqCikp6cDAPr7+5GTkwOn04k1a9bg2rVrRhw3iFrwEZSk9/HSXQ7YbFB3Fxt9pABTRh8cHMTdu3dRWFgY9LWBgQGkpaVh3rx5AAC3242ioiJ4PB6cPn0aO3fuxPPnz6d9DUVRdP1qbW2d0X+DoihQ0tcAjx9D/e53oERFhfwcra2tus8ZCtNGB4DFixdPefzZs2doa2sL3NofPXqE69evo7S0FACQnZ2NJUuW4OrVq7N74Lfweb3Qvvgr1KJCaPVfwPf110YfKcCU0ePj4wEAHo9nyuPHjh3D0NBQ4E3cwMAAEhISAlc9ACxbtgwPHjyY9jV8Pp+uXy6XK+Tz+56/ePX7+A8KYCstgfLhB5g4fgI+TQvpeVwul+5zhsKUb+QcDgfS09NRXV2NRYsWYenSpWhsbERTUxMAmPadu5/26RkodjvU4h8CAGxlH+OluxzauQuwFW4z+HQmvdJVVcXZs2eRlpaGffv2oaSkBPHx8SgvL4fNZgu8iUtKSsLDhw8xPj4e+He9Xi+Sk5ONOjq0rm5oTc2wHToIxf7qmlIWLIDtFwegfV4Pn9dr2Nn8FCtNfxcXF6Onpwd37twJPJabm4uCggKUlZXh5s2b2L59O/r7+xEdHR2W19za+XfD/rjU5vcWoWXDprA/rylv7+/S2dmJrKysKY+dOnUKe/bsQW1tLaKjo9HQ0BC24HOVZaKPjY3B4/GgrKxsyuMOhwPt7e0GncqaLBM9NjYWExMTRh9jTjDlGzmKLEYXiNEFYnSBGF0gRheI0QWyzPfpRslYuHDOvbalPnun8ODtXSBGF4jRBWJ0gRhdIEYXiNEFYnSBGF0gRheI0QVidIEYXSD+aHUaXIEWiCvQNCcwukCMLhCjC8ToAjG6QIwuEKMLZOroVhwEtgJTR7fiILAVVqBN+zGsVQeB/SvQL3+6H8rGb0PNXDu5Av2bX5tiBdq0V7reQeAjR47A6XRCVVU0NjYacdQgr69A+0ZGuAKtRyiDwHl5eWhubsbmzZtn+5j/E1egQ6R3EBh4NQLscDhCfg2uQJuM3kFgMzPzCrQp38jNxiCw3r+hPZOZ0Ckr0CW74RsdxcTxE7D9rhqKqv86c7lcaInA3yQ35ZWudxDYrN62Au37z0No5y4YfLJXTHmlA4DT6Qwa6y8uLsaqVaswf/58g041Pf8KtP1PfwxagZ745a+gblgHZdkyQ89oyiv9XTo7O4Nu7YcPH0ZiYiJu3boFt9uNxMRE9PX1GXRCQM1ci6gvz0N5P3Hq46vTEPXVBcODAxaK7h8EfvNDmaqqKgwODmJ8fBzDw8MYHBxESkqKQae0BtPe3t/EQeDwscyVTuHD6AIxukCMLhCjC8ToAjG6QJb5Pt0oHASmOYG3d4EYXSBGF4jRBWJ0gRhdIEYXiNEFYnSBGF0gRheI0QVidIH4o9VpcBBYIA4C05zA6AIxukCMLhCjC8ToAjG6QIwukKmjcxA4MkwdnYPAkWHaj2E5CBw5pr3S9QwCj46OIj8/H06nExkZGcjNzUVvb69BJ57EQeAZ0DsIrCgK9u/fD4/Hg56eHuTn56OkpMSAEwfjIHCI9A4Cx8XFYevWrYGvZ2dnw+v16noNDgKbzEwHgWtra1FQUBDp4+nCQeAQzWQQuLKyEr29vbhy5Yqu1+AgsMmEOgh89OhRXLx4Ec3NzViwYIFBp57EQeAZ0jsIXFlZiaamJly+fBlxcXGzfMpgVhgEttQSRWpqKrKysnDmzBkAwL1797B69WqkpKQg9rXvf7u7u8P2mjO5vYfL5vcWoWXDprA/r2mv9Df5B4HLysoCj6Wlpen+vZkmWSY6B4HDx5Rv5CiyGF0gRheI0QVidIEYXSBGF4jRBbLMhzNG4Qo0zQm8vQvE6AIxukCMLhCjC8ToAjG6QIwuEKMLxOgCMbpAjC4QowvEH61OgyvQAnEFmuYERheI0QVidIEYXSBGF4jRBWJ0gUwd3Yor0FYYBDZ1dCuuQPsHgbWLTdC6ugFgchD44AFTDAKb9mNYq65AA1MHgZVPTnIQWC89K9AAUFBQgPT0dGRmZmLjxo1oaWkx4rhBzDwIbMor3b8CXVFREfS111egAaCuri4wGtjV1QWXy4WRkRHYbLbZPHIQ/yCw7x+3oe7cMaNB4Egx5ZWudwUawJSVyMePH0NRFF3bcpFegQ7HIDBXoPHuFejy8nI4HA5s27YN586dg91u7A1syiBwaQmUDz/AxPET8GmaoefyM+VfVdY0DZmZmRgaGkJNTc2UFeiBgQF0dHRg06bg+cy2tjZUVFSgvb19ymzo/2MmM6ETp/4M3z/vwvaHE1Dsdvi++QYv3eVQv58PW+E23c8TqZlQU17poa5A++Xk5EBVVdy4cWOWTzzJPwhsO3QwaBBY+7wePp3/E4JIMuUbOUDfCvTY2BiGh4eRnJwM4NUbub6+PqSmps76ef3UzLVQvzwf/PjqNKhfcfo7ZJ2dncjKygr889OnT1FUVISxsTHY7XbExMSgvr4eSUlJBp7S/CwT/W0r0AkJCejo6DDwVNZkmehcgQ4fU76Ro8hidIEYXSBGF4jRBWJ0gRhdIMt8n24UDgLTnMDbu0CMLhCjC8ToAjG6QIwuEKMLxOgCMbpAjC4QowvE6AIxukCMLhCjC8ToAjG6QIwu0H8BN7lrlL7p9tMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 146.797x264.88 with 1 Axes>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit_nature.circuit.library import HartreeFock\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "init_state = HartreeFock(num_SO, [2, 2], qubit_converter)\n",
    "\n",
    "##############################\n",
    "\n",
    "init_state.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccbaee31",
   "metadata": {},
   "source": [
    "### 5. Ansatz\n",
    "\n",
    "One of the most important choices is the quantum circuit that you choose to approximate your ground state.\n",
    "Here is the example of qiskit circuit library that contains many possibilities for making your own circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f30d3bd3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/sympy/core/expr.py:3949: SymPyDeprecationWarning: \n",
      "\n",
      "expr_free_symbols method has been deprecated since SymPy 1.9. See\n",
      "https://github.com/sympy/sympy/issues/21494 for more info.\n",
      "\n",
      "  SymPyDeprecationWarning(feature=\"expr_free_symbols method\",\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAADWCAYAAAAaVxFlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxZ0lEQVR4nO3dd1hUZxr38e/MgCJgAVGJaBRFEFCJvQsaE8u6iSbKRo0by8YEMPa0tcUSzWts775R0nTdNVEj6EbXaDZGAQuaWKIRS7ATI8auoFhgeP+YiNJmKDPznIH7c11e0TNnzvzi9XjfnDPPeY4uOzs7GyGEEEJoil51ACGEEELkJw1aCCGE0CBp0EIIIYQGSYMWQgghNEgatBBCCKFB0qCFEEIIDZIGLYQQQmiQNGghhBBCg6RBCyGEEBokDVoIIYTQIGnQQgghhAZJgxZCCCE0SBq0EEIIoUHSoIUQQggNkgYthBBCaJA0aCGEEEKDpEELIYQQGiQNWgghhNAgadBCCCGEBjmpDiBEURw/ftziPh999BGjRo0yu0/jxo2tFUmIUpNxLcyRM2hRZixevFh1BCGsTsZ1+SUNWgghhNAgadBCCCGEBkmDFmVGbGys6ghCWJ2M6/JLGrQQQgihQdKgRZnRv39/1RGEsDoZ1+WX3GZlwS/bIO2Sms+uXBMCupXsvapylyazEEIbpO4Vj63qnjRoC9IuwY3zqlMUn6PmFkKo56j1w1FzF0YucYsyIyoqSnUEIaxOxnX5JQ1alBmWVlsSwhHJuC6/pEGLMqNLly6qIwhhdTKuyy9p0KLMuHz5suoIQlidjOvySyaJWcmE6DCOnduNweCMXm/A28OXQU9PIjRkgOpohXLEzEII7XDEGuJImaVBW9Hg7lMY3H0yWVmZrE/8iDkrB+Hn0xwfLz/V0QrliJkLExQUpDqCEFan9XHtiDXEUTLLJW4bMBic6NX2VbKMmZy6cFB1nCJxxMx5rV27VnUEIazOUca1I9YQrWeWBm0DDzLvszExGoA6Xv6K0xSNI2bOa+rUqaojCGF1jjKuHbGGaD2zXOK2opVb3ycmYR4Z99IwGJwZP+BzGtRuBsDsLwfRrfkg2gX1AWDa8r78uX0krQKeVRnZbObNPy7l+/0rcvZNvXaapr6deXfQl6rimhUTE8OMGTNUxyiW7GzQ6VSnEFqm9XEtdc92NH0GbTQamTdvHo0aNcLFxYWQkBASEhIICAhg5MiRquPlM+jpSXw98wax712hTePeHDoZl/NaxPOLWP6/KWTcS2fH4XW4uVRVPkjBfOZebUYwPyKe+RHxTBq8GpcKbgzr+b7CtGXD7zdhzY/wzhoYtxKmrYPNP0P6XdXJhCg+qXu2o+kGPWLECGbOnMlrr73G5s2bCQ8PZ+DAgZw+fZqWLVuqjleoyq4ejB/wOT8c/4bEpPUAeLjXpF+nMSxeP5qVW2fx+nMLFafMraDMDxmNRuasGsyIXnPw9qyvJmAZcewCfLgJdp+Auw9M225mwP8Om7ZfSVObT4iSkrpnfZpt0KtWrWL58uVs2LCBiRMn0rVrVyZNmkT79u3JzMykRYsWqiOaVcXVkxc7j2fZt3/HaDQC0KP1UM5fTqZvx9FUcfVUnDC/gjIDrNgyHV/vpnRs0ldduCJISEhQHcGs9LuwbDtkGSG7gNdv3YWl202XvYV4SOvj+nFS96xLsw169uzZ9OzZk9DQ0Fzb/fz8cHZ2plkz0/cFZ8+eJTQ0FH9/f5o2bcqOHTtUxC1Qv85juHYrlS37/52zrXZ1P81N5X9c3swHTmxlf/J3vPqnuYqTWXbkyBHVEczacwoeZBXcnMHUmFNvwGlZl0I8RuvjOi+pe9ajy87W3s/r58+fp27duixdupThw4fnem3gwIEcP36cn376CYAePXrw/PPPExkZSWJiIgMGDODMmTNUqFDB7GfoijgzZ97rcYQ0DCvR/0dB5q4eSu+2f6OJbyeL+x46Fc/Ej7uW6HOskfvarYtM/KQrs0dsLvIlntJkNmfcuHEW91m4cKHF/RYuVHeJ7cVJcfgEdEGnL/zn4uxsI/v++wGJaybZMZlQRavjWuqebeteUduuJmdxnz9vel6Yt7d3ru0ZGRkkJCTQq1cvAK5cucLOnTvZsGEDAB06dKB27drExcXRo0cP+4Yug774fia3797kw6+G5myrWyOAsf0/URfKgRmcKpptzgBkZ6M3mP/hUghhO1qqe5o8gz558iSNGjVi4cKFjB07Nmf79OnTee+991i8eDGRkZEcOHCA/v37c/r06Zx9wsPD6d69u9Vmee9bre75otXqQKuXSvZeVblLk9mc48ePW9wnMDCQY8eOmd2ncePG1opUbKv3wA+nCr/E/dDAdtC2oV0iCcW0Oq6l7hWPreqeJs+gGzRoQLNmzZg9ezaenp74+PgQGxvLpk2bADQ9g1uoM336dNURzOrYyPQ9tDkVnaB5PfvkEY5B6+Na2I4mJ4np9XpiYmIIDg4mIiKCYcOG4eXlRVRUFAaDIWeC2JNPPsnvv//OvXv3ct575swZ6tWTClcehYeHq45gVt3q0LmQxYoezogY0AYqaPLHZqGK1se1sB1NNmgAf39/4uLiuH37NikpKcycOZPDhw8TFBREpUqVAPDy8qJjx44sXboUgMTERH777Te6drX+JCWhfYGBgaojWNSvFfR5ClzzfM1c3R2Gd4FWvkpiCQ1zhHEtbMOhflbft28f7dq1y7Xt448/ZujQoSxatIgKFSqwatUqizO4hVBFr4PuwRDWGCauNm0b/Qz41pAlP4UQuTlMg05PTyc5OZnIyMhc2xs0aMD27dsVpTJZuuldjpzdRXD9jtSpEcDquDmMffFTQhqGsib+QxKPrKeWRz3e/MtyHmTe461Pu+NT3Y93Bn2h2dwAOw6vI3r9GFZO/pWMe+mayV0WOBke/b5BTXU5hCiNwmpIU9/OzF39Vy7dSMHJUIFJL6+mgpOLshpSWM7g+h0Yv6QLZy4e5uNxB/Hx8iP12hnmrv4rOnR4Va3D2wNXYNAbmLysD+kZN1gUtdNuuTV7iTsvd3d3srKyeOONN1RHyeXMxSRu373Fgsjt3Lpzlbv3bzMg9E1CGoZyPf0SB0/FsShqJ75PNGNX0tdUqujOpMGrVcc2m/uhHT/HUqNaXQDN5DYnLCxMdQQhrE6r49pcDTl14SBOThVYELmdHq2HsfXAl8pqiLmcBr0T04d+Teem/XP2d3epxqxhG1kQuR1vT19+PG6anDxr+Ea7Z3eYBq1VSWd20srftPh7i0bPoNc/OjVK/nUfIQ3C/nitO8fO7VYRsUDmcgP8cGwTLRp1R6dznCESHR2tOoIQVqfVcW2uhnhV9cFozAIgPeMGVVyrK8kI5nPqdDo8KtfKtX9lVw/cKlUFwMngjF6Xuzbak+NUX41Ku3ONf303jQnRYazc+j5pd67lvHb77g1cXaoA4OZSlfS7NxSlzM9cboAt+//F0y1eVpSuZCIiIlRHEMLqtDquzdWQKm5e3HuQwfAPA9m4O5pOTV/QZE5zrty8wP7kLTnNXQWH+Q5aqyq7evJKjxl0CH6OPUc3cvnmo7vk3VyqcvmPu+bv3L2Fu0s1RSnzM5f7p5PbCKrXHmcnx5psFx8frzqCEFan1XFtrobsT/6Oqm41WPbmMbb/HEtMwjyGPDNVczkLcz/zHh9+9QrjB3yGwaCuTcoZdCk18e3E4dOmSWqHTsXnXNYB8K/bmp9Pm55Ec+DE9wTWa1fgMVQwl/vsxSR2H9nAu5/15NzvR/jnt5NVxRRCaJS5GpKdnZ3z5Kqqbl7cvntTSUYwn7Mwi2JH8lyHKOrVCrJ1PLOkQZeSr3cTnAzOTIgOw8ngjEsFt5zXPNxr0rRBF8Yu7sSpCwfpENxXXdA8zOXu12k0H76+jTmvfku9WsEM6zlLYVIhhBaZqyGt/J8l5dIxJkSH8a//TeW59pFmjqQuJ8DMFeHsP/Edc1e/QmLSeo6e3c3OpHWs27GICdFh7Dz8H0XJ5RK3VYzoPSfn99t/jmV13Af4eDUipGEoL3V9m5e6vp3zesa9dD5Y9TIBdVuriJqLudwPPbylQEu5C2NpvWIhHJGWx7W5GjL1r7G59lVZQ8zlnDJkTb79N8xKy7dt8rI+eFZ5wqY589LkwzK0RBaNLx6VD8tYs2aNxWURVT4sI6+xX5r+u2iw2hxCHa2Oa6l7xVOuHpahJZUVLiJRms9WlVvl39e0adNk3WJR5qgY11L3tPG50qAtCOimOkHJOGpuIYR6jlo/HDV3YWSSmBBCCKFB0qBFmbFkyRLVEYSwOhnX5Zc0aFFmBAcHq44ghNXJuC6/pEGLMiM0NNTyTkI4GBnX5Zc0aCGEEEKDpEELIYQQGiQNWpQZrVtrd5UzIUpKxnX5JQ1alBl79+5VHUEIq5NxXX5JgxZCCCE0SBq0EEIIoUHSoEWZERsba3knIRyMjOvySxq0EEIIoUHSoEWZ0b9/f9URhLA6GdfllzzNyoJftkHaJTWfXblmyZ/Ooip3aTILIbRB6l7x2KruSYO2IO2SugeXl4aj5hZCqOeo9cNRcxdGLnGLMiMqKkp1BCGsTsZ1+SUNWpQZo0aNUh1BCKuTcV1+SYMWZUaXLl1UR7Do7gM4mAL//enRtuU74LvDcDwVsozqsgltcoRxLWxDvoO2kgnRYRw7txuDwRm93oC3hy+Dnp5EaMgA1dEK5YiZzbl8+bLqCIW6cQe+S4J9Z+B+Zu7XDqaYfgFUrQQdGkHXQKgg/zoF2h7XjlhDHCmzlAArGtx9CoO7TyYrK5P1iR8xZ+Ug/Hya4+PlpzpaoRwxs6P58TSs22c6e7bkZgZs/tnUyAe1B98ats8nRGk4Yg1xlMxyidsGDAYnerV9lSxjJqcuHFQdp0gcMXNeQUFBqiPks+kQrNxdtOb8uMtp8NH3cPhX2+QSjkOL47ogjlhDtJ5ZGrQNPMi8z8bEaADqePkrTlM0jpg5r7Vr16qOkMv2X0yXtUsqywjLd8JpRfejCm3Q2rgujCPWEK1nlkvcVrRy6/vEJMwj414aBoMz4wd8ToPazQCY/eUgujUfRLugPgBMW96XP7ePpFXAsyojm828+celfL9/Rc6+qddO09S3M+8O+lJVXLOmTp3KjBkzVMcA4PebsOGA+X0WDTb9d6yZv84sI6zcA2/1lu+kyystjeuCSN2zHU2fQRuNRubNm0ejRo1wcXEhJCSEhIQEAgICGDlypOp4+Qx6ehJfz7xB7HtXaNO4N4dOxuW8FvH8Ipb/bwoZ99LZcXgdbi5VlQ9SMJ+5V5sRzI+IZ35EPJMGr8alghvDer6vMK15MTExqiPkWLcfMq00I/tKGmw7Zp1jlcaDLNN34//9yXTp/vQlyM5Wnars09K4LojUPdvRdIMeMWIEM2fO5LXXXmPz5s2Eh4czcOBATp8+TcuWLVXHK1RlVw/GD/icH45/Q2LSegA83GvSr9MYFq8fzcqts3j9uYWKU+ZWUOaHjEYjc1YNZkSvOXh71lcT0IH8fgt+SbXuMRNPqL0F6+A5mLoWvkiErUdNl+7/sQXmbTb9ACGE1D3r02yDXrVqFcuXL2fDhg1MnDiRrl27MmnSJNq3b09mZiYtWrRQHdGsKq6evNh5PMu+/TtGo6my9mg9lPOXk+nbcTRVXD0VJ8yvoMwAK7ZMx9e7KR2b9FUXzoHsP2P9Y97KgOSL1j9uURw5b/ouPKOAiW4XbsD/2wJpd+0eS2iQ1D3r0myDnj17Nj179iQ0NDTXdj8/P5ydnWnWzPR9wdSpU/H390ev12vuuan9Oo/h2q1Utuz/d8622tX9NDeV/3F5Mx84sZX9yd/x6p/mKk5mWUJCguoIAJy76ljHNSc7G9b/BDozr9/MgJ3Jdo1VrmhlXBeV1D3r0eS0k/Pnz5OUlMS4cePyvZaSkkJwcDAVK1YEoGfPngwdOpThw4fbO2Yu8yPi821zc6nCuhnX7B+miCxlvnbrIh99PYrZIzbj7FTBzumK78iRI9SsWVN1DFJv2Oi4121zXHPOXYVLtyzvl3gCejWzfZ7ySCvjuiBS92xLsw0awNvbO9f2jIwMEhIS6NWrV862Dh06lOgzdLrCzglym/d6HCENw0r0GaWVkBBP64FdS/Rea+T+4vuZ3L57kw+/GpqzrW6NAMb2/6TQ95QmszkF/bCW18KFCy3ut3Ch7b8De/2zm1SsVCXnzw9naxemsNfzzu5e/823DA/tVfDONuLX+kX+NMbylam0u6DT62XWWDFpdVxL3bNt3csu4r8TTTZoLy8vAJKTk+ndu3fO9rlz55KamqrpCWKWvPXSctURimz0C4sZ/cJi1TEcjjHzvkMd15z7GTeLtN+De7elOYtCSd0rGV12UVu5HRmNRpo3b05qairz5s3Dx8eH2NhYNm3aREpKCnv27KFt27a53hMWFsaoUaPo37+/VbPsW63u+aLV6kCrl0r2XlW5S5PZnOPHj1vcJzAwkGPHzN+P1LhxY2tFKtSi/8HZK0XYrwj3QT+uayA8b+e5kZlZMO0/cPte4fvogLYN4aV2dotVZmh1XEvdKx5b1T1NThLT6/XExMQQHBxMREQEw4YNw8vLi6ioKAwGQ84EMSEeN336dNURAKhro4mqtjquOU4G0w8GhdEBej2E2v7nnnJLK+Na2J8mL3ED+Pv7ExcXl2vbkCFDCAoKolKlSopSCS0LDw9XHQGAZk/CDivPanY2QOPa1j1mUXULMj2Na2eyqSE/fsnNoIdhneGJamqylQdaGdfC/jR5Bl2Yffv25fv+ecqUKdSpU4fdu3fz2muvUadOHU6dOqUooVApMNDMqZ4d+dWEWlUs71ccLeuDq6IJpXod9G8N43tC6waPtvduBlP7QnAdNbnKC62Ma2F/mj2Dzis9PZ3k5GQiIyNzbZ85cyYzZ85UlMpk6aZ3OXJ2F8H1O1KnRgCr4+Yw9sVPCa7fgfFLunDm4mE+HncQHy8/Mu6l89an3fGp7sc7g77QZO6QhqE8P6UqfrWbAzDtlXVUcfVk8rI+pGfcYFHUTqW5tU6nM31X/Gm8dY5X0Ql6NLXOsUrjyeqmR2D+eNr052c1kEmoY65+fPT1G5xJPcwT1Rswrv9nGPQGpfWjODX65u0rTP3ncxgMzri5VGXyy19hNGYpqdsOcwbt7u5OVlYWb7zxhuoouZy5mMTtu7dYELmdW3eucvf+bQaEvklIw1AMeiemD/2azk0fTVyrVNGdSYNXK0xsYi43gK9305z1aB+u/jNr+EaVkR1KkA+0aWB5v6Lo2xI83KxzLCGswVz9+OXXvWRm3md+RDz1agXzw1FT3VBVP4pbo90rebAwcicLIhLw92nJnqMbldVth2nQWpV0Ziet/E2Lv7do9Ax6vSHnNZ1Oh0flWqqimWUuN0DKpWOMW9KZzze9U+R79lQLCwtTHSGXAW3A37vw18d+aXkGd7cgaNfQurmEY9HauAbz9SP16ml8nzBN5G1Y+ymOnEtUkvGh4tZog96AXm9qjVnZWfh4NbJf2DykQZdS2p1r/Ou7aUyIDmPl1vdJu6PdFXQeZyn38rdPsCBiO+l3rrP76H8VpSye6Oho1RFycTbAq2ElO5M26E2Xyf/8lOmSuSi/tDauwXz9qFMjgJ9Pm5YnPXhyG7czbihKaVKSGn085Uci/28rDp7cxhOevnZIWTBp0KVU2dWTV3rMYH5EPMN6zqKyBheDL4il3FVcPdHpdHRo0pezF5MUpSyeiIgI1RHycTaYvrd9NQxqFnHiWKNaMLGX6fYmac5Ci+PaXP3w83mK+t5NmPhxV+7cu0U1xVcRS1KjGz/ZhiVj9tGxST++3bvMDikLJg26lJr4duLw6e0AHDoVj9GYpThR0ZjLnXH/Nll//PnI2V3Uru4Y11jj4+NVRyhUsA+82weinoZOjaBedahUwdTA3SuaLoU/HQTv9IGo7nLbknhEi+PaUt0b8sxU5r0eRxXX6rQN/JOKiDmKW6MfPLZin5tLFSo4q7ut12FmcWuVr3cTnAzOTIgOI6hee1wquJFlzMx5feaKcJLO7uS3Kyf4S9hbdGjyvMK0j5jL/dvlE8yPGU6lCu54e/ry12dloQRr0OmgkbfplxCOzFz9MBqNvPlJN/R6A839nibwybYWjqYuK+Sv0Z5VnuDTb95Er9NTuZInbw9coSy7NGgrGNF7Ts7vt/8cy+q4D/DxakRIw1CmDFmTa9+Me+l8sOplAuq2tnfMfMzljh57IN/+k5f1wbPKE/aMKITQKHP1o6AnRqmsH8Wp0QALInI/4lNV3dbkWtxaImvSFo/KtbiLwh5rcZd1D2eeW3pKl7BMq+Na6l7x2KruyRm0BZUVPoa1NJ+tKrfKv681a9bIsoiizFExrqXuaeNz5QxaOAStPvWnPJIzaOuRcS3MkVncQgghhAZJgxZCCCE0SBq0KDOWLFmiOoIQVifjuvySBi3KjODgYNURhLA6GdfllzRoUWaEhoaqjiCE1cm4Lr+kQQshhBAaJA1alBmtW6tfnU0Ia5NxXX5JgxZlxt69e1VHEMLqZFyXX9KghRBCCA2SBi2EEEJokDRoUWbExsaqjiCE1cm4Lr+kQQshhBAaJA1alBn9+/dXHUEIq5NxXX7J4yYt+GUbpF1S89mVa0JAt5K9V1Xu0mQWQmiD1L3isVXdkwZtQdoldQ8uLw1HzS2EUM9R64ej5i6MXOIWZUZUVJTqCEJYnYzr8ksatCgzRo0apTqCEFYn47r8kgYtyowuXbqojlCmZRnhwvVHf76aDtnZ6vKUFzKuyy/5DlqUGZcvX1Ydocx5kAU/nYM9JyHlKmQaH702cz24VYQAb+jkD741QKdTl7WsknFdfkmDtpIJ0WEcO7cbg8EZvd6At4cvg56eRGjIANXRCuWImYX9HP0NvvoBbmYUvs/te3DgnOmXvze81BY83e2XUajliDXEkTLLJW4rGtx9Cv99P511712lR+uhzFk5iN+unFQdyyxHzFyYoKAg1RHKBGM2rNsHn8abb855JV+E//MNHClDs2i1QOvj2hFriKNklgZtAwaDE73avkqWMZNTFw6qjlMkjpg5r7Vr16qO4PCysyH2R9j+S8nefy8Tlm6HI79ZN1d55ijj2hFriNYzS4O2gQeZ99mYGA1AHS9/xWmKxhEz5zV16lTVERze3jOQaOFEYtFg06/CGLNhxS64cce62corRxnXjlhDtJ5Z099BG41GFixYwCeffMKvv/5KQEAA//jHPxg5ciShoaF8+umnqiPmsnLr+8QkzCPjXhoGgzPjB3xOg9rNAJj95SC6NR9Eu6A+AExb3pc/t4+kVcCzKiObzbz5x6V8v39Fzr6p107T1Lcz7w76UlVcs2JiYpgxY4bqGA7rVgb8Z591jnX3Aaz5AUZ2tc7xSuvuA3iQCa4VweBgpyVaH9dS92xH0w16xIgRrFu3jilTptCyZUsSExMZOHAgly9fZvz48arj5TPo6UkM7j6ZtDvXmR8zgkMn4+jVZgQAEc8v4t3PehDSMIx9yd/h5lJV+SAF85l7tRmR8/trty4y8ZOuDOv5vsq4woZ2nYCMB9Y73tEL8Nt18PGw3jGL69gF2HYUTvxu+rOLM7RtCN2CoGoldbnKEql7tqPZnyVXrVrF8uXL2bBhAxMnTqRr165MmjSJ9u3bk5mZSYsWLVRHLFRlVw/GD/icH45/Q2LSegA83GvSr9MYFq8fzcqts3j9uYWKU+ZWUOaHjEYjc1YNZkSvOXh71lcTUNhUlhF2n7D+cXfZ4JhFlXAcPomDk4+tzXz3gWn7gs2m+7iF9Ujdsz7NNujZs2fTs2dPQkNDc2338/PD2dmZZs2acf36dfr06YO/vz8hISE8++yznDypjZl4VVw9ebHzeJZ9+3eMRtPNoz1aD+X85WT6dhxNFVdPxQnzKygzwIot0/H1bkrHJn3VhSuChIQE1REc1u834dZd6x/3xEXrH7Mozl+D/+w3/b6gxVRuZcAXifbNVFKONK6l7lmXJhv0+fPnSUpKYsCA/PelpaSkEBwcTMWKFdHpdIwdO5bk5GQOHTpEnz59GDZsmILEBevXeQzXbqWyZf+/c7bVru6Hj5efwlTm5c184MRW9id/x6t/mqs4mWVHjhxRHcFh/XrNNse9nAYZ921zbHN2JoO5NVOygTOXTZfgtc7RxrXUPevRZWdrb7G+PXv20L59e7755ht69+6dsz0jI4OGDRvSq1cvli5dmu99+/bto2/fvpw/b/lGTF0Rlzya93ocIQ3Dipzdkrmrh9K77d9o4tvJ4r6HTsUz8eOSzbKxRu6H37/MHrG5yJd4SpPZnHHjxlncZ+HChRb3W7hQW5fYtKJN38m07z8z1zZzM7XNGZtnLs2/3wrk+oXjJUxWMkMXnqZqDV+L+yWsGMPB//3DDokKptVxLXXPtnWvqG1Xk5PEvLy8AEhOTs7VoOfOnUtqaiotW7Ys8H2LFi2ib9++9ohYLnzx/Uxu373Jh18NzdlWt0YAY/t/oi6UsBHbrdGps+GxC6PXG4q0n66I+4nyQ0t1T5Nn0EajkebNm5Oamsq8efPw8fEhNjaWTZs2kZKSwp49e2jbtm2u90yfPp3Nmzezbds2XF1drZZl32p1zxetVgdavVSy96rKXZrM5hw/bvkMLDAwkGPHjpndp3HjxtaKVKbsSoaYvUXb9+GZdd4z5cLMehHcXUqWq6T+uR1+Pm/5YR6juoNfLftkKohWx7XUveKxVd3T5HfQer2emJgYgoODiYiIYNiwYXh5eREVFYXBYKBZs2a59p81axYbN27k22+/tWpzFo5l+vTpqiM4rDo2mrtTzdX+zRmgo7/55qzTQc0q0LCm/TKVlIzr8kuTl7gB/P39iYuLy7VtyJAhBAUFUanSoxsYp0+fzqZNm9iyZQvVqlWzc0qhJeHh4aojOKzaHqZ7hO9a8T5oUNcAG9UyPWFrZ3L+13Q6cNLDyx0c4+lbMq7LL02eQRdm3759ub5/PnLkCO+99x5Xr14lLCyMp556iqeeekpdQKFUYGCg6ggOy9kArRtY/7jtFU3c1engxVbQr6XpLP5xjZ+AsT3gyepqshWXjOvyS7Nn0Hmlp6eTnJxMZGRkzrbg4OAiz4YTQpjXxR8ST5gWLbGGup5qLyHrdBDaGDr7w/hVpm3T+oKHm7pMQhSHwzRod3d3srKyVMco0NJN73Lk7C6C63ekTo0AVsfNYeyLn1LT40nmrv4rOnR4Va3D2wNXYNAbmLysD+kZN1gUtVOTuUMahrI/eQurt83BmG3ktT/Pp26NAN76tDs+1f14Z9AXSnML26hRBXo2hW8Olf5YBj0MbKeNS8j6x64TSnMuvcLqRnD9Doxf0oUzFw/z8biDOfc9r4n/kMQj66nlUY83/7KcB5n37F5LipsZYMfhdUSvH8PKyb+ScS9dSf1zqEvcWnTmYhK3795iQeR2bt25yt37txkQ+iYhDUNxd6nGrGEbWRC5HW9PX348vgmAWcM3Kk5tPve9Bxl8s+cTPhi5hfkR8fjXaUmliu5MGrxadWyzwsLCVEdweN2CTJeAzRn7peUZ3P1amr7XFqWnpXFtrm4Y9E5MH/o1nZv2z9n/evolDp6KY1HUTnyfaMaupK/tXkuKm/mhHT/HUqNaXQBl9U8adCklndlJK3/T4u8tGj2T6/7Lyq4euFWqCoCTwRm9Tjv3XJrLffTcbnQ6PX//vBcfrBpCxv3bqmIWS3R0tOoIDs+gh+FdIKh2yd6vw9ScO2nvyX0OS0vj2lzd0Ol0eFTOfc9a8q/7CGkQ9sf+3Tl2brfdsj5U3MwAPxzbRItG3dHp1LZIadCllHbnGv/6bhoTosNYufV90u7kXzPxys0L7E/ekjNItMBc7utpv3MtLZXZf9tMcL0OfLPbMRYmiYiIUB2hTKjgBH8Lhb4twKkYP1N6VYY3njF97yusR0vjuij17nG3797A1aUKAG4uVUm/e8MOKXMrbmaALfv/xdMtXrZDOvMc5jtorars6skrPWbQIfg59hzdyOWbue+Sv595jw+/eoXxAz7DYNDOX7e53G4uVWlSvxMGvYGn/LoRkzBPYdKii4+PVx2hzNDrISwQQp40TRzbfQrSC3mYRh1P6NQIWtQ3NXdhXVoa15bqXV5uLlW5/MfKIXfu3sLdpZodUuZW3Mw/ndxGUL32ODtVsFPCwskZdCk18e3E4dPbAdN6rEZj7olsi2JH8lyHKOrVClIRr1DmcgfUbU3KJdPKRacuHMTb0/KaxqJs8nCDPz0FM16ASX+GVzrBgNYQ3gZe62paJWxiL2jnJ825PLBU7/Lyr9uan0+bnsZ14MT3BNZrZ/OMeRU389mLSew+soF3P+vJud+P8M9vJ9sjZoHkn1Qp+Xo3wcngzIToMILqtcelghtZxkwAjp7dzc6kdfx+/RzrdiyiX6cxdGraT3FiE3O5q7nXoFmDUMYv6UJFZ1feHbxScVqhml5nmuVdo4rqJEIlc3UDYOaKcJLO7uS3Kyf4S9hbdGjyPE0bdGHs4k7UrPYkL3Qeq/nM/TqNpl+n0QCMXdyJYT1n2T3zQ9KgrWBE7zk5v9/+cyyr4z7Ax6sRIQ1D2TArLd/+k5f1wbOKhamydmAu94tdxvFil0dP0Mm4l84Hq14moG5rFVGLxNJ6xUI4Iq2Na3N1Y8qQNfn2f6nr27zU9e2cP6uoJcXN/NDDW2FV1T9NPixDS2TR+OJR+bCMNWvWWFwWUR6WUT49vC2spI/QtBWtjmupe8Vjq7onZ9AWVFa4ElJpPltVbpV/X9OmTZN1i0WZo2JcS93TxudKg7YgoJvqBCXjqLmFEOo5av1w1NyFkVncQgghhAZJgxZlxpIlS1RHEMLqZFyXX9KgRZkRHBysOoIQVifjuvySBi3KjNDQUNURhLA6GdfllzRoIYQQQoOkQQshhBAaJA1alBmtW2t3lTMhSkrGdfklDVqUGXv37lUdQQirk3FdfkmDFkIIITRIGrQQQgihQdKgRZkRGxurOoIQVifjuvySBi2EEEJokDRoUWb0799fdQQhrE7GdfklT7Oy4JdtkHZJzWdXrlnyp7Ooyl2azEIIbZC6Vzy2qnvSoC1Iu6TuweWl4ai5hRDqOWr9cNTchZFL3KLMiIqKUh1BCKuTcV1+SYMWZcaoUaNURxDC6mRcl1/SoEWZ0aVLF9URhIZcugVxx+DfOx9t++h7WLcPDpyF+5nKohWLjOvyS76DtpIJ0WEcO7cbg8EZvd6At4cvg56eRGjIANXRCuWImc25fPmy6ghCA85egc2H4JeL+V87+bvp1/ZfoFIFaO8HzzYBF2f75ywqLY9rR6whjpRZGrQVDe4+hcHdJ5OVlcn6xI+Ys3IQfj7N8fHyUx2tUI6YWYiCZBnhm0MQdxSyi7B/xn3YdhR+OgeD24NfLZtHLJMcsYY4Sma5xG0DBoMTvdq+SpYxk1MXDqqOUySOmDmvoKAg1RGEIllG+OcOU8MtSnN+3PXbEL0VDv9qk2il5ijj2hFriNYzS4O2gQeZ99mYGA1AHS9/xWmKxhEz57V27VrVEYQisXshqRS312Rlw/KdkHLVepmsxVHGtSPWEK1nlkvcVrRy6/vEJMwj414aBoMz4wd8ToPazQCY/eUgujUfRLugPgBMW96XP7ePpFXAsyojm828+celfL9/Rc6+qddO09S3M+8O+lJVXLOmTp3KjBkzVMcQdnb0N9h90vw+iwab/jvWzNDNMsLK3TChFzgbrJevtLQ+rqXu2Y6mz6CNRiPz5s2jUaNGuLi4EBISQkJCAgEBAYwcOVJ1vHwGPT2Jr2feIPa9K7Rp3JtDJ+NyXot4fhHL/zeFjHvp7Di8DjeXqsoHKZjP3KvNCOZHxDM/Ip5Jg1fjUsGNYT3fV5jWvJiYGNURhJ0Zs02zsq3l4k3YlWy941mD1se11D3b0XSDHjFiBDNnzuS1115j8+bNhIeHM3DgQE6fPk3Lli1VxytUZVcPxg/4nB+Of0Ni0noAPNxr0q/TGBavH83KrbN4/bmFilPmVlDmh4xGI3NWDWZErzl4e9ZXE1CIAvySClfSrXvMXSdMjV8Uj9Q969Nsg161ahXLly9nw4YNTJw4ka5duzJp0iTat29PZmYmLVq0UB3RrCqunrzYeTzLvv07RqMRgB6th3L+cjJ9O46miqun4oT5FZQZYMWW6fh6N6Vjk77qwglRgP1nrX/My2nwqwa/i3YEUvesS7MNevbs2fTs2ZPQ0NBc2/38/HB2dqZZM9P3BX379qVZs2Y0b96cNm3a8P3336uIW6B+ncdw7VYqW/b/O2db7ep+mpvK/7i8mQ+c2Mr+5O949U9zFSezLCEhQXUEYWe2mtT16zXbHLckHG1cS92zHk1OEjt//jxJSUmMGzcu32spKSkEBwdTsWJFAJYvX061atUA+OmnnwgLC+PatWsYDPad5TE/Ij7fNjeXKqyboaF/6XlYynzt1kU++noUs0dsxtmpgp3TFd+RI0eoWbOm6hjCTrKMcPmWbY6desM2xy0JLY9rqXu2pdkGDeDt7Z1re0ZGBgkJCfTq1Stn28PmDHDz5k10Oh3Z2Za/QNLpdEXKMu/1OEIahhVpX2tLSIin9cCuJXqvNXJ/8f1Mbt+9yYdfDc3ZVrdGAGP7f1Loe0qT2ZyCfljLa+HChRb3W7hQW9+BiZJzquhK1NLbubY9nK1dmMJezzu7+7Ol/yS87fBSpCsarY5rqXu2rXtF6VGg0Qbt5eUFQHJyMr17987ZPnfuXFJTU/NNEIuKimLz5s3cvHmTtWvX4uSkyf8tAN56abnqCEU2+oXFjH5hseoYQhTImHkfMBW7ov7AXVSZ9+9a9XjlndS9ktFlF7WV25HRaKR58+akpqYyb948fHx8iI2NZdOmTaSkpLBnzx7atm2b730JCQmMGzeO7du34+7ubpUs+1are75otTrQ6qWSvVdV7tJkNuf48eMW9wkMDOTYsWNm92ncuLG1IgkNeH+DaVKXJUW5D/pxL7SCLgElz1VUWh3XUveKx1Z1T5OTxPR6PTExMQQHBxMREcGwYcPw8vIiKioKg8GQM0Esr9DQUPR6Pbt27bJzYqEF06dPVx1B2FldG00KttVxS0LGdfml2WvB/v7+xMXF5do2ZMgQgoKCqFSpEgDp6elcvXqVevXqAaZJYqdOnSIwMNDueYV64eHhqiMIO2teDw6cs+4xPdygXnXrHrM0ZFyXX5pt0AXZt28f7dq1y/nz7du3+ctf/kJ6ejpOTk64uLjwxRdf8OSTTypMKVQpyqVAUbYE+YCHK1y/Y71jdmwEeg1dW5RxXX45TINOT08nOTmZyMjInG21atViz549ClOZLN30LkfO7iK4fkfq1Ahgddwcxr74KfW9g5n6z+cwGJxxc6nK5Je/wmjM4q1Pu+NT3Y93Bn2hydxuLlWJ3jAWgEvXz9Gv8xhe6DyWycv6kJ5xg0VRO5XmFuIhgx76tjQ9ycoavNyhsx2+e3ZkhdWN4PodGL+kC2cuHubjcQfx8fIjM+tBvm0Z99LtXgOLkxng+SlV8avdHIBpr6yjiqunkvqnoZ8TzXN3dycrK4s33nhDdZRczlxM4vbdWyyI3M6tO1e5e/82A0LfJKRhKO6VPFgYuZMFEQn4+7Rkz9GNVKrozqTBq1XHNpvbz+epnLVofZ9oRttA00L3s4ZvVJxaiPxCnoSW9c3vM/ZLyxPE9DoY2B4qOsxpi/2ZqxsGvRPTh35N56b9c/YvaJu9a2BxMwP4ejfNqYEPVz9TUf8cpkFrVdKZnbTyNy3+3qLRM+j1jxZIMegN6P+4VpaVnYWPVyMlGQtiLvdDGfdvcz3toqZXAHpcWFiY6ghCkZfagb+35f0KowMGtYeGGlwPREvj2lzd0Ol0eFSulWv/grbZW3EzA6RcOsa4JZ35fNM7Rb5n2RbkZ8VSSrtzjY27P2btjoWkZ9wgNCScau6P/pUfT/mRf/wnkgpOLgzoMkFh0tws5QbYe3wzrQJ6KkpYfNHR0aojCEWcDfBqGPxnHyRaePRkXpVdYGA70/fZWqSlcV2UuqE1Jcm8/O0TVK7kwf9d+zq7j/6XDsHP2SltbnIGXUqVXT15pccM5kfEM6znLCrnWQy+8ZNtWDJmHx2b9OPbvcsUpczPUm6AXUn/oVPTFxSkK5mIiAjVEYRCzgYIbwsR3Yo2C9vZAB384J0+2m3OoK1xXZS6oTUlyVzF1ROdTkeHJn05ezHJDikLJg26lJr4duLw6e0AHDoVj9GYlfPagz9WOgLTWq8VnCvZPV9hzOUGyMx6QMqlYzSsHaIiXonEx8erjiA0IOAJGNcTJvSEHk0hqDbUrALV3aGOB7RpAOFtYHo/U0N3q6g6sXlaGteW6oYWFTdzxv3bZP2xz5Gzu6hdvaHNMxZGLnGXkq93E5wMzkyIDiOoXntcKriRZcwE4NSFg3z6zZvodXoqV/Lk7YErFKd9xFxugJ9ObuOpht0UJhSidOpWN/0S1mOpbsxcEU7S2Z38duUEfwl7iw5Nni9wm5Yz1/Sox/yY4VSq4I63py9/fVbdQjHSoK1gRO85Ob/f/nMsq+M+wMerESENQ1kQkftRcRn30vlg1csE1G1t75j5mMvdOqAHrQN65Np/8rI+eFZ5wt4xhRAaYq5uTBmyJt/+ebepqIHFzRw99kC+bSrqnybX4tYSWZO2eFSuxV0Usha30BKtjmupe8Vjq7onZ9AWVFY4QbE0n60qt8q/rzVr1siyiKLMUTGupe5p43PlDFo4BK0+9UeI0pBxLcyRWdxCCCGEBkmDFkIIITRIGrQoM5YsWaI6ghBWJ+O6/JIGLcqM4OBg1RGEsDoZ1+WXNGhRZoSGhqqOIITVybguv6RBCyGEEBok90ELh1CU20imTZsmt5sIhyLjWpgj90ELIYQQGiSXuIUQQggNkgYthBBCaJA0aCGEEEKDpEELIYQQGiQNWgghhNAgadBCCCGEBkmDFkIIITRIGrQQQgihQdKghRBCCA2SBi2EEEJo0P8Hmq9lyFc0ZKoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 628.397x264.88 with 1 Axes>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.circuit.library import EfficientSU2, TwoLocal, NLocal, PauliTwoDesign\n",
    "from qiskit_nature.circuit.library import UCCSD, PUCCD, SUCCD\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "ansatz = EfficientSU2(num_qubits=qubit_op.num_qubits, reps=1, entanglement='linear', insert_barriers=True)\n",
    "\n",
    "##############################\n",
    "\n",
    "ansatz.decompose().draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5c8f3f2",
   "metadata": {},
   "source": [
    "## Ground state energy calculation \n",
    "\n",
    "### Calculation using NumPy\n",
    "\n",
    "For learning purposes, we can solve the problem exactly with the exact diagonalization of the Hamiltonian matrix so we know where to aim with VQE. Of course, the dimensions of this matrix scale exponentially in the number of molecular orbitals so you can try doing this for a large molecule of your choice and see how slow this becomes. For very large systems you would run out of memory trying to store their wavefunctions.\n",
    "\n",
    "<center><img src=\"resources/vqe.png\" width=\"600\"></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "38cc3c99",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact electronic energy: -0.766115 Hartree\n",
      "\n",
      "=== GROUND STATE ENERGY ===\n",
      " \n",
      "* Electronic ground state energy (Hartree): -4043.632437386865\n",
      "  - computed part:      -0.766114826773\n",
      "  - ActiveSpaceTransformer extracted energy part: -4042.866322560092\n",
      " \n",
      "=== MEASURED OBSERVABLES ===\n",
      " \n",
      "  0:  # Particles: 2.000 S: 0.000 S^2: 0.000 M: 0.000\n"
     ]
    }
   ],
   "source": [
    "from qiskit.algorithms import NumPyMinimumEigensolver\n",
    "from qiskit_nature.algorithms import GroundStateEigensolver\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "numpy_solver = NumPyMinimumEigensolver()\n",
    "numpy_ground_state_solver = GroundStateEigensolver(qubit_converter, numpy_solver)\n",
    "numpy_results = numpy_ground_state_solver.solve(es_problem)\n",
    "\n",
    "##############################\n",
    "\n",
    "exact_energy = numpy_results.computed_energies[0]\n",
    "print(f\"Exact electronic energy: {exact_energy:.6f} Hartree\\n\")\n",
    "print(numpy_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "196be28b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitting your answer for 2c. Please wait...\n",
      "Congratulations 🎉! Your answer is correct and has been submitted.\n"
     ]
    }
   ],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2c\n",
    "grade_ex2c(numpy_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1656df6",
   "metadata": {},
   "source": [
    "### Calculation using VQE\n",
    "\n",
    "The next step would be to use VQE to calculate this ground state energy and you would have found the solution to one half of your electronic problem!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "092e0e70",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit.providers.aer import StatevectorSimulator, QasmSimulator\n",
    "from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA, SLSQP\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "backend = StatevectorSimulator()\n",
    "optimizer = SLSQP()\n",
    " \n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bd90631f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit.algorithms import VQE\n",
    "from qiskit_nature.algorithms import VQEUCCFactory, GroundStateEigensolver\n",
    "from jupyterplot import ProgressPlot\n",
    "import numpy as np\n",
    "\n",
    "error_threshold = 10 # mHartree\n",
    "\n",
    "np.random.seed(5)  # fix seed for reproducibility\n",
    "initial_point = np.random.random(ansatz.num_parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f90f9121",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# for live plotting\n",
    "pp = ProgressPlot(plot_names=['Energy'],\n",
    "                  line_names=['Runtime VQE', f'Target + {error_threshold}mH', 'Target']) \n",
    "\n",
    "intermediate_info = {\n",
    "    'nfev': [],\n",
    "    'parameters': [],\n",
    "    'energy': [],\n",
    "    'stddev': []\n",
    "}\n",
    "\n",
    "def callback(nfev, parameters, energy, stddev):\n",
    "    intermediate_info['nfev'].append(nfev)\n",
    "    intermediate_info['parameters'].append(parameters)\n",
    "    intermediate_info['energy'].append(energy)\n",
    "    intermediate_info['stddev'].append(stddev)\n",
    "    pp.update([[energy, exact_energy+error_threshold/1000, exact_energy]])\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "vqe = VQEUCCFactory(backend)\n",
    "vqe_ground_state_solver = GroundStateEigensolver(qubit_converter, vqe)\n",
    "vqe_results = vqe_ground_state_solver.solve(es_problem)\n",
    "\n",
    "##############################\n",
    "\n",
    "print(vqe_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33c9fc09",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "error = (vqe_results.computed_energies[0] - exact_energy) * 1000 # mHartree\n",
    "print(f'Error is: {error:.3f} mHartree')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50843931",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2d\n",
    "grade_ex2d(vqe_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144bf8c7",
   "metadata": {},
   "source": [
    "## Excited state calculation \n",
    "\n",
    "### Calculation using QEOM\n",
    "\n",
    "For the molecule of our interest we also need to compute the same but this time for the excited state of our molecular hamiltonian. Since we've already defined the system, we would now need to access the excitation energy using the quantum Equation of Motion (qEOM) algorithm which does this by solving the following pseudo-eigenvalue problem\n",
    "\n",
    "<center><img src=\"resources/math-1.svg\" width=\"400\"></center>\n",
    "\n",
    "with\n",
    "\n",
    "<center><img src=\"resources/math-2.svg\" width=\"300\"></center>\n",
    "\n",
    "where each corresponding matrix element must be measured on our quantum computer with its corresponding ground state.\n",
    "\n",
    "To learn more, you can read up about excited state calculation with [**this tutorial**](https://qiskit.org/documentation/nature/tutorials/04_excited_states_solvers.html), and about qEOM itself from the [**corresponding paper by Ollitrault et al., 2019**](https://arxiv.org/abs/1910.12890)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10089c13",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit_nature.algorithms import QEOM\n",
    "\n",
    "##############################\n",
    "# Provide your code here  \n",
    "\n",
    "qeom_excited_state_solver = \n",
    "qeom_results = \n",
    "\n",
    "##############################\n",
    "\n",
    "print(qeom_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "966911a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2e\n",
    "grade_ex2e(qeom_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeb75226",
   "metadata": {},
   "source": [
    "Finally, you just need to calculate the band gap or energy gap (which is the minimum amount of energy required by an electron to break free of its ground state into its excited state) by computing the difference of the two sets of energies that you have calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd6c33cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bandgap = qeom_results.computed_energies[1] - qeom_results.computed_energies[0]\n",
    "bandgap # in Hartree"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e635bb",
   "metadata": {},
   "source": [
    "## Running VQE on the cloud using Qiskit Runtime\n",
    "\n",
    "Qiskit Runtime is a new architecture offered by IBM Quantum that streamlines computations requiring many iterations. These experiments will execute significantly faster within this improved hybrid quantum/classical process.\n",
    "\n",
    "Qiskit Runtime allows authorized users to upload their Qiskit quantum programs for themselves or others to use. A Qiskit quantum program, also called a Qiskit Runtime program, is a piece of Python code that takes certain inputs, performs quantum and maybe classical computation, interactively provides intermediate results if desired, and returns the processing results. The same or other authorized users can then invoke these quantum programs by simply passing in the required input parameters.\n",
    "\n",
    "\n",
    "<center><img src=\"resources/qiskit-runtime1.gif\" width=\"600\"></center>\n",
    "\n",
    "<center><img src=\"resources/runtime_arch.png\" width=\"600\"></center>\n",
    "\n",
    "\n",
    "To run the VQE using Qiskit Runtime, we only have to do very few changes from the local VQE run and mainly have to replace the VQE class by the VQEProgram class. Both follow the same MinimumEigensolver interface and thus share the compute_minimum_eigenvalue method to execute the algorithm and return the same type of result object. Merely the signature of the initializer differs sligthly.\n",
    "\n",
    "We start by choosing the provider with access to the Qiskit Runtime service and the backend to execute the circuits on.\n",
    "\n",
    "For more information about Qiskit Runtime, please refer to [**VQEProgram**](https://qiskit.org/documentation/partners/qiskit_runtime/tutorials/vqe.html#Runtime-VQE:-VQEProgram) and [**Leveraging Qiskit Runtime**](https://qiskit.org/documentation/nature/tutorials/07_leveraging_qiskit_runtime.html) tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a34f3c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qc_grader.util import get_challenge_provider\n",
    "provider = get_challenge_provider()\n",
    "if provider:\n",
    "    backend = provider.get_backend('ibmq_qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3835b7fc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit_nature.runtime import VQEProgram\n",
    "\n",
    "error_threshold = 10 # mHartree\n",
    "\n",
    "# for live plotting\n",
    "pp = ProgressPlot(plot_names=['Energy'],\n",
    "                  line_names=['Runtime VQE', f'Target + {error_threshold}mH', 'Target']) \n",
    "\n",
    "intermediate_info = {\n",
    "    'nfev': [],\n",
    "    'parameters': [],\n",
    "    'energy': [],\n",
    "    'stddev': []\n",
    "}\n",
    "\n",
    "def callback(nfev, parameters, energy, stddev):\n",
    "    intermediate_info['nfev'].append(nfev)\n",
    "    intermediate_info['parameters'].append(parameters)\n",
    "    intermediate_info['energy'].append(energy)\n",
    "    intermediate_info['stddev'].append(stddev)\n",
    "    pp.update([[energy,exact_energy+error_threshold/1000, exact_energy]])\n",
    "\n",
    "##############################\n",
    "# Provide your code here \n",
    "\n",
    "optimizer = {\n",
    "    'name': 'QN-SPSA',  # leverage the Quantum Natural SPSA\n",
    "    # 'name': 'SPSA',  # set to ordinary SPSA\n",
    "    'maxiter': 100,\n",
    "}\n",
    "    \n",
    "runtime_vqe = \n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71bb7a75",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 2f grading**\n",
    "    \n",
    "The grading for this exercise is slightly different from the previous exercises. \n",
    "    \n",
    "1. You will first need to use `prepare_ex2f` to submit a runtime job to IBM Quantum (to run on a simulator), using `runtime_vqe (VQEProgram)`, `qubit_converter (QubitConverter)`, `es_problem (ElectronicStructureProblem)`. Depending on the queue, the job can take up to a few minutes to complete. Under the hood, the `prepare_ex2f` does the following:\n",
    "```python\n",
    "runtime_vqe_groundstate_solver = GroundStateEigensolver(qubit_converter, runtime_vqe)\n",
    "runtime_vqe_result = runtime_vqe_groundstate_solver.solve(es_problem)\n",
    "```\n",
    "    \n",
    "2. After the job has completed, you can use `grade_ex2f` to check the answer and submit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "348bc7bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submit a runtime job using the following code\n",
    "from qc_grader import prepare_ex2f\n",
    "runtime_job = prepare_ex2f(runtime_vqe, qubit_converter, es_problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b49aa7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex2f\n",
    "grade_ex2f(runtime_job)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75774d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(runtime_job.result().get(\"eigenvalue\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe571521",
   "metadata": {},
   "source": [
    "Congratulations! You have submitted your first Qiskit Runtime program and passed the exercise.\n",
    "\n",
    "But the fun is not over! We have reserved a dedicated quantum system for the quantum challenge. As bonus exercise (not graded), you can try your hands on submitting a VQE runtime job to a real quantum system! \n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Running VQE on a real quantum system (Optional)**\n",
    "    \n",
    "We have reserved a dedicated quantum system [`ibm_perth`](https://quantum-computing.ibm.com/services?services=systems&system=ibm_perth) for this challenge. Please follow the steps below to submit runtime job on the real quantum system.\n",
    "    \n",
    "1. Update backend selection to `ibm_perth` and pass it to `runtime_vqe` again\n",
    "    ```python\n",
    "    backend = provider.get_backend('ibm_perth')\n",
    "    runtime_vqe = VQEProgram(...\n",
    "                             backend=backend, \n",
    "                             ...)\n",
    "    ```\n",
    "2. Set `real_device` flag in `prepare_ex2f` to `True`.\n",
    "3. Run `prepare_ex2f` to submit a runtime job to `ibm_perth`.\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "Note: Qiskit runtime speeds up VQE by up to 5 times. However, each runtime job can still take 30 ~ 60 minutes of quantum processor time. Therefore, **the queue time for completing a job can be hours or even days**, depending on how many participants are submitting jobs. \n",
    "\n",
    "To ensure a pleasant experience for all participants, please only submit a job to the real quantum system after trying with these settings using the simulator:\n",
    "1. Consider using `PartiyMapper` and set `two_qubit_reduction=True` to reduce number of qubits to 2 and make the VQE program converge to ground state energy faster (with lower number of iterations).\n",
    "1. Limit optimizer option `maxiter=100` or less. Use the simulator runs to find an optimal low number of iterations.\n",
    "1. Verify your runtime program is correct by passing `grade_ex2f` with simulator as backend.\n",
    "1. Limit your jobs to only 1 job per participant to allow more participants to try runtime on a real quantum system.\n",
    "    \n",
    "Don't worry if your job is getting too long to execute or it can't be executed before the challenge ends. This is an optional exercise. You can still pass all challenge exercises and get a digital badge without running a job on the real quantum system.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7b579e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please change backend to ibm_perth before running the following code\n",
    "runtime_job_real_device = prepare_ex2f(runtime_vqe, qubit_converter, es_problem, real_device=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b914ba92",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(runtime_job_real_device.result().get(\"eigenvalue\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d50b0c",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "**Created by:** Junye Huang, Samanvay Sharma\n",
    "\n",
    "**Version:** 1.0.0"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Qiskit v0.31.0 (ipykernel)",
   "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"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}