{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to OpenFermion\n",
    "\n",
    "In this tutorial we will learn how to use OpenFermion to represent and manipulate strongly-correlated Hamiltonians, and to interface it with Cirq.\n",
    "To install OpenFermion, please see the instructions at https://github.com/quantumlib/OpenFermion\n",
    "To install the OpenFermion-Cirq interface, please see the instructions at https://github.com/quantumlib/OpenFermion-Cirq.\n",
    "\n",
    "(Some of the code in this tutorial is adapted from tutorial 4 of OpenFermion-Cirq.) \n",
    "\n",
    "1. What is OpenFermion?\n",
    "1. Operators\n",
    "    1. The QubitOperator class\n",
    "    1. The FermionOperator class\n",
    "    1. The Jordan-Wigner transformation\n",
    "1. Quantum Chemistry\n",
    "    1. The MolecularData class\n",
    "    1. The InteractionOperator class\n",
    "1. The OpenFermion-Cirq interface\n",
    "1. Exercise: the unitary coupled cluster ansatz\n",
    "    \n",
    "\n",
    "# What is OpenFermion?\n",
    "\n",
    "OpenFermion is a package that allows for the easy representation and manipulation of operators on strongly-correlated fermionic, bosonic, and qubit systems. This in turn allows it to act as an interface between computational chemistry software and quantum circuit design packages such as cirq. It also contains various methods for generating quantum circuits of use for quantum chemistry/strongly correlated materials.\n",
    "\n",
    "# Operators\n",
    "\n",
    "OpenFermion provides many convenient methods for representing the various operators that appear in strongly-correlated physics and chemistry problems. In particular, OpenFermion aims for these operators to be human-readable, and to not require storing an exponentially large matrix.\n",
    "\n",
    "## The QubitOperator class\n",
    "\n",
    "We have repeatedly encountered the Pauli operators on $N$ qubits, $\\{I,X,Y,Z\\}^{\\otimes N}$. OpenFermion represents these, and linear combinations of Pauli operators, using the QubitOperator class. Internally, the data about the operator is stored as a dictionary, with the names of individual Pauli operators used as keys.\n",
    "\n",
    "QubitOperators are initialized as single Pauli terms separated by space, and each Pauli term is written index-last. An optional coefficient may be provided to multiply the term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0 [X1 Y2]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import QubitOperator\n",
    "qop1 = QubitOperator(\"X1 Y2\", 3.0)\n",
    "print(qop1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QubitOperators may be added, multiplied, and multiplied by constants to produce other QubitOperators "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0 [X1 Y2] +\n",
      "3.0 [Z1 Z2]\n",
      "(3+0j) [Y1 X2]\n"
     ]
    }
   ],
   "source": [
    "qop2 = QubitOperator(\"Z1 Z2\")\n",
    "print(qop1 + 3 * qop2)\n",
    "print(qop1 * qop2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1:** Verify that QubitOperators on a single qubit obey the Pauli operator commutation relations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "qopx0 = QubitOperator('X0', 1.0)\n",
    "qopy0 = QubitOperator('Y0', 1.0)\n",
    "qopz0 = QubitOperator('Z0', 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Square to identity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 []"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopx0 * qopx0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Commutation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2j [X0]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopy0 * qopz0 - qopz0 * qopy0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anti-commutation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopy0 * qopz0 + qopz0 * qopy0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Paulis in acting on different qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "qopx1 = QubitOperator('X1', 1.0)\n",
    "qopy1 = QubitOperator('Y1', 1.0)\n",
    "qopz1 = QubitOperator('Z1', 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Product of Paulis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 [X0 X1]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopx0 * qopx1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Commutator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopz0 * qopy1 - qopy1 * qopz0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anti-commutator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0 [Z0 Y1]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qopz0 * qopy1 + qopy1 * qopz0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenFermion additionally has methods to convert QubitOperators into sparse matrix form, or to calculate the eigenspectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-6. -6.  0.  0.  0.  0.  6.  6.]\n",
      "\n",
      "[[ 3.+0.j  0.+0.j  0.+0.j  0.-3.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j -3.+0.j  0.+3.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.-3.j -3.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+3.j  0.+0.j  0.+0.j  3.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  3.+0.j  0.+0.j  0.+0.j  0.-3.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j -3.+0.j  0.+3.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.-3.j -3.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+3.j  0.+0.j  0.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import eigenspectrum, get_sparse_operator\n",
    "print(eigenspectrum(qop1 + 3*qop2))\n",
    "matrix = get_sparse_operator(qop1 + 3*qop2).todense()\n",
    "print()\n",
    "print(matrix.round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** Verify that the matrix produced is identical to the matrix obtained using numpy's kron function, and check that the eigenvalues match those produced by eigenspectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "px = numpy.array([[0, 1],\n",
    "                  [1, 0]])\n",
    "py= numpy.array([[0, -1j],\n",
    "                 [1j, 0]])\n",
    "pz = numpy.array([[1,  0],\n",
    "                  [0, -1]])\n",
    "p0 = numpy.identity(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "qop1_mat = 3 * numpy.kron(p0, numpy.kron(px, py))\n",
    "qop2_mat = numpy.kron(p0, numpy.kron(pz, pz))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix2 = qop1_mat + 3 * qop2_mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OpenFermion and numpy matrices are the same:  True\n"
     ]
    }
   ],
   "source": [
    "print('OpenFermion and numpy matrices are the same: ', numpy.allclose(matrix, matrix2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues of numpy matrix:  [-6. -6.  0.  0.  0.  0.  6.  6.]\n"
     ]
    }
   ],
   "source": [
    "print('Eigenvalues of numpy matrix: ', numpy.linalg.eigvalsh(matrix2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The FermionOperator class\n",
    "\n",
    "OpenFermion can also represent operators on fermionic systems via the FermionicOperator class. This stores a representation of an operator as a sum of creation and annihilation operators. Individual FermionOperators are instantiated in a similar way to QubitOperators, but instead of needing to write a Pauli operator, we simply write the indices, and represent creation operators with a '^' symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 [2 1]\n",
      "\n",
      "6 [2 1 2^ 3]\n",
      "\n",
      "2 [2 1] +\n",
      "9.0 [2^ 3]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import FermionOperator\n",
    "\n",
    "fop1 = FermionOperator('2 1', 2)\n",
    "print(fop1)\n",
    "fop2 = FermionOperator('2^ 3', 3)\n",
    "print()\n",
    "print(fop1 * fop2)\n",
    "print()\n",
    "print(fop1 + 3 * fop2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FermionOperators introduce a new issue - individual creation and annihilation operators no longer commute, so there is not an immediately obvious choice for the 'standard form' that we would always want to write a product of these operators in --- nor is it immediately obvious whether two FermionOperators are equal!\n",
    "\n",
    "One representation that is standard is 'normal ordering', where we first write all creation operators in decreasing order of their index and then all annihilation operators in decreasing order of their index. To convert a FermionOperator into its normal-ordered form, OpenFermion provides the normal_ordered function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 [2 1 2^ 3]\n",
      "\n",
      "6.0 [2^ 3 2 1] +\n",
      "6.0 [3 1]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import normal_ordered\n",
    "print(fop1 * fop2)\n",
    "print()\n",
    "print(normal_ordered(fop1 * fop2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that what was a single term in its original order becomes a sum of two terms when converting to normal order! This is why normal order is not always preferable - the number of terms can grow quite large. Normal-ordering is also quite computationally costly to perform, and so openfermion does not perform it at each step. This implies that equality testing between FermionicOperators requires you to perform normal ordering first!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n",
      "1.0 [1 2] +\n",
      "1.0 [2 1]\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "fop1 = FermionOperator('1 2', 1.0)\n",
    "fop2 = FermionOperator('2 1', -1.0)\n",
    "print(fop1 == fop2)\n",
    "print(normal_ordered(fop1) == normal_ordered(fop2))\n",
    "print(fop1 - fop2)\n",
    "print(normal_ordered(fop1 - fop2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** Verify that FermionOperators obey the correct anti-commutation rules on a system with 2 indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion import anticommutator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "fop3 = FermionOperator('0^ 1', 1)\n",
    "fop4 = FermionOperator('1 2^', 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1 [0^ 1 1 2^] +\n",
       "1.0 [1 2^ 0^ 1]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "anticommutator(fop3, fop4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1 [0^ 1 1 2^] +\n",
       "1.0 [1 2^ 0^ 1]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fop3 * fop4 + fop4 * fop3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normal_ordered(fop3 * fop4 + fop4 * fop3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "fop5 = FermionOperator('0^ 1')\n",
    "fop6 = FermionOperator('1^ 0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 [0^ 1 1^ 0] +\n",
       "1.0 [1^ 0 0^ 1]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(fop5 * fop6 + fop6 * fop5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 [0^ 0] +\n",
       "2.0 [1^ 0^ 1 0] +\n",
       "1.0 [1^ 1]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normal_ordered(fop5 * fop6 + fop6 * fop5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above operations to generate a sparse matrix and eigenspectrum work just as well for fermion operators as they do qubit operators. However, remember that creation and annihilation operators are not themselves hermitian. To create an observable operator with real eigenspectrum, we should sum our fermion operators with their hermitian conjugates. OpenFermion provides the functions hermitian_conjugated and is_hermitian to help with this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "\n",
      "1.0 [1 2] +\n",
      "1.0 [2^ 1^]\n",
      "\n",
      "True\n",
      "\n",
      "[[ 0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [-1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j]]\n",
      "\n",
      "[-1. -1.  0.  0.  0.  0.  1.  1.]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import hermitian_conjugated, is_hermitian\n",
    "\n",
    "print(is_hermitian(fop1))\n",
    "print()\n",
    "hermitian_op = fop1 + hermitian_conjugated(fop1)\n",
    "print(hermitian_op)\n",
    "print()\n",
    "print(is_hermitian(hermitian_op))\n",
    "print()\n",
    "matrix = get_sparse_operator(hermitian_op)\n",
    "print(matrix.todense().round(2))\n",
    "print()\n",
    "print(eigenspectrum(fop1 + hermitian_conjugated(fop1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4:** Investigate the matrix representation of the sum of a single creation and single annihilation operator. What does this remind you of?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "ferop0 = FermionOperator('0') + FermionOperator('0^')\n",
    "ferop1 = FermionOperator('1') + FermionOperator('1^')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_ferop0 = get_sparse_operator(ferop0)\n",
    "matrix_ferop1 = get_sparse_operator(ferop1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.+0.j 1.+0.j]\n",
      " [1.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "print(matrix_ferop0.todense().round(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]\n",
      " [ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j -1.+0.j]\n",
      " [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "print(matrix_ferop1.todense().round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those matrices look a lot like Pauli operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is operator 1.0 [0] +\n",
      "1.0 [0^] \n",
      " the same matrix as X: True\n"
     ]
    }
   ],
   "source": [
    "print('Is operator {} \\n the same matrix as X: {}'.format(ferop0, \n",
    "                                                          numpy.allclose(matrix_ferop0.todense(), px)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is operator 1.0 [1] +\n",
      "1.0 [1^] \n",
      " the same matrix as Z0X1: True\n"
     ]
    }
   ],
   "source": [
    "print('Is operator {} \\n the same matrix as Z0X1: {}'.format(ferop1, \n",
    "                                                          numpy.allclose(matrix_ferop1.todense(),\n",
    "                                                                         numpy.kron(pz, px))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Jordan-Wigner transformation\n",
    "\n",
    "Transforming between fermionic and qubit operators is a somewhat complicated task, as one needs to preserve the commutation relations of individual operators. Many choices for this transformation exist, all with various advantages and disadvantages. We will study this shortly in both the lectures and tutorials, but for the sake of this tutorial we will just use one possible transformation - the Jordan-Wigner transformation.\n",
    "\n",
    "In short, the Jordan-Wigner transformation maps $$\\hat{c}_i\\rightarrow\\frac{1}{2}\\otimes_{j<i}Z_j\\otimes(X_i+iY_i)$$\n",
    "$$\\hat{c}^{\\dagger}_i\\rightarrow\\frac{1}{2}\\otimes_{j<i}Z_j\\otimes(X_i-iY_i)$$\n",
    "This is implemented in OpenFermion by the jordan_wigner function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mapping annihilation operator\n",
      "0.5 [Z0 Z1 Z2 Z3 Z4 X5] +\n",
      "0.5j [Z0 Z1 Z2 Z3 Z4 Y5]\n",
      "\n",
      "Mapping creation operator\n",
      "0.5 [Z0 Z1 Z2 Z3 Z4 X5] +\n",
      "-0.5j [Z0 Z1 Z2 Z3 Z4 Y5]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import jordan_wigner\n",
    "print('Mapping annihilation operator')\n",
    "print(jordan_wigner(FermionOperator('5')))\n",
    "print()\n",
    "print('Mapping creation operator')\n",
    "print(jordan_wigner(FermionOperator('5^')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5:** Check that the Jordan-Wigner transform preserves the fermionic commutation relations on a $2$ site system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "fop0 = FermionOperator('0', 1.0)\n",
    "fop0dag = FermionOperator('0^', 1.0)\n",
    "fop1 = FermionOperator('1', 1.0)\n",
    "fop1dag = FermionOperator('1^', 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "qop0 = jordan_wigner(fop0)\n",
    "qop0dag = jordan_wigner(fop0dag)\n",
    "qop1 = jordan_wigner(fop1)\n",
    "qop1dag = jordan_wigner(fop1dag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fermionic anticommutator in site 0:  1.0 []\n",
      "Jordan-Wigner transformed anticommutator:  (1+0j) []\n",
      "\n",
      "Fermionic anticommutator in site 1:  1.0 []\n",
      "Jordan-Wigner transformed anticommutator:  (1+0j) []\n"
     ]
    }
   ],
   "source": [
    "# The anticommutator of creation annihilation on same site\n",
    "print('Fermionic anticommutator in site 0: ', normal_ordered(anticommutator(fop0, fop0dag)))\n",
    "print('Jordan-Wigner transformed anticommutator: ', anticommutator(qop0, qop0dag))\n",
    "print()\n",
    "print('Fermionic anticommutator in site 1: ', normal_ordered(anticommutator(fop1, fop1dag)))\n",
    "print('Jordan-Wigner transformed anticommutator: ', anticommutator(qop1, qop1dag))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fermionic anticommutator 0 1^:  0\n",
      "Jordan-Wigner transformed anticommutator:  0\n",
      "\n",
      "Fermionic anticommutator 0^ 1^:  0\n",
      "Jordan-Wigner transformed anticommutator:  0\n"
     ]
    }
   ],
   "source": [
    "# The anticommutator of creation annihilation on different sites\n",
    "print('Fermionic anticommutator 0 1^: ', normal_ordered(anticommutator(fop0, fop1dag)))\n",
    "print('Jordan-Wigner transformed anticommutator: ', anticommutator(qop0, qop1dag))\n",
    "print()\n",
    "print('Fermionic anticommutator 0^ 1^: ', normal_ordered(anticommutator(fop0dag, fop1dag)))\n",
    "print('Jordan-Wigner transformed anticommutator: ', anticommutator(qop0dag, qop1dag))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We leave as an exercise to check the rest of anticommutation relations between 0 and 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 6:** Recall the conditions for a QubitOperator to be Hermitian, and check that this is preserved under the Jordan-Wigner transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion import hermitian_conjugated, is_hermitian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is 0.5 [X0] +\n",
      "0.5j [Y0] hermitian: False\n"
     ]
    }
   ],
   "source": [
    "# A QubitOperator is Hermitian if its coefficients are real\n",
    "print('Is {} hermitian: {}'.format(qop0, is_hermitian(qop0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is 1.0 [X0] hermitian: True\n"
     ]
    }
   ],
   "source": [
    "print('Is {} hermitian: {}'.format(qop0+hermitian_conjugated(qop0),\n",
    "                                   is_hermitian(qop0+hermitian_conjugated(qop0))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 7:** Check that the Jordan-Wigner transform preserves the spectrum of the operator $\\hat{c}^{\\dagger}_0\\hat{c}_1+\\hat{c}^{\\dagger}_1\\hat{c}_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 [0^ 1] +\n",
      "1.0 [1^ 0]\n",
      "\n",
      "(0.5+0j) [X0 X1] +\n",
      "(0.5+0j) [Y0 Y1]\n"
     ]
    }
   ],
   "source": [
    "toy_hamiltonian = fop0dag*fop1 + fop1dag*fop0\n",
    "toy_ham_jw = jordan_wigner(toy_hamiltonian)\n",
    "print(toy_hamiltonian)\n",
    "print()\n",
    "print(toy_ham_jw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spectrum of Fermionic toy hamiltonian: [-1.  0.  0.  1.]\n"
     ]
    }
   ],
   "source": [
    "print('Spectrum of Fermionic toy hamiltonian: {}'.format(eigenspectrum(toy_hamiltonian)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spectrum of Fermionic toy hamiltonian: [-1.  0.  0.  1.]\n"
     ]
    }
   ],
   "source": [
    "print('Spectrum of Fermionic toy hamiltonian: {}'.format(eigenspectrum(toy_ham_jw)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Chemistry\n",
    "\n",
    "One of the most critical functions of OpenFermion is its ability to interface with computational chemistry packages such as psi4, pySCF, and DIRAC, which are then used to provide data (e.g. Hamiltonians) on molecular systems. This is done through various interface packages; due to the difficulty in installing these packages we will not use them here. However, OpenFermion ships with some molecular system data already, which is enough for us to use.\n",
    "\n",
    "## The MolecularData class\n",
    "\n",
    "OpenFermion uses the MolecularData class to store data from computational chemistry packages. This class has additionally a load and save function, as this data often takes a significant amount of time to prepare. To load stored data, we need to state the geometry of the molecule (in cartesian co-ordinates), the basis in which the molecule is stored, and its charge and multiplicity (i.e. $2n+1$, where $n$ is the number of unpaired spins). OpenFermion also requests a description of the molecule studied"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H2_sto-3g_singlet_0.7414\n"
     ]
    }
   ],
   "source": [
    "from openfermion import MolecularData\n",
    "\n",
    "diatomic_bond_length = .7414\n",
    "\n",
    "geometry = [('H', (0., 0., 0.)), \n",
    "            ('H', (0., 0., diatomic_bond_length))]\n",
    "\n",
    "basis = 'sto-3g'\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "description = format(diatomic_bond_length)\n",
    "\n",
    "molecule = MolecularData(\n",
    "    geometry,\n",
    "    basis,\n",
    "    multiplicity,\n",
    "    description=description)\n",
    "molecule.load()\n",
    "\n",
    "print(molecule.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the above data, and being able to generate the Hamiltonian of the target molecule, the MolecularData class contains various information about classical computational approximations of the ground state energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond Length in Angstroms: 0.7414\n",
      "Hartree Fock (mean-field) energy in Hartrees: -1.116684386906734\n",
      "FCI (Exact) energy in Hartrees: -1.137270174625328\n"
     ]
    }
   ],
   "source": [
    "print(\"Bond Length in Angstroms: {}\".format(diatomic_bond_length))\n",
    "print(\"Hartree Fock (mean-field) energy in Hartrees: {}\".format(molecule.hf_energy))\n",
    "print(\"FCI (Exact) energy in Hartrees: {}\".format(molecule.fci_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The InteractionOperator class\n",
    "\n",
    "The Hamiltonian of the above system may be obtained via the get_molecular_hamiltonian function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'openfermion.ops._interaction_operator.InteractionOperator'>\n",
      "() 0.713753990544915\n",
      "((0, 1), (0, 0)) -1.2524635715927757\n",
      "((1, 1), (1, 0)) -1.2524635715927757\n",
      "((2, 1), (2, 0)) -0.4759487172683097\n",
      "((3, 1), (3, 0)) -0.4759487172683097\n",
      "((0, 1), (0, 1), (0, 0), (0, 0)) 0.3372443828669511\n",
      "((0, 1), (0, 1), (2, 0), (2, 0)) 0.09064440419713085\n",
      "((0, 1), (1, 1), (1, 0), (0, 0)) 0.3372443828669511\n",
      "((0, 1), (1, 1), (3, 0), (2, 0)) 0.09064440419713085\n",
      "((0, 1), (2, 1), (0, 0), (2, 0)) 0.09064440419713081\n",
      "((0, 1), (2, 1), (2, 0), (0, 0)) 0.33173404792821903\n",
      "((0, 1), (3, 1), (1, 0), (2, 0)) 0.09064440419713081\n",
      "((0, 1), (3, 1), (3, 0), (0, 0)) 0.33173404792821903\n",
      "((1, 1), (0, 1), (0, 0), (1, 0)) 0.3372443828669511\n",
      "((1, 1), (0, 1), (2, 0), (3, 0)) 0.09064440419713085\n",
      "((1, 1), (1, 1), (1, 0), (1, 0)) 0.3372443828669511\n",
      "((1, 1), (1, 1), (3, 0), (3, 0)) 0.09064440419713085\n",
      "((1, 1), (2, 1), (0, 0), (3, 0)) 0.09064440419713081\n",
      "((1, 1), (2, 1), (2, 0), (1, 0)) 0.33173404792821903\n",
      "((1, 1), (3, 1), (1, 0), (3, 0)) 0.09064440419713081\n",
      "((1, 1), (3, 1), (3, 0), (1, 0)) 0.33173404792821903\n",
      "((2, 1), (0, 1), (0, 0), (2, 0)) 0.33173404792821914\n",
      "((2, 1), (0, 1), (2, 0), (0, 0)) 0.0906444041971308\n",
      "((2, 1), (1, 1), (1, 0), (2, 0)) 0.33173404792821914\n",
      "((2, 1), (1, 1), (3, 0), (0, 0)) 0.0906444041971308\n",
      "((2, 1), (2, 1), (0, 0), (0, 0)) 0.09064440419713085\n",
      "((2, 1), (2, 1), (2, 0), (2, 0)) 0.34869688341114263\n",
      "((2, 1), (3, 1), (1, 0), (0, 0)) 0.09064440419713085\n",
      "((2, 1), (3, 1), (3, 0), (2, 0)) 0.34869688341114263\n",
      "((3, 1), (0, 1), (0, 0), (3, 0)) 0.33173404792821914\n",
      "((3, 1), (0, 1), (2, 0), (1, 0)) 0.0906444041971308\n",
      "((3, 1), (1, 1), (1, 0), (3, 0)) 0.33173404792821914\n",
      "((3, 1), (1, 1), (3, 0), (1, 0)) 0.0906444041971308\n",
      "((3, 1), (2, 1), (0, 0), (1, 0)) 0.09064440419713085\n",
      "((3, 1), (2, 1), (2, 0), (3, 0)) 0.34869688341114263\n",
      "((3, 1), (3, 1), (1, 0), (1, 0)) 0.09064440419713085\n",
      "((3, 1), (3, 1), (3, 0), (3, 0)) 0.34869688341114263\n",
      "\n"
     ]
    }
   ],
   "source": [
    "hamiltonian = molecule.get_molecular_hamiltonian()\n",
    "print(hamiltonian.__class__)\n",
    "print(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above operator is returned neither as a FermionOperator nor a QubitOperator, but an InteractionOperator. This is similar to a FermionOperator, but has the restriction that it contains only one-body and two-body terms that conserve particle number and spin. This allows for various speedups in calculation not available to the more general class.\n",
    "\n",
    "Convering an InteractionOperator to a FermionOperator is performed using the get_fermion_operator function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.713753990544915 [] +\n",
      "-1.2524635715927757 [0^ 0] +\n",
      "0.3372443828669511 [0^ 0^ 0 0] +\n",
      "0.09064440419713085 [0^ 0^ 2 2] +\n",
      "0.3372443828669511 [0^ 1^ 1 0] +\n",
      "0.09064440419713085 [0^ 1^ 3 2] +\n",
      "0.09064440419713081 [0^ 2^ 0 2] +\n",
      "0.33173404792821903 [0^ 2^ 2 0] +\n",
      "0.09064440419713081 [0^ 3^ 1 2] +\n",
      "0.33173404792821903 [0^ 3^ 3 0] +\n",
      "0.3372443828669511 [1^ 0^ 0 1] +\n",
      "0.09064440419713085 [1^ 0^ 2 3] +\n",
      "-1.2524635715927757 [1^ 1] +\n",
      "0.3372443828669511 [1^ 1^ 1 1] +\n",
      "0.09064440419713085 [1^ 1^ 3 3] +\n",
      "0.09064440419713081 [1^ 2^ 0 3] +\n",
      "0.33173404792821903 [1^ 2^ 2 1] +\n",
      "0.09064440419713081 [1^ 3^ 1 3] +\n",
      "0.33173404792821903 [1^ 3^ 3 1] +\n",
      "0.33173404792821914 [2^ 0^ 0 2] +\n",
      "0.0906444041971308 [2^ 0^ 2 0] +\n",
      "0.33173404792821914 [2^ 1^ 1 2] +\n",
      "0.0906444041971308 [2^ 1^ 3 0] +\n",
      "-0.4759487172683097 [2^ 2] +\n",
      "0.09064440419713085 [2^ 2^ 0 0] +\n",
      "0.34869688341114263 [2^ 2^ 2 2] +\n",
      "0.09064440419713085 [2^ 3^ 1 0] +\n",
      "0.34869688341114263 [2^ 3^ 3 2] +\n",
      "0.33173404792821914 [3^ 0^ 0 3] +\n",
      "0.0906444041971308 [3^ 0^ 2 1] +\n",
      "0.33173404792821914 [3^ 1^ 1 3] +\n",
      "0.0906444041971308 [3^ 1^ 3 1] +\n",
      "0.09064440419713085 [3^ 2^ 0 1] +\n",
      "0.34869688341114263 [3^ 2^ 2 3] +\n",
      "-0.4759487172683097 [3^ 3] +\n",
      "0.09064440419713085 [3^ 3^ 1 1] +\n",
      "0.34869688341114263 [3^ 3^ 3 3]\n",
      "\n",
      "(-0.09886397351781583+0j) [] +\n",
      "(-0.04532220209856541+0j) [X0 X1 Y2 Y3] +\n",
      "(0.04532220209856541+0j) [X0 Y1 Y2 X3] +\n",
      "(0.04532220209856541+0j) [Y0 X1 X2 Y3] +\n",
      "(-0.04532220209856541+0j) [Y0 Y1 X2 X3] +\n",
      "(0.17119774853325848+0j) [Z0] +\n",
      "(0.16862219143347554+0j) [Z0 Z1] +\n",
      "(0.12054482186554413+0j) [Z0 Z2] +\n",
      "(0.16586702396410954+0j) [Z0 Z3] +\n",
      "(0.1711977485332586+0j) [Z1] +\n",
      "(0.16586702396410954+0j) [Z1 Z2] +\n",
      "(0.12054482186554413+0j) [Z1 Z3] +\n",
      "(-0.22278592890107018+0j) [Z2] +\n",
      "(0.17434844170557132+0j) [Z2 Z3] +\n",
      "(-0.22278592890107013+0j) [Z3]\n"
     ]
    }
   ],
   "source": [
    "from openfermion import get_fermion_operator\n",
    "fermionic_hamiltonian = get_fermion_operator(hamiltonian)\n",
    "print(fermionic_hamiltonian)\n",
    "print()\n",
    "qubit_hamiltonian = jordan_wigner(fermionic_hamiltonian)\n",
    "print(qubit_hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8:** Check that the above operations preserve the eigenspectrum of the H2 Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenspectrum of H2 in FermionOperator: \n",
      " [-1.13727017 -0.53870958 -0.53870958 -0.53247901 -0.53247901 -0.53247901\n",
      " -0.44698572 -0.44698572 -0.16990139  0.23780527  0.23780527  0.35243413\n",
      "  0.35243413  0.47983611  0.71375399  0.92010671]\n"
     ]
    }
   ],
   "source": [
    "print('Eigenspectrum of H2 in FermionOperator: \\n {}'.format(eigenspectrum(fermionic_hamiltonian)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenspectrum of H2 in QubitOperator: \n",
      " [-1.13727017 -0.53870958 -0.53870958 -0.53247901 -0.53247901 -0.53247901\n",
      " -0.44698572 -0.44698572 -0.16990139  0.23780527  0.23780527  0.35243413\n",
      "  0.35243413  0.47983611  0.71375399  0.92010671]\n"
     ]
    }
   ],
   "source": [
    "print('Eigenspectrum of H2 in QubitOperator: \\n {}'.format(eigenspectrum(qubit_hamiltonian)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Are the spectrum equal: True\n"
     ]
    }
   ],
   "source": [
    "print('Are the spectrum equal: {}'.format(numpy.allclose(eigenspectrum(fermionic_hamiltonian),\n",
    "                                                         eigenspectrum(qubit_hamiltonian))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The OpenFermion-Cirq interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the VQE context, the Hamiltonians $H$ we created above above generate cost functions $f(\\vec{\\theta})=\\langle\\Psi(\\vec{\\theta})|H|\\Psi(\\vec{\\theta})\\rangle$ for classical optimization. The openfermioncirq package provides an interface between cirq (for generating and simulating circuits) and openfermion (for generating Hamiltonians) to provide precisely this functionality.\n",
    "\n",
    "As a simple example (before coming back to H$_2$), let us study a trivial VQE on a single qubit - optimize the ansatz $|\\Psi(\\theta)\\rangle=e^{i\\theta Y}|0\\rangle$ for the Hamiltonian $H=Z+X$.\n",
    "\n",
    "In openfermioncirq, a VQE experiment is contained within a openfermioncirq.VariationalStudy class. This requires as input a name, an ansatz, and a cost function.\n",
    "\n",
    "The cost function can be generated from our hamiltonian using the openfermioncirq.HamiltonianObjective class.\n",
    "\n",
    "The variational ansatz itself must extend the openfermioncirq.VariationalAnsatz class (i.e. be a subclass). In doing so, it needs to contain the following methods:\n",
    "\n",
    "- params: a function that returns the parameters to be optimized (in a list/tuple)\n",
    "- operations: a generator for the variational circuit that takes as input a list of qubits (note - a generator, not the final cirq.Circuit - this is typically cleaner)\n",
    "- _generate_qubits: a function that produces a list/tuple containing the qubits in order.\n",
    "(Note, the order of the qubits will correspond to the index order in the openfermion Hamiltonian)\n",
    "\n",
    "**Exercise 9:** use cirq, openfermion, and openfermioncirq to create the above variational ansatz (remember to include a free parameter using sympy!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0): ───Ry(theta)───\n"
     ]
    }
   ],
   "source": [
    "# To describe: VariationalStudy, HamiltonianObjective, OptimizationParams\n",
    "# Toy example - VQE on a single qubit!\n",
    "import cirq\n",
    "import sympy\n",
    "import openfermioncirq\n",
    "\n",
    "theta = sympy.Symbol('theta')\n",
    "\n",
    "class HelloWorldsAnsatz(openfermioncirq.VariationalAnsatz):\n",
    "    \n",
    "    def params(self):\n",
    "        return [theta]\n",
    "\n",
    "    def operations(self, qubits):\n",
    "        yield cirq.ry(theta).on(qubits[0])\n",
    "        \n",
    "    def _generate_qubits(self):\n",
    "        return [cirq.GridQubit(0,0)]\n",
    "\n",
    "ansatz = HelloWorldsAnsatz()\n",
    "\n",
    "hamiltonian = QubitOperator('X0') + QubitOperator('Z0')\n",
    "\n",
    "objective = openfermioncirq.HamiltonianObjective(hamiltonian)\n",
    "study = openfermioncirq.VariationalStudy(\n",
    "    name='Hello, worlds!',\n",
    "    ansatz=ansatz,\n",
    "    objective=objective)\n",
    "\n",
    "print(study.circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To actually simulate (or run) the VQE, we need to define the optimization strategy. Typically for digital quantum simulation, this is performed via gradient descent or gradient-free optimization methods, of which there are many popular choices. All typically require an initial guess --- quite often in a VQE the guess should be slightly off-set to prevent starting in a local maximum or saddle point.\n",
    "\n",
    "Openfermion summarizes these in a neat form in the openfermioncirq.optimization package, which has the OptimizationParams class (basically a container for optimization metaparameters) and various wrappers for calls to scipy routines. For example, let's optimize the above VQE using the COBYLA algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized VQE result: -1.4142135527847444\n",
      "Target Hamiltonian eigenvalues: [-1.41421356  1.41421356]\n"
     ]
    }
   ],
   "source": [
    "from openfermioncirq.optimization import COBYLA, OptimizationParams\n",
    "optimization_params = OptimizationParams(\n",
    "    algorithm=COBYLA,\n",
    "    initial_guess=[0.01])\n",
    "result = study.optimize(optimization_params)\n",
    "print('Optimized VQE result: {}'.format(result.optimal_value))\n",
    "print('Target Hamiltonian eigenvalues: {}'.format(eigenspectrum(hamiltonian)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The unitary coupled cluster ansatz\n",
    "\n",
    "The unitary coupled cluster ansatz is probably the most well-known variational ansatz. We will study this ansatz in more detail later in class, but after going through the above you can hopefully implement it for H$_2$.\n",
    "\n",
    "In short, the unitary coupled cluster ansatz takes the form\n",
    "\n",
    "$$ e^{T(\\vec{\\theta})-T^{\\dagger}(\\vec{\\theta})}|111\\ldots000\\ldots\\rangle $$,\n",
    "\n",
    "where in the starting state the first $\\eta$ spin-orbitals contain an electron and the remaining $N-\\eta$ sites are empty, and the cluster operator $T$ contains terms which excite from empty orbitals to filled. This operator is commonly truncated to single operators $\\theta^i_j\\hat{c}_i^{\\dagger}\\hat{c}_j$ (where $i$ is empty and $j$ is filled) and double operators $\\theta^{i,j}_{k,l}\\hat{c}_i^{\\dagger}\\hat{c}_j^{\\dagger}\\hat{c}_k\\hat{c}_l$. (where $i$ and $j$ are empty and $k$ and $l$ are filled).\n",
    "\n",
    "In the H$_2$ example above (because of the small basis size), we only have two empty and two filled orbitals. Moreover, due to symmetry constraints it turns out that the ground state contains no single excitations, and so the ansatz takes the form of a single term (we can drop the indices on $\\theta$ for brevity)\n",
    "\n",
    "$$ e^{\\theta\\hat{c}^{\\dagger}_3\\hat{c}^{\\dagger}_2\\hat{c}_1\\hat{c}_0-\\mathrm{h.c.}}|1100\\rangle $$\n",
    "\n",
    "**Exercise 9:** convert the cluster operator $T-T^{\\dagger}$ for H$_2$ into qubit form using the Jordan-Wigner transform via openfermion. (Check your result with pen and paper.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster operators:  1.0 [0^ 1^ 2 3] 1.0 [3^ 2^ 1 0]\n"
     ]
    }
   ],
   "source": [
    "# We just create a cluster operator and its Hermitian conjugate with value 1\n",
    "# later we will add a parameter to it.\n",
    "t_fop = FermionOperator('0^ 1^ 2 3')\n",
    "t_fopdag = hermitian_conjugated(t_fop)\n",
    "print('Cluster operators: ', t_fop, t_fopdag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unitary coupled-cluster JW: \n",
      " -0.125j [X0 X1 X2 Y3] +\n",
      "-0.125j [X0 X1 Y2 X3] +\n",
      "0.125j [X0 Y1 X2 X3] +\n",
      "-0.125j [X0 Y1 Y2 Y3] +\n",
      "0.125j [Y0 X1 X2 X3] +\n",
      "-0.125j [Y0 X1 Y2 Y3] +\n",
      "0.125j [Y0 Y1 X2 Y3] +\n",
      "0.125j [Y0 Y1 Y2 X3]\n"
     ]
    }
   ],
   "source": [
    "ucc_qop = jordan_wigner(t_fop - t_fopdag)\n",
    "print('Unitary coupled-cluster JW: \\n', ucc_qop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10:** Write a function to generate the circuit for $\\exp(\\theta c_0^{\\dagger}c_1^{\\dagger}c_2c_3-\\mathrm{h.c.})$ with a free angle $\\theta$. Use this to create a openfermioncirq.VariationalAnsatz.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion import count_qubits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubitop_exponent_to_cirq(operator_list, *qubits, parameter):\n",
    "    '''\n",
    "    Exponentiate a list of QubitOperator to a cirq.\n",
    "    \n",
    "    Exponentiate a list of QubitOperators to yield a cirq circuit generator.\n",
    "    \n",
    "    Args:\n",
    "        operator_list(QubitOperator, list): Operators to transform to a circuit.\n",
    "        qubits(cirq.qubits): List of cirq qubits.\n",
    "    \n",
    "    Yields:\n",
    "        cirq circuit generator\n",
    "    '''\n",
    "    if not isinstance(operator_list, (list, QubitOperator)):\n",
    "        raise TypeError('Input must be a list or QubitOperator')\n",
    "    if isinstance(operator_list, QubitOperator):\n",
    "        operator_list = list(operator_list)\n",
    "        \n",
    "    rot_dic = {'X':lambda q,s: cirq.ry(-1*s*numpy.pi/2).on(q),\n",
    "               'Y':lambda q,s: cirq.rx(-1*s*numpy.pi/2).on(q),\n",
    "               'Z':lambda q,s: cirq.I.on(q)}\n",
    "    \n",
    "    qubit_list=[*qubits]\n",
    "    \n",
    "    for i, operator in enumerate(operator_list):\n",
    "        paulis = list(operator.terms.keys())[0]\n",
    "        sign = numpy.sign(list(operator.terms.values())[0])\n",
    "        \n",
    "        for qbt, pau in paulis:\n",
    "            yield rot_dic[pau](qubit_list[qbt], 1)\n",
    "\n",
    "        if len(paulis) > 1:\n",
    "            for j in range(len(paulis) - 1):\n",
    "                yield cirq.CNOT(qubit_list[paulis[j][0]],\n",
    "                                qubit_list[paulis[j+1][0]])\n",
    "\n",
    "        yield cirq.rz(sign*parameter[0]).on(qubit_list[paulis[-1][0]])\n",
    "        \n",
    "        if len(paulis) > 1:\n",
    "            for j in range(len(paulis)-1, 0, -1):\n",
    "                yield cirq.CNOT(qubit_list[paulis[j-1][0]],\n",
    "                                qubit_list[paulis[j][0]])\n",
    "        \n",
    "        for qbt, pau in paulis:\n",
    "            yield rot_dic[pau](qubit_list[qbt], -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "class UCCAnstaz(openfermioncirq.VariationalAnsatz):\n",
    "    def __init__(self, operators, num_qubits, num_params):\n",
    "        '''\n",
    "        Start variational ansatz class.\n",
    "        \n",
    "        Args:\n",
    "            operators(QubitOperator): UCCSD Pauli operators.\n",
    "            num_qubits(int): Number of qubits to create the circuit.\n",
    "            num_params(int): Number of parameters.\n",
    "            \n",
    "        Note:\n",
    "            If one wants to be accurate when using CCSD amplitudes\n",
    "            as parameters for UCCSD, the number of free parameters\n",
    "            is much smaller than the number of Pauli operators.\n",
    "            Hence one needs to specify how many parameters are allowed.\n",
    "            Beawere that parameters are added by module num_params.\n",
    "        '''\n",
    "        self.operators = operators\n",
    "        self.num_qubits = num_qubits\n",
    "        self.num_params = num_params\n",
    "        super().__init__(qubits = cirq.LineQubit.range(num_qubits))\n",
    "    \n",
    "    def params(self):\n",
    "        return [sympy.Symbol('theta_'+str(i))\n",
    "                for i in range(self.num_params)]\n",
    "        \n",
    "    def _generate_qubits(self):\n",
    "        pass\n",
    "    \n",
    "    def operations(self, qubits):\n",
    "        yield qubitop_exponent_to_cirq(self.operators, *qubits,\n",
    "                                       parameter=self.params())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz=UCCAnstaz(list(ucc_qop), count_qubits(ucc_qop), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"overflow: auto; white-space: pre;\">0: ───Rx(-0.5π)───@─────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@────────────────────────────────────────────@──────────Ry(0.5π)───\n",
       "                  │                                             │                                   │                                            │                                   │                                            │                                   │                                            │                                   │                                             │                                   │                                             │                                   │                                             │                                   │                                            │\n",
       "1: ───Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@─────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@──────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X──────────Rx(0.5π)───\n",
       "                      │                              │                                                  │                             │                                                  │                             │                                                  │                             │                                                  │                              │                                                  │                              │                                                  │                              │                                                  │                             │\n",
       "2: ───Rx(-0.5π)───────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@─────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@──────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)──────────────\n",
       "                          │                      │                                                          │                     │                                                          │                     │                                                          │                     │                                                          │                      │                                                          │                      │                                                          │                      │                                                          │                     │\n",
       "3: ───Rx(-0.5π)───────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)─────────────────────────</pre>"
      ],
      "text/plain": [
       "0: ───Rx(-0.5π)───@─────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@────────────────────────────────────────────@──────────Ry(0.5π)───\n",
       "                  │                                             │                                   │                                            │                                   │                                            │                                   │                                            │                                   │                                             │                                   │                                             │                                   │                                             │                                   │                                            │\n",
       "1: ───Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@─────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@──────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X──────────Rx(0.5π)───\n",
       "                      │                              │                                                  │                             │                                                  │                             │                                                  │                             │                                                  │                              │                                                  │                              │                                                  │                              │                                                  │                             │\n",
       "2: ───Rx(-0.5π)───────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@─────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@──────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)──────────────\n",
       "                          │                      │                                                          │                     │                                                          │                     │                                                          │                     │                                                          │                      │                                                          │                      │                                                          │                      │                                                          │                     │\n",
       "3: ───Rx(-0.5π)───────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)─────────────────────────"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ansatz.circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11:** Write a function to prepare the computational basis state $c_1^{\\dagger}c_0^{\\dagger}|0000\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_state_prepartion(qubits, positions):\n",
    "    return cirq.Circuit(cirq.X(qubits[i]) for i in positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"overflow: auto; white-space: pre;\">0: ───X───\n",
       "\n",
       "1: ───X───</pre>"
      ],
      "text/plain": [
       "0: ───X───\n",
       "\n",
       "1: ───X───"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initial_state_prepartion(ansatz.qubits, positions=[0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 12:** Create a study of the H$_2$ molecule using OpenFermion-Cirq, optimize the resulting VQE, and compare your results to the Full-CI energy reported above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = openfermioncirq.HamiltonianObjective(qubit_hamiltonian)\n",
    "study = openfermioncirq.VariationalStudy(\n",
    "    name='UCC for H2',\n",
    "    ansatz=ansatz,\n",
    "    preparation_circuit=initial_state_prepartion(ansatz.qubits, [0,1]),\n",
    "    objective=objective)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"overflow: auto; white-space: pre;\">0: ───X───Rx(-0.5π)───@─────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@────────────────────────────────────────────@──────────Ry(0.5π)───\n",
       "                      │                                             │                                   │                                            │                                   │                                            │                                   │                                            │                                   │                                             │                                   │                                             │                                   │                                             │                                   │                                            │\n",
       "1: ───X───Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@─────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@──────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X──────────Rx(0.5π)───\n",
       "                          │                              │                                                  │                             │                                                  │                             │                                                  │                             │                                                  │                              │                                                  │                              │                                                  │                              │                                                  │                             │\n",
       "2: ───────Rx(-0.5π)───────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@─────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@──────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)──────────────\n",
       "                              │                      │                                                          │                     │                                                          │                     │                                                          │                     │                                                          │                      │                                                          │                      │                                                          │                      │                                                          │                     │\n",
       "3: ───────Rx(-0.5π)───────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)─────────────────────────</pre>"
      ],
      "text/plain": [
       "0: ───X───Rx(-0.5π)───@─────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Rx(-0.5π)───@────────────────────────────────────────────@───────────Rx(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@─────────────────────────────────────────────@───────────Ry(0.5π)────Ry(-0.5π)───@────────────────────────────────────────────@──────────Ry(0.5π)───\n",
       "                      │                                             │                                   │                                            │                                   │                                            │                                   │                                            │                                   │                                             │                                   │                                             │                                   │                                             │                                   │                                            │\n",
       "1: ───X───Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@─────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X───────────Rx(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Ry(-0.5π)───X───@──────────────────────────────@──────────X───────────Ry(0.5π)────Rx(-0.5π)───X───@──────────────────────────────@──────────X───────────Rx(0.5π)────Rx(-0.5π)───X───@─────────────────────────────@──────────X──────────Rx(0.5π)───\n",
       "                          │                              │                                                  │                             │                                                  │                             │                                                  │                             │                                                  │                              │                                                  │                              │                                                  │                              │                                                  │                             │\n",
       "2: ───────Rx(-0.5π)───────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@─────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@──────────────────────@───X──────────Ry(0.5π)────Rx(-0.5π)───────────────────X───@──────────────────────@───X──────────Rx(0.5π)────Ry(-0.5π)───────────────────X───@─────────────────────@───X──────────Ry(0.5π)──────────────\n",
       "                              │                      │                                                          │                     │                                                          │                     │                                                          │                     │                                                          │                      │                                                          │                      │                                                          │                      │                                                          │                     │\n",
       "3: ───────Rx(-0.5π)───────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Ry(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Rx(-0.5π)───────────────────────────────────X───Rz(-1.0*theta_0)───X───Rx(0.5π)───Ry(-0.5π)───────────────────────────────────X───Rz(1.0*theta_0)───X───Ry(0.5π)─────────────────────────"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study.circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized VQE result: -1.1372700216580982\n"
     ]
    }
   ],
   "source": [
    "optimization_params = OptimizationParams(\n",
    "    algorithm=COBYLA,\n",
    "    initial_guess=numpy.random.rand(ansatz.num_params))\n",
    "result = study.optimize(optimization_params)\n",
    "print('Optimized VQE result: {}'.format(result.optimal_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full configuration-interaction energy: -1.137270174625328\n"
     ]
    }
   ],
   "source": [
    "print('Full configuration-interaction energy: {}'.format(molecule.fci_energy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized parameters: \n",
      " [0.75705424]\n"
     ]
    }
   ],
   "source": [
    "print('Optimized parameters: \\n',result.optimal_parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Repeat the same study for different bond distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13:** Openfermion also contains data for the H$_2$ molecule at diatomic bond lengths from 0.3 to 2.5 Angstrom in increments of 0.1 Angstrom. Repeat the above study for all these curves, and plot the bond dissociation curve (a plot of the ground state energy as a function of the bond length) for the H$_2$ molecule. Compare to the FCI results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 14:** The UCC circuit that you generated is horribly long - how much can you optimize it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is left as an exercise for the students"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}