{
 "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here!"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here!"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here!"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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.\n",
    "\n",
    "**Exercise 6:** Recall the conditions for a QubitOperator to be Hermitian, and check that this is preserved under the Jordan-Wigner transformation.\n",
    "\n",
    "**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": "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": null,
   "metadata": {},
   "outputs": [],
   "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)"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "        ### Insert your code here.\n",
    "\n",
    "    def operations(self, qubits):\n",
    "        ### Insert your code here\n",
    "        \n",
    "    def _generate_qubits(self):\n",
    "        ### Insert your code here\n",
    "\n",
    "ansatz = HelloWorldsAnsatz()\n",
    "\n",
    "### Insert your code here.\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": null,
   "metadata": {},
   "outputs": [],
   "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.)\n",
    " \n",
    "**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",
    "\n",
    "**Exercise 11:** Write a function to prepare the computational basis state $c_1^{\\dagger}c_0^{\\dagger}|0000\\rangle$.\n",
    "\n",
    "**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.\n",
    "\n",
    "**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.\n",
    "\n",
    "**Exercise 14:** The UCC circuit that you generated is horribly long - how much can you optimize it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Insert your code here!"
   ]
  }
 ],
 "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.5rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}