content/ex5/ex5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 5 - Variational quantum eigensolver\n",
    "\n",
    "\n",
    "## Historical background\n",
    "\n",
    "During the last decade, quantum computers matured quickly and began to realize Feynman's initial dream of a computing system that could simulate the laws of nature in a quantum way. A 2014 paper first authored by Alberto Peruzzo introduced the **Variational Quantum Eigensolver (VQE)**, an algorithm meant for finding the ground state energy (lowest energy) of a molecule, with much shallower circuits than other approaches.[1] And, in 2017, the IBM Quantum team used the VQE algorithm to simulate the ground state energy of the lithium hydride molecule.[2]\n",
    "\n",
    "VQE's magic comes from outsourcing some of the problem's processing workload to a classical computer. The algorithm starts with a parameterized quantum circuit called an ansatz (a best guess) then finds the optimal parameters for this circuit using a classical optimizer. The VQE's advantage over classical algorithms comes from the fact that a quantum processing unit can represent and store the problem's exact wavefunction, an exponentially hard problem for a classical computer. \n",
    "\n",
    "This exercise 5 allows you to realize Feynman's dream yourself, setting up a variational quantum eigensolver to determine the ground state and the energy of a molecule. This is interesting because the ground state can be used to calculate various molecular properties, for instance the exact forces on nuclei than can serve to run molecular dynamics simulations to explore what happens in chemical systems with time.[3]\n",
    "\n",
    "\n",
    "### References\n",
    "\n",
    "1. Peruzzo, Alberto, et al. \"A variational eigenvalue solver on a photonic quantum processor.\" Nature communications 5.1 (2014): 1-7.\n",
    "2. Kandala, Abhinav, et al. \"Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.\" Nature 549.7671 (2017): 242-246.\n",
    "3. Sokolov, Igor O., et al. \"Microcanonical and finite-temperature ab initio molecular dynamics simulations on quantum computers.\" Physical Review Research 3.1 (2021): 013125.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "For the implementation of VQE, you will be able to make choices on how you want to compose your simulation, in particular focusing on the ansatz quantum circuits.\n",
    "This is motivated by the fact that one of the important tasks when running VQE on noisy quantum computers is to reduce the loss of fidelity (which introduces errors) by finding the most compact quantum circuit capable of representing the ground state.\n",
    "Practically, this entails to minimizing the number of two-qubit gates (e.g. CNOTs) while not loosing accuracy.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "<b>Goal</b> \n",
    "\n",
    "Find the shortest ansatz circuits for representing accurately the ground state of given problems. Be creative!\n",
    "    \n",
    "<b>Plan</b> \n",
    "    \n",
    "First you will learn how to compose a VQE simulation for the smallest molecule and then apply what you have learned to a case of a larger one.\n",
    "    \n",
    "**1. Tutorial -  VQE for H$_2$:** familiarize yourself with VQE and select the best combination of ansatz/classical optimizer by running statevector simulations.\n",
    "\n",
    "**2. Final Challenge - VQE for LiH:** perform similar investigation as in the first part but restricting to statevector simulator only. Use the qubit number reduction schemes available in Qiskit and find the optimal circuit for this larger system. Optimize the circuit and use your imagination to find ways to select the best building blocks of parameterized circuits and compose them to construct the most compact ansatz circuit for the ground state, better than the ones already available in Qiskit.  \n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "Below is an introduction to the theory behind VQE simulations. You don't have to understand the whole thing before moving on. Don't be scared!\n",
    "\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "Here below is the general workflow representing how the molecular simulations using VQE are performed on quantum computers.\n",
    "\n",
    "<img src=\"resources/workflow.png\" width=800 height= 1400/>\n",
    "\n",
    "The core idea hybrid quantum-classical approach is to outsource to **CPU (classical processing unit)** and **QPU (quantum processing unit)** the parts that they can do best. The CPU takes care of listing the terms that need to be measured to compute the energy and also optimizing the circuit parameters. The QPU implements a quantum circuit representing the quantum state of a system and measures the energy. Some more details are given below:\n",
    "\n",
    "**CPU** can compute efficiently the energies associated to electron hopping and interactions (one-/two-body integrals by means of a Hartree-Fock calculation) that serve to represent the total energy operator, Hamiltonian. The [Hartree–Fock (HF) method](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method#:~:text=In%20computational%20physics%20and%20chemistry,system%20in%20a%20stationary%20state.) efficiently computes an approximate grounds state wavefunction by assuming that the latter can be represented by a single Slater determinant (e.g. for H$_2$ molecule in STO-3G basis with 4 spin-orbitals and qubits, $|\\Psi_{HF} \\rangle = |0101 \\rangle$ where electrons occupy the lowest energy spin-orbitals). What QPU does later in VQE is finding a quantum state (corresponding circuit and its parameters) that can also represent other states associated missing electronic correlations (i.e. $\\sum_i c_i |i\\rangle$ states in  $|\\Psi \\rangle = c_{HF}|\\Psi_{HF} \\rangle + \\sum_i c_i |i\\rangle $ where $i$ is a bitstring). \n",
    "\n",
    "After a HF calculation, operators in the Hamiltonian are mapped to measurements on a QPU using fermion-to-qubit transformations (see Hamiltonian section below). One can further analyze the properties of the system to reduce the number of qubits or shorten the ansatz circuit:\n",
    "\n",
    "- For Z2 symmetries and two-qubit reduction, see [Bravyi *et al*, 2017](https://arxiv.org/abs/1701.08213v1).\n",
    "- For entanglement forging, see [Eddins *et al.*, 2021](https://arxiv.org/abs/2104.10220v1).\n",
    "- For the adaptive ansatz see, [Grimsley *et al.*,2018](https://arxiv.org/abs/1812.11173v2), [Rattew *et al.*,2019](https://arxiv.org/abs/1910.09694), [Tang *et al.*,2019](https://arxiv.org/abs/1911.10205). You may use the ideas found in those works to find ways to shorten the quantum circuits.\n",
    "\n",
    "**QPU** implements quantum circuits (see Ansatzes section below), parameterized by angles $\\vec\\theta$, that would represent the ground state wavefunction by placing various single qubit rotations and entanglers (e.g. two-qubit gates). The quantum advantage lies in the fact that QPU can efficiently represent and store the exact wavefunction, which becomes intractable on a classical computer for systems that have more than a few atoms.  Finally, QPU measures the operators of choice (e.g. ones representing a Hamiltonian).\n",
    "\n",
    "Below we go slightly more in mathematical details of each component of the VQE algorithm. It might be also helpful if you watch our [video episode about VQE](https://www.youtube.com/watch?v=Z-A6G0WVI9w).\n",
    "\n",
    "\n",
    "### Hamiltonian \n",
    "\n",
    "Here we explain how we obtain the operators that we need to measure to obtain the energy of a given system.\n",
    "These terms are included in the molecular Hamiltonian defined as:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{H} &=\\sum_{r s} h_{r s} \\hat{a}_{r}^{\\dagger} \\hat{a}_{s} \\\\\n",
    "&+\\frac{1}{2} \\sum_{p q r s} g_{p q r s} \\hat{a}_{p}^{\\dagger} \\hat{a}_{q}^{\\dagger} \\hat{a}_{r} \\hat{a}_{s}+E_{N N}\n",
    "\\end{aligned}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "h_{p q}=\\int \\phi_{p}^{*}(r)\\left(-\\frac{1}{2} \\nabla^{2}-\\sum_{I} \\frac{Z_{I}}{R_{I}-r}\\right) \\phi_{q}(r)\n",
    "$$\n",
    "$$\n",
    "g_{p q r s}=\\int \\frac{\\phi_{p}^{*}\\left(r_{1}\\right) \\phi_{q}^{*}\\left(r_{2}\\right) \\phi_{r}\\left(r_{2}\\right) \\phi_{s}\\left(r_{1}\\right)}{\\left|r_{1}-r_{2}\\right|} \n",
    "$$\n",
    "\n",
    "where the $h_{r s}$  and $g_{p q r s}$ are the one-/two-body integrals (using the Hartree-Fock method) and $E_{N N}$ the nuclear repulsion energy. \n",
    "The one-body integrals represent the kinetic energy of the electrons and their interaction with nuclei. \n",
    "The two-body integrals represent the electron-electron interaction.\n",
    "The $\\hat{a}_{r}^{\\dagger}, \\hat{a}_{r}$ operators represent creation and annihilation of electron in spin-orbital $r$ and require mappings to operators, so that we can measure them on a quantum computer.\n",
    "Note that VQE minimizes the electronic energy so you have to retrieve and add the nuclear repulsion energy $E_{NN}$ to compute the total energy. \n",
    "    \n",
    "\n",
    "\n",
    "So, for every non-zero matrix element in the $ h_{r s}$ and $g_{p q r s}$ tensors, we can construct corresponding Pauli string (tensor product of Pauli operators) with the following fermion-to-qubit transformation. \n",
    "For instance, in Jordan-Wigner mapping for an orbital $r = 3$, we obtain the following Pauli string:\n",
    "$$\n",
    "\\hat a_{3}^{\\dagger}= \\hat \\sigma_z \\otimes \\hat \\sigma_z \\otimes\\left(\\frac{ \\hat \\sigma_x-i \\hat \\sigma_y}{2}\\right) \\otimes 1 \\otimes \\cdots \\otimes 1\n",
    "$$\n",
    "where $\\hat \\sigma_x, \\hat \\sigma_y, \\hat \\sigma_z$ are the well-known Pauli operators. The tensor products of $\\hat \\sigma_z$ operators are placed to enforce the fermionic anti-commutation relations.\n",
    "A representation of the Jordan-Wigner mapping between the 14 spin-orbitals of a water molecule and some 14 qubits is given below:\n",
    "\n",
    "<img src=\"resources/mapping.png\" width=600 height= 1200/>\n",
    "\n",
    "\n",
    "Then, one simply replaces the one-/two-body excitations (e.g. $\\hat{a}_{r}^{\\dagger} \\hat{a}_{s}$, $\\hat{a}_{p}^{\\dagger} \\hat{a}_{q}^{\\dagger} \\hat{a}_{r} \\hat{a}_{s}$) in the Hamiltonian by corresponding Pauli strings (i.e. $\\hat{P}_i$, see picture above). The resulting operator set is ready to be measured on the QPU.\n",
    "For additional details see [Seeley *et al.*, 2012](https://arxiv.org/abs/1208.5986v1).\n",
    "\n",
    "### Ansatzes\n",
    "\n",
    "There are mainly 2 types of ansatzes you can use for chemical problems. \n",
    "\n",
    "- **q-UCC ansatzes** are physically inspired, and roughly map the electron excitations to quantum circuits. The q-UCCSD ansatz (`UCCSD`in Qiskit) possess all possible single and double electron excitations. The paired double q-pUCCD (`PUCCD`) and singlet q-UCCD0 (`SUCCD`) just consider a subset of such excitations (meaning significantly shorter circuits) and have proved to provide good results for dissociation profiles. For instance, q-pUCCD doesn't have single excitations and the double excitations are paired as in the image below.\n",
    "- **Heuristic ansatzes (`TwoLocal`)** were invented to shorten the circuit depth but still be able to represent the ground state. \n",
    "As in the figure below, the R gates represent the parametrized single qubit rotations and $U_{CNOT}$ the entanglers (two-qubit gates). The idea is that after repeating certain $D$-times the same block (with independent parameters) one can reach the ground state. \n",
    "\n",
    "For additional details refer to [Sokolov *et al.* (q-UCC ansatzes)](https://arxiv.org/abs/1911.10864v2) and [Barkoutsos *et al.* (Heuristic ansatzes)](https://arxiv.org/pdf/1805.04340.pdf).\n",
    "\n",
    "<img src=\"resources/ansatz.png\" width=700 height= 1200/>\n",
    "\n",
    "\n",
    "\n",
    "### VQE\n",
    "\n",
    "Given a Hermitian operator $\\hat H$ with an unknown minimum eigenvalue $E_{min}$, associated with the eigenstate $|\\psi_{min}\\rangle$, VQE provides an estimate $E_{\\theta}$, bounded by $E_{min}$:\n",
    "\n",
    "\\begin{align*}\n",
    "    E_{min} \\le E_{\\theta} \\equiv \\langle \\psi(\\theta) |\\hat H|\\psi(\\theta) \\rangle\n",
    "\\end{align*}  \n",
    "\n",
    "where $|\\psi(\\theta)\\rangle$ is the trial state associated with $E_{\\theta}$. By applying a parameterized circuit, represented by $U(\\theta)$, to some arbitrary starting state $|\\psi\\rangle$, the algorithm obtains an estimate $U(\\theta)|\\psi\\rangle \\equiv |\\psi(\\theta)\\rangle$ on $|\\psi_{min}\\rangle$. The estimate is iteratively optimized by a classical optimizer by changing the parameter $\\theta$ and minimizing the expectation value of $\\langle \\psi(\\theta) |\\hat H|\\psi(\\theta) \\rangle$. \n",
    "\n",
    "As applications of VQE, there are possibilities in molecular dynamics simulations, see [Sokolov *et al.*, 2021](https://arxiv.org/abs/2008.08144v1), and excited states calculations, see [Ollitrault *et al.*, 2019](https://arxiv.org/abs/1910.12890) to name a few.\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    " \n",
    "<b> References for additional details</b> \n",
    "\n",
    "For the qiskit-nature tutorial that implements this algorithm see [here](https://qiskit.org/documentation/nature/tutorials/01_electronic_structure.html)\n",
    "but this won't be sufficient and you might want to look on the [first page of github repository](https://github.com/Qiskit/qiskit-nature) and the [test folder](https://github.com/Qiskit/qiskit-nature/tree/main/test) containing tests that are written for each component, they provide the base code for the use of each functionality.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Tutorial - VQE for H$_2$ molecule \n",
    "\n",
    "\n",
    "\n",
    "In this part, you will simulate H$_2$ molecule using the STO-3G basis with the PySCF driver and Jordan-Wigner mapping.\n",
    "We will guide you through the following parts so then you can tackle harder problems.\n",
    "    \n",
    "\n",
    "\n",
    "#### 1. Driver\n",
    "\n",
    "The interfaces to the classical chemistry codes that are available in Qiskit are called drivers.\n",
    "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers import PySCFDriver\n",
    "\n",
    "molecule = \"H .0 .0 .0; H .0 .0 0.739\"\n",
    "driver = PySCFDriver(atom=molecule)\n",
    "qmolecule = driver.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "<b> Tutorial questions 1</b> \n",
    "    \n",
    "Look into the attributes of `qmolecule` and answer the questions below.\n",
    "\n",
    "    \n",
    "1. We need to know the basic characteristics of our molecule. What is the total number of electrons in your system?\n",
    "2. What is the number of molecular orbitals?\n",
    "3. What is the number of spin-orbitals?\n",
    "3. How many qubits would you need to simulate this molecule with Jordan-Wigner mapping?\n",
    "5. What is the value of the nuclear repulsion energy?\n",
    "\n",
    "You can find the answers at the end of this notebook.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WRITE YOUR CODE BETWEEN THESE LINES - START\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# WRITE YOUR CODE BETWEEN THESE LINES - END"
   ]
  },
  {
   "cell_type": "markdown",
   "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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem\n",
    "problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "# Generate the second-quantized operators\n",
    "second_q_ops = problem.second_q_ops()\n",
    "\n",
    "# Hamiltonian\n",
    "main_op = second_q_ops[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. QubitConverter\n",
    "\n",
    "Allows to define the mapping that you will use in the simulation. You can try different mapping but \n",
    "we will stick to `JordanWignerMapper` as allows a simple correspondence: a qubit represents a spin-orbital in the molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.mappers.second_quantization import ParityMapper, BravyiKitaevMapper, JordanWignerMapper\n",
    "from qiskit_nature.converters.second_quantization.qubit_converter import QubitConverter\n",
    "\n",
    "# Setup the mapper and qubit converter\n",
    "mapper_type = 'JordanWignerMapper'\n",
    "\n",
    "if mapper_type == 'ParityMapper':\n",
    "    mapper = ParityMapper()\n",
    "elif mapper_type == 'JordanWignerMapper':\n",
    "    mapper = JordanWignerMapper()\n",
    "elif mapper_type == 'BravyiKitaevMapper':\n",
    "    mapper = BravyiKitaevMapper()\n",
    "\n",
    "converter = QubitConverter(mapper=mapper, two_qubit_reduction=False)\n",
    "\n",
    "# The fermionic operators are mapped to qubit operators\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "qubit_op = converter.convert(main_op, num_particles=num_particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. Initial state\n",
    "As we described in the Theory section, a good initial state in chemistry is the HF state (i.e. $|\\Psi_{HF} \\rangle = |0101 \\rangle$). We can initialize it as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.circuit.library import HartreeFock\n",
    "\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals\n",
    "init_state = HartreeFock(num_spin_orbitals, num_particles, converter)\n",
    "print(init_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5. Ansatz\n",
    "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.circuit.library import TwoLocal\n",
    "from qiskit_nature.circuit.library import UCCSD, PUCCD, SUCCD\n",
    "\n",
    "# Choose the ansatz\n",
    "ansatz_type = \"TwoLocal\"\n",
    "\n",
    "# Parameters for q-UCC antatze\n",
    "num_particles = (problem.molecule_data_transformed.num_alpha,\n",
    "             problem.molecule_data_transformed.num_beta)\n",
    "num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals\n",
    "\n",
    "# Put arguments for twolocal\n",
    "if ansatz_type == \"TwoLocal\":\n",
    "    # Single qubit rotations that are placed on all qubits with independent parameters\n",
    "    rotation_blocks = ['ry', 'rz']\n",
    "    # Entangling gates\n",
    "    entanglement_blocks = 'cx'\n",
    "    # How the qubits are entangled \n",
    "    entanglement = 'full'\n",
    "    # Repetitions of rotation_blocks + entanglement_blocks with independent parameters\n",
    "    repetitions = 3\n",
    "    # Skip the final rotation_blocks layer\n",
    "    skip_final_rotation_layer = True\n",
    "    ansatz = TwoLocal(qubit_op.num_qubits, rotation_blocks, entanglement_blocks, reps=repetitions, \n",
    "                      entanglement=entanglement, skip_final_rotation_layer=skip_final_rotation_layer)\n",
    "    # Add the initial state\n",
    "    ansatz.compose(init_state, front=True, inplace=True)\n",
    "elif ansatz_type == \"UCCSD\":\n",
    "    ansatz = UCCSD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"PUCCD\":\n",
    "    ansatz = PUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"SUCCD\":\n",
    "    ansatz = SUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)\n",
    "elif ansatz_type == \"Custom\":\n",
    "    # Example of how to write your own circuit\n",
    "    from qiskit.circuit import Parameter, QuantumCircuit, QuantumRegister\n",
    "    # Define the variational parameter\n",
    "    theta = Parameter('a')\n",
    "    n = qubit_op.num_qubits\n",
    "    # Make an empty quantum circuit\n",
    "    qc = QuantumCircuit(qubit_op.num_qubits)\n",
    "    qubit_label = 0\n",
    "    # Place a Hadamard gate\n",
    "    qc.h(qubit_label)\n",
    "    # Place a CNOT ladder\n",
    "    for i in range(n-1):\n",
    "        qc.cx(i, i+1)\n",
    "    # Visual separator\n",
    "    qc.barrier()\n",
    "    # rz rotations on all qubits\n",
    "    qc.rz(theta, range(n))\n",
    "    ansatz = qc\n",
    "    ansatz.compose(init_state, front=True, inplace=True)\n",
    "\n",
    "print(ansatz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6. Backend\n",
    "This is where you specify the simulator or device where you want to run your algorithm.\n",
    "We will focus on the `statevector_simulator` in this challenge.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import Aer\n",
    "backend = Aer.get_backend('statevector_simulator')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7. Optimizer\n",
    "\n",
    "The optimizer guides the evolution of the parameters of the ansatz so it is very important to investigate the energy convergence as it would define the number of measurements that have to be performed on the QPU.\n",
    "A clever choice might reduce drastically the number of needed energy evaluations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA, SLSQP\n",
    "\n",
    "optimizer_type = 'COBYLA'\n",
    "\n",
    "# You may want to tune the parameters \n",
    "# of each optimizer, here the defaults are used\n",
    "if optimizer_type == 'COBYLA':\n",
    "    optimizer = COBYLA(maxiter=500)\n",
    "elif optimizer_type == 'L_BFGS_B':\n",
    "    optimizer = L_BFGS_B(maxfun=500)\n",
    "elif optimizer_type == 'SPSA':\n",
    "    optimizer = SPSA(maxiter=500)\n",
    "elif optimizer_type == 'SLSQP':\n",
    "    optimizer = SLSQP(maxiter=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8. Exact eigensolver\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.\n",
    "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. \n",
    "For very large systems you would run out of memory trying to store their wavefunctions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.algorithms.ground_state_solvers.minimum_eigensolver_factories import NumPyMinimumEigensolverFactory\n",
    "from qiskit_nature.algorithms.ground_state_solvers import GroundStateEigensolver\n",
    "import numpy as np \n",
    "\n",
    "def exact_diagonalizer(problem, converter):\n",
    "    solver = NumPyMinimumEigensolverFactory()\n",
    "    calc = GroundStateEigensolver(converter, solver)\n",
    "    result = calc.solve(problem)\n",
    "    return result\n",
    "\n",
    "result_exact = exact_diagonalizer(problem, converter)\n",
    "exact_energy = np.real(result_exact.eigenenergies[0])\n",
    "print(\"Exact electronic energy\", exact_energy)\n",
    "print(result_exact)\n",
    "\n",
    "# The targeted electronic energy for H2 is -1.85336 Ha\n",
    "# Check with your VQE result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 9. VQE and initial parameters for the ansatz\n",
    "Now we can import the VQE class and run the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms import VQE\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# Print and save the data in lists\n",
    "def callback(eval_count, parameters, mean, std):  \n",
    "    # Overwrites the same line when printing\n",
    "    display(\"Evaluation: {}, Energy: {}, Std: {}\".format(eval_count, mean, std))\n",
    "    clear_output(wait=True)\n",
    "    counts.append(eval_count)\n",
    "    values.append(mean)\n",
    "    params.append(parameters)\n",
    "    deviation.append(std)\n",
    "\n",
    "counts = []\n",
    "values = []\n",
    "params = []\n",
    "deviation = []\n",
    "\n",
    "# Set initial parameters of the ansatz\n",
    "# We choose a fixed small displacement \n",
    "# So all participants start from similar starting point\n",
    "try:\n",
    "    initial_point = [0.01] * len(ansatz.ordered_parameters)\n",
    "except:\n",
    "    initial_point = [0.01] * ansatz.num_parameters\n",
    "\n",
    "algorithm = VQE(ansatz,\n",
    "                optimizer=optimizer,\n",
    "                quantum_instance=backend,\n",
    "                callback=callback,\n",
    "                initial_point=initial_point)\n",
    "\n",
    "result = algorithm.compute_minimum_eigenvalue(qubit_op)\n",
    "\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 9. Scoring function \n",
    "We need to judge how good are your VQE simulations, your choice of ansatz/optimizer.\n",
    "For this, we implemented the following simple scoring function:\n",
    "\n",
    "$$ score = N_{CNOT}$$\n",
    "\n",
    "where $N_{CNOT}$ is the number of CNOTs. \n",
    "But you have to reach the chemical accuracy which is $\\delta E_{chem} = 0.004$ Ha $= 4$ mHa, which may be hard to reach depending on the problem. \n",
    "You have to reach the accuracy we set in a minimal number of CNOTs to win the challenge. \n",
    "The lower the score the better!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store results in a dictionary\n",
    "from qiskit.transpiler import PassManager\n",
    "from qiskit.transpiler.passes import Unroller\n",
    "\n",
    "# Unroller transpile your circuit into CNOTs and U gates\n",
    "pass_ = Unroller(['u', 'cx'])\n",
    "pm = PassManager(pass_)\n",
    "ansatz_tp = pm.run(ansatz)\n",
    "cnots = ansatz_tp.count_ops()['cx']\n",
    "score = cnots\n",
    "\n",
    "accuracy_threshold = 4.0 # in mHa\n",
    "energy = result.optimal_value\n",
    "\n",
    "if ansatz_type == \"TwoLocal\":\n",
    "    result_dict = {\n",
    "        'optimizer': optimizer.__class__.__name__,\n",
    "        'mapping': converter.mapper.__class__.__name__,\n",
    "        'ansatz': ansatz.__class__.__name__,\n",
    "        'rotation blocks': rotation_blocks,\n",
    "        'entanglement_blocks': entanglement_blocks,\n",
    "        'entanglement': entanglement,\n",
    "        'repetitions': repetitions,\n",
    "        'skip_final_rotation_layer': skip_final_rotation_layer,\n",
    "        'energy (Ha)': energy,\n",
    "        'error (mHa)': (energy-exact_energy)*1000,\n",
    "        'pass': (energy-exact_energy)*1000 <= accuracy_threshold,\n",
    "        '# of parameters': len(result.optimal_point),\n",
    "        'final parameters': result.optimal_point,\n",
    "        '# of evaluations': result.optimizer_evals,\n",
    "        'optimizer time': result.optimizer_time,\n",
    "        '# of qubits': int(qubit_op.num_qubits),\n",
    "        '# of CNOTs': cnots,\n",
    "        'score': score}\n",
    "else:\n",
    "    result_dict = {\n",
    "        'optimizer': optimizer.__class__.__name__,\n",
    "        'mapping': converter.mapper.__class__.__name__,\n",
    "        'ansatz': ansatz.__class__.__name__,\n",
    "        'rotation blocks': None,\n",
    "        'entanglement_blocks': None,\n",
    "        'entanglement': None,\n",
    "        'repetitions': None,\n",
    "        'skip_final_rotation_layer': None,\n",
    "        'energy (Ha)': energy,\n",
    "        'error (mHa)': (energy-exact_energy)*1000,\n",
    "        'pass': (energy-exact_energy)*1000 <= accuracy_threshold,\n",
    "        '# of parameters': len(result.optimal_point),\n",
    "        'final parameters': result.optimal_point,\n",
    "        '# of evaluations': result.optimizer_evals,\n",
    "        'optimizer time': result.optimizer_time,\n",
    "        '# of qubits': int(qubit_op.num_qubits),\n",
    "        '# of CNOTs': cnots,\n",
    "        'score': score}\n",
    "\n",
    "# Plot the results\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.set_xlabel('Iterations')\n",
    "ax.set_ylabel('Energy')\n",
    "ax.grid()\n",
    "fig.text(0.7, 0.75, f'Energy: {result.optimal_value:.3f}\\nScore: {score:.0f}')\n",
    "plt.title(f\"{result_dict['optimizer']}-{result_dict['mapping']}\\n{result_dict['ansatz']}\")\n",
    "ax.plot(counts, values)\n",
    "ax.axhline(exact_energy, linestyle='--')\n",
    "fig_title = f\"\\\n",
    "{result_dict['optimizer']}-\\\n",
    "{result_dict['mapping']}-\\\n",
    "{result_dict['ansatz']}-\\\n",
    "Energy({result_dict['energy (Ha)']:.3f})-\\\n",
    "Score({result_dict['score']:.0f})\\\n",
    ".png\"\n",
    "fig.savefig(fig_title, dpi=300)\n",
    "\n",
    "# Display and save the data\n",
    "import pandas as pd\n",
    "import os.path\n",
    "filename = 'results_h2.csv'\n",
    "if os.path.isfile(filename):\n",
    "    result_df = pd.read_csv(filename)\n",
    "    result_df = result_df.append([result_dict])\n",
    "else:\n",
    "    result_df = pd.DataFrame.from_dict([result_dict])\n",
    "result_df.to_csv(filename)\n",
    "result_df[['optimizer','ansatz', '# of qubits', '# of parameters','rotation blocks', 'entanglement_blocks',\n",
    "    'entanglement', 'repetitions', 'error (mHa)', 'pass', 'score']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "<b>Tutorial questions 2</b> \n",
    "\n",
    "Experiment with all the parameters and then:\n",
    "\n",
    "1. Can you find your best (best score) heuristic ansatz (by modifying parameters of `TwoLocal` ansatz) and optimizer?\n",
    "2. Can you find your best q-UCC ansatz (choose among `UCCSD, PUCCD or SUCCD` ansatzes) and optimizer?\n",
    "3. In the cell where we define the ansatz, \n",
    "   can you modify the `Custom` ansatz by placing gates yourself to write a better circuit than your `TwoLocal` circuit? \n",
    "\n",
    "For each question, give `ansatz` objects.\n",
    "Remember, you have to reach the chemical accuracy $|E_{exact} - E_{VQE}| \\leq 0.004 $ Ha $= 4$ mHa.\n",
    "    \n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WRITE YOUR CODE BETWEEN THESE LINES - START\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# WRITE YOUR CODE BETWEEN THESE LINES - END"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Final Challenge - VQE for LiH molecule \n",
    "\n",
    "\n",
    "In this part, you will simulate LiH molecule using the STO-3G basis with the PySCF driver.\n",
    "\n",
    "</div>\n",
    "    \n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "<b>Goal</b> \n",
    "\n",
    "Experiment with all the parameters and then find your best ansatz. You can be as creative as you want!\n",
    "\n",
    "For each question, give `ansatz` objects as for Part 1. Your final score will be based only on Part 2.\n",
    "    \n",
    "</div>\n",
    "\n",
    "Be aware that the system is larger now. Work out how many qubits you would need for this system by retrieving the number of spin-orbitals. \n",
    "\n",
    "### Reducing the problem size\n",
    "\n",
    "You might want to reduce the number of qubits for your simulation:\n",
    "- you could freeze the core electrons that do not contribute significantly to chemistry and consider only the valence electrons. Qiskit  already has this functionality implemented. So inspect the different transformers in `qiskit_nature.transformers`  and find the one that performs the freeze core approximation.\n",
    "- you could use `ParityMapper` with `two_qubit_reduction=True` to eliminate 2 qubits.\n",
    "- you could reduce the number of qubits by inspecting the symmetries of your Hamiltonian. Find a way to use `Z2Symmetries` in Qiskit.\n",
    "\n",
    "### Custom ansatz \n",
    "\n",
    "You might want to explore the ideas proposed in [Grimsley *et al.*,2018](https://arxiv.org/abs/1812.11173v2), [H. L. Tang *et al.*,2019](https://arxiv.org/abs/1911.10205), [Rattew *et al.*,2019](https://arxiv.org/abs/1910.09694), [Tang *et al.*,2019](https://arxiv.org/abs/1911.10205). \n",
    "You can even get try machine learning algorithms to generate best ansatz circuits.\n",
    "\n",
    "### Setup the simulation\n",
    "\n",
    "Let's now run the Hartree-Fock calculation and the rest is up to you!\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "<b>Attention</b> \n",
    "\n",
    "We give below the `driver`, the `initial_point`, the `initial_state` that should remain as given.\n",
    "You are free then to explore all other things available in Qiskit.\n",
    "So you have to start from this initial point (all parameters set to 0.01):\n",
    "    \n",
    "`initial_point = [0.01] * len(ansatz.ordered_parameters)`\n",
    "    or\n",
    "`initial_point = [0.01] * ansatz.num_parameters`\n",
    "\n",
    "and your initial state has to be the Hartree-Fock state:\n",
    "    \n",
    "`init_state = HartreeFock(num_spin_orbitals, num_particles, converter)`\n",
    "    \n",
    "For each question, give `ansatz` object.\n",
    "Remember you have to reach the chemical accuracy $|E_{exact} - E_{VQE}| \\leq 0.004 $ Ha $= 4$ mHa.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers import PySCFDriver\n",
    "\n",
    "molecule = 'Li 0.0 0.0 0.0; H 0.0 0.0 1.5474'\n",
    "driver = PySCFDriver(atom=molecule)\n",
    "qmolecule = driver.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WRITE YOUR CODE BETWEEN THESE LINES - START\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# WRITE YOUR CODE BETWEEN THESE LINES - END"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer using following code\n",
    "from qc_grader import grade_ex5\n",
    "freeze_core = False # change to True if you freezed core electrons\n",
    "grade_ex5(ansatz,qubit_op,result,freeze_core)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submit your answer. You can re-submit at any time.\n",
    "from qc_grader import submit_ex5\n",
    "submit_ex5(ansatz,qubit_op,result,freeze_core)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answers for Part 1\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "<b>Questions</b> \n",
    "    \n",
    "Look into the attributes of `qmolecule` and answer the questions below.\n",
    "\n",
    "    \n",
    "1. We need to know the basic characteristics of our molecule. What is the total number of electrons in your system?\n",
    "2. What is the number of molecular orbitals?\n",
    "3. What is the number of spin-orbitals?\n",
    "3. How many qubits would you need to simulate this molecule with Jordan-Wigner mapping?\n",
    "5. What is the value of the nuclear repulsion energy?\n",
    "    \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "<b>Answers </b> \n",
    "\n",
    "1. `n_el = qmolecule.num_alpha + qmolecule.num_beta`\n",
    "    \n",
    "2. `n_mo = qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "3. `n_so = 2 * qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "4. `n_q = 2* qmolecule.num_molecular_orbitals`\n",
    "    \n",
    "5. `e_nn = qmolecule.nuclear_repulsion_energy`\n",
    "    \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "**Created by:** Igor Sokolov, Junye Huang, Rahul Pratap Singh\n",
    "\n",
    "**Version:** 1.0.1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}