OpenFermion hands-on tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OpenFermion hands-on tutorial\n",
    "\n",
    "OpenFermion is a library for obtaining and manipulating representations of fermionic and qubit Hamiltonians for the quantum simulation of chemistry and materials. This tutorial introduces the core data structures of OpenFermion and some basic functionality.\n",
    "\n",
    "## Background\n",
    "\n",
    "A system of $N$ fermionic modes is\n",
    "described by a set of fermionic *annihilation operators*\n",
    "$\\{a_p\\}_{p=0}^{N-1}$ satisfying the *canonical anticommutation relations*\n",
    "$$\\begin{aligned}\n",
    "    \\{a_p, a_q\\} &= 0, \\\\\n",
    "    \\{a_p, a^\\dagger_q\\} &= \\delta_{pq},\n",
    "  \\end{aligned}$$ where $\\{A, B\\} := AB + BA$. The adjoint\n",
    "$a^\\dagger_p$ of an annihilation operator $a_p$ is called a *creation\n",
    "operator*, and we refer to creation and annihilation operators as\n",
    "fermionic *ladder operators*.\n",
    "    \n",
    "The canonical anticommutation relations impose a number of consequences on the structure of the vector space on which the ladder operators act; see [Michael Nielsen's notes](http://michaelnielsen.org/blog/archive/notes/fermions_and_jordan_wigner.pdf) for a good discussion.\n",
    "\n",
    "## FermionOperator and QubitOperator\n",
    "\n",
    "### FermionOperator\n",
    "\n",
    "- Stores a weighted sum (linear combination) of fermionic terms\n",
    "- A fermionic term is a product of ladder operators\n",
    "- Examples of FermionOperators:\n",
    "$$\n",
    "\\begin{align}\n",
    "& a_1 \\nonumber \\\\\n",
    "& 1.7 a^\\dagger_3 \\nonumber \\\\\n",
    "&-1.7 \\, a^\\dagger_3 a_1 \\nonumber \\\\\n",
    "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 \\nonumber \\\\\n",
    "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 - 1.7 \\, a^\\dagger_3 a_1 \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "- A fermionic term is internally represented as a tuple of tuples\n",
    "- Each inner tuple represents a single ladder operator as (index, action)\n",
    "- Examples of fermionic terms:\n",
    "$$\n",
    "\\begin{align}\n",
    "I & \\mapsto () \\nonumber \\\\\n",
    "a_1 & \\mapsto ((1, 0),) \\nonumber \\\\\n",
    "a^\\dagger_3 & \\mapsto ((3, 1),) \\nonumber \\\\\n",
    "a^\\dagger_3 a_1 & \\mapsto ((3, 1), (1, 0)) \\nonumber \\\\\n",
    "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto ((4, 1), (3, 1), (9, 0), (1, 0)) \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "- FermionOperator is a sum of terms, represented as a dictionary from term to coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{((4, 1), (3, 1), (9, 0), (1, 0)): (1+2j), ((3, 1), (1, 0)): -1.7}\n"
     ]
    }
   ],
   "source": [
    "import openfermion as of\n",
    "\n",
    "op = of.FermionOperator(((4, 1), (3, 1), (9, 0), (1, 0)), 1+2j) + of.FermionOperator(((3, 1), (1, 0)), -1.7)\n",
    "\n",
    "print(op.terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative notation, useful when playing around:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "I & \\mapsto \\textrm{\"\"} \\nonumber \\\\\n",
    "a_1 & \\mapsto \\textrm{\"1\"} \\nonumber \\\\\n",
    "a^\\dagger_3 & \\mapsto \\textrm{\"3^\"} \\nonumber \\\\\n",
    "a^\\dagger_3 a_1 & \\mapsto \\textrm{\"3^}\\;\\textrm{1\"} \\nonumber \\\\\n",
    "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto \\textrm{\"4^}\\;\\textrm{3^}\\;\\textrm{9}\\;\\textrm{1\"} \\nonumber\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{((4, 1), (3, 1), (9, 0), (1, 0)): (1+2j), ((3, 1), (1, 0)): -1.7}\n"
     ]
    }
   ],
   "source": [
    "op = of.FermionOperator('4^ 3^ 9 1', 1+2j) + of.FermionOperator('3^ 1', -1.7)\n",
    "\n",
    "print(op.terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just print the operator for a nice readable representation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.7 [3^ 1] +\n",
      "(1+2j) [4^ 3^ 9 1]\n"
     ]
    }
   ],
   "source": [
    "print(op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `normal_ordered` function to normal-order FermionOperators. Normal-ordering a term means using the anticommutation relations to write the term into an equivalent form where all creation operators precede annihilation operators and higher indices precede lower indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.0 [3^ 2^ 1 0]\n",
      "\n",
      "1.0 [] +\n",
      "-1.0 [0^ 0]\n"
     ]
    }
   ],
   "source": [
    "op = of.FermionOperator('0 1 3^ 2^')\n",
    "\n",
    "print(of.normal_ordered(op))\n",
    "print()\n",
    "\n",
    "op = of.FermionOperator('0 0^')\n",
    "print(of.normal_ordered(op))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "- Use the `*`, `+`, and `**` operators and the `normal_ordered` function to verify the following:\n",
    "$$\n",
    "\\begin{align}\n",
    "    (a_2)^2 &= 0 \\\\\n",
    "    a_2 a_7 + a_7 a_2 &= 0 \\\\\n",
    "    a_2 a_7^\\dagger + a_7^\\dagger a_2 &= 0\\\\\n",
    "    a_2 a_2^\\dagger + a_2^\\dagger a_2 &= 1\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Use the `==` operator and the `normal_ordered` function to verify that\n",
    "$$\n",
    "a_1^\\dagger a_2 a_2^\\dagger a_4 = a_1^\\dagger a_4 - a_1^\\dagger a_2^\\dagger a_2 a_4\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QubitOperator\n",
    "\n",
    "Same as FermionOperator, but the possible actions are 'X', 'Y', and 'Z' instead of 1 and 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 [X1 Y2 Z3] +\n",
      "3.0 [X3 Z4]\n"
     ]
    }
   ],
   "source": [
    "op = of.QubitOperator(((1, 'X'), (2, 'Y'), (3, 'Z')))\n",
    "op += of.QubitOperator('X3 Z4', 3.0)\n",
    "\n",
    "print(op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FermionOperator and QubitOperator actually inherit from the same parent class, SymbolicOperator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "- Use the `commutator` function to verify the following:\n",
    "$$\n",
    "\\begin{align}\n",
    "[X_0, Y_0] &= 2 i Z_0 \\\\\n",
    "[Z_0, Y_0] &= -2 i X_0 \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Jordan-Wigner and Bravyi-Kitaev transforms\n",
    "\n",
    "A fermionic transform maps FermionOperators to QubitOperators in a way that preserves the canonical anticommutation relations. The most basic transforms are the Jordan-Wigner transform (JWT) and Bravyi-Kitaev transform (BKT). Note that the BKT requires the total number of qubits to be predetermined. Whenever a fermionic transform is being applied implicitly, it is the JWT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.25+0j) [X2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 Z12 Z13 Z14 X15] +\n",
      "0.25j [X2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 Z12 Z13 Z14 Y15] +\n",
      "-0.25j [Y2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 Z12 Z13 Z14 X15] +\n",
      "(0.25+0j) [Y2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 Z12 Z13 Z14 Y15]\n",
      "\n",
      "(-0.25+0j) [Z1 X2 X3 X7 Z15] +\n",
      "-0.25j [Z1 X2 X3 Y7 Z11 Z13 Z14] +\n",
      "0.25j [Z1 Y2 X3 X7 Z15] +\n",
      "(-0.25+0j) [Z1 Y2 X3 Y7 Z11 Z13 Z14]\n"
     ]
    }
   ],
   "source": [
    "op = of.FermionOperator('2^ 15')\n",
    "\n",
    "print(of.jordan_wigner(op))\n",
    "print()\n",
    "print(of.bravyi_kitaev(op, n_qubits=16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "- Below are some examples of how FermionOperators are mapped to QubitOperators by the Jordan-Wigner transform (the notation 'h.c.' stands for 'hermitian conjugate'):\n",
    "$$\n",
    "\\begin{align*}\n",
    "    a_p^\\dagger &\\mapsto \\frac12 (X_p - i Y_p) Z_0 \\cdots Z_{p-1}\\\\\n",
    "    a_p^\\dagger a_p &\\mapsto \\frac12 (I - Z_p)\\\\\n",
    "    (\\beta a_p^\\dagger a_q + \\text{h.c.}) &\\mapsto \\frac12 [\\text{Re}(\\beta) (X_p ZZ \\cdots ZZ X_q + Y_p ZZ \\cdots ZZ Y_q) + \\text{Im}(\\beta) (Y_p ZZ \\cdots ZZ X_q - X_p ZZ \\cdots ZZ Y_q)]\n",
    "\\end{align*}\n",
    "$$\n",
    "Verify these mappings for $p=2$ and $q=7$. The `hermitian_conjugated` function may be useful here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Verify that after applying the JWT to ladder operators, the resulting QubitOperators still satisfy\n",
    "$$\n",
    "\\begin{align}\n",
    "    (a_2)^2 &= 0 \\\\\n",
    "    a_2 a_7 + a_7 a_2 &= 0 \\\\\n",
    "    a_2 a_7^\\dagger + a_7^\\dagger a_2 &= 0\\\\\n",
    "    a_2 a_2^\\dagger + a_2^\\dagger a_2 &= 1\n",
    "\\end{align}\n",
    "$$\n",
    "You can use the `anticommutator` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sparse matrices for numerical computations\n",
    "\n",
    "- QubitOperators can be converted to Scipy sparse matrices using `get_sparse_operator`\n",
    "- Sparse matrices can be used to obtain eigenvalues, expectation values, etc.\n",
    "- As a demonstration, we will compute the ground state energy of a Hubbard model Hamiltonian:\n",
    "$$\n",
    "H = - t \\sum_{\\langle i,j \\rangle} \\sum_{\\sigma}\n",
    "             (a^\\dagger_{i, \\sigma} a_{j, \\sigma} +\n",
    "              a^\\dagger_{j, \\sigma} a_{i, \\sigma})\n",
    "     + U \\sum_{i} a^\\dagger_{i, \\uparrow} a_{i, \\uparrow}\n",
    "                  a^\\dagger_{j, \\downarrow} a_{j, \\downarrow}\n",
    "$$\n",
    "where $\\langle i, j \\rangle$ ranges over the edges in a 2-dimensional grid and $\\sigma$ is either spin up or spin down.\n",
    "- The Hubbard model is an idealized Hamiltonian model that may capture qualitative aspects of high-temperature superconductors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.0 [0^ 0 1^ 1] +\n",
      "-1.0 [0^ 2] +\n",
      "-1.0 [0^ 4] +\n",
      "-1.0 [1^ 3] +\n",
      "-1.0 [1^ 5] +\n",
      "-1.0 [2^ 0] +\n",
      "4.0 [2^ 2 3^ 3] +\n",
      "-1.0 [2^ 6] +\n",
      "-1.0 [3^ 1] +\n",
      "-1.0 [3^ 7] +\n",
      "-1.0 [4^ 0] +\n",
      "4.0 [4^ 4 5^ 5] +\n",
      "-1.0 [4^ 6] +\n",
      "-1.0 [5^ 1] +\n",
      "-1.0 [5^ 7] +\n",
      "-1.0 [6^ 2] +\n",
      "-1.0 [6^ 4] +\n",
      "4.0 [6^ 6 7^ 7] +\n",
      "-1.0 [7^ 3] +\n",
      "-1.0 [7^ 5]\n"
     ]
    }
   ],
   "source": [
    "# Set Hubbard model parameters\n",
    "x_dim = 2\n",
    "y_dim = 2\n",
    "tunneling = 1.0\n",
    "coulomb = 4.0\n",
    "\n",
    "# Create Hubbard model Hamiltonian as FermionOperator\n",
    "hubbard_model = of.fermi_hubbard(x_dim, y_dim, tunneling, coulomb)\n",
    "\n",
    "print(hubbard_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4+0j) [] +\n",
      "(-0.5+0j) [X0 Z1 X2] +\n",
      "(-0.5+0j) [X0 Z1 Z2 Z3 X4] +\n",
      "(-0.5+0j) [Y0 Z1 Y2] +\n",
      "(-0.5+0j) [Y0 Z1 Z2 Z3 Y4] +\n",
      "(-1+0j) [Z0] +\n",
      "(1+0j) [Z0 Z1] +\n",
      "(-0.5+0j) [X1 Z2 X3] +\n",
      "(-0.5+0j) [X1 Z2 Z3 Z4 X5] +\n",
      "(-0.5+0j) [Y1 Z2 Y3] +\n",
      "(-0.5+0j) [Y1 Z2 Z3 Z4 Y5] +\n",
      "(-1+0j) [Z1] +\n",
      "(-0.5+0j) [X2 Z3 Z4 Z5 X6] +\n",
      "(-0.5+0j) [Y2 Z3 Z4 Z5 Y6] +\n",
      "(-1+0j) [Z2] +\n",
      "(1+0j) [Z2 Z3] +\n",
      "(-0.5+0j) [X3 Z4 Z5 Z6 X7] +\n",
      "(-0.5+0j) [Y3 Z4 Z5 Z6 Y7] +\n",
      "(-1+0j) [Z3] +\n",
      "(-0.5+0j) [X4 Z5 X6] +\n",
      "(-0.5+0j) [Y4 Z5 Y6] +\n",
      "(-1+0j) [Z4] +\n",
      "(1+0j) [Z4 Z5] +\n",
      "(-0.5+0j) [X5 Z6 X7] +\n",
      "(-0.5+0j) [Y5 Z6 Y7] +\n",
      "(-1+0j) [Z5] +\n",
      "(-1+0j) [Z6] +\n",
      "(1+0j) [Z6 Z7] +\n",
      "(-1+0j) [Z7]\n"
     ]
    }
   ],
   "source": [
    "# Map to QubitOperator using the JWT\n",
    "hubbard_jw = of.jordan_wigner(hubbard_model)\n",
    "\n",
    "print(hubbard_jw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3.418550718873853\n"
     ]
    }
   ],
   "source": [
    "# Convert to Scipy sparse matrix\n",
    "hubbard_jw_sparse = of.get_sparse_operator(hubbard_jw)\n",
    "\n",
    "# Compute ground state energy\n",
    "ground_energy_jw, _ = of.get_ground_state(hubbard_jw_sparse)\n",
    "\n",
    "print(ground_energy_jw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "- Compute the ground energy of the same Hamiltonian, but via the Bravyi-Kitaev transform. Verify that you get the same value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array data structures\n",
    "\n",
    "- When FermionOperators have specialized structure we can store coefficients in numpy arrays, enabling easier numerical manipulation.\n",
    "- Array data structures can always be converted to FermionOperator using `get_fermion_operator`.\n",
    "\n",
    "### InteractionOperator\n",
    "\n",
    "- Stores the one- and two-body tensors $T_{pq}$ and $V_{pqrs}$ of the molecular Hamiltonian\n",
    "$$\n",
    "\\sum_{pq} T_{pq} a_p^\\dagger a_q + \\sum_{pqrs} V_{pqrs} a_p^\\dagger a_q^\\dagger a_r a_s\n",
    "$$\n",
    "- Default data structure for molecular Hamiltonians\n",
    "- Convert from FermionOperator using `get_interaction_operator`\n",
    "\n",
    "### DiagonalCoulombHamiltonian\n",
    "\n",
    "- Stores the one- and two-body coefficient matrices $T_{pq}$ and $V_{pq}$ of a Hamiltonian with a diagonal Coulomb term:\n",
    "$$\n",
    "\\sum_{pq} T_{pq} a_p^\\dagger a_q + \\sum_{pq} V_{pq} a_p^\\dagger a_p a_q^\\dagger a_q\n",
    "$$\n",
    "- Leads to especially efficient algorithms for quantum simulation\n",
    "- Convert from FermionOperator using `get_diagonal_coulomb_hamiltonian`\n",
    "\n",
    "### QuadraticHamiltonian\n",
    "\n",
    "- Stores the Hermitian matrix $M_{pq}$ and antisymmetric matrix $\\Delta_{pq}$ describing a general quadratic Hamiltonian\n",
    "$$\n",
    "\\sum_{p, q} M_{pq} a^\\dagger_p a_q\n",
    "+ \\frac12 \\sum_{p, q}\n",
    "    (\\Delta_{pq} a^\\dagger_p a^\\dagger_q + \\text{h.c.})\n",
    "$$\n",
    "- Routines included for efficient diagonalization (can handle thousands of fermionic modes)\n",
    "- Convert from FermionOperator using `get_quadratic_hamiltonian`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hubbard model has the form of a DiagonalCoulombHamiltonian, so it can be converted to one. OpenFermion-Cirq contains specialized low-depth Hamiltonian simulation algorithms which take this data structure as input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.+0.j  0.+0.j -1.+0.j  0.+0.j -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  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  0.+0.j  0.+0.j -1.+0.j  0.+0.j]\n",
      " [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j]\n",
      " [-1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j]\n",
      " [ 0.+0.j -1.+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 -1.+0.j  0.+0.j -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  0.+0.j -1.+0.j  0.+0.j  0.+0.j]]\n",
      "\n",
      "[[0. 2. 0. 0. 0. 0. 0. 0.]\n",
      " [2. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 2. 0. 0. 0. 0.]\n",
      " [0. 0. 2. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 2. 0. 0.]\n",
      " [0. 0. 0. 0. 2. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 2.]\n",
      " [0. 0. 0. 0. 0. 0. 2. 0.]]\n",
      "\n",
      "4.0 [0^ 0 1^ 1] +\n",
      "-1.0 [0^ 2] +\n",
      "-1.0 [0^ 4] +\n",
      "-1.0 [1^ 3] +\n",
      "-1.0 [1^ 5] +\n",
      "-1.0 [2^ 0] +\n",
      "4.0 [2^ 2 3^ 3] +\n",
      "-1.0 [2^ 6] +\n",
      "-1.0 [3^ 1] +\n",
      "-1.0 [3^ 7] +\n",
      "-1.0 [4^ 0] +\n",
      "4.0 [4^ 4 5^ 5] +\n",
      "-1.0 [4^ 6] +\n",
      "-1.0 [5^ 1] +\n",
      "-1.0 [5^ 7] +\n",
      "-1.0 [6^ 2] +\n",
      "-1.0 [6^ 4] +\n",
      "4.0 [6^ 6 7^ 7] +\n",
      "-1.0 [7^ 3] +\n",
      "-1.0 [7^ 5]\n"
     ]
    }
   ],
   "source": [
    "hubbard_hamiltonian = of.get_diagonal_coulomb_hamiltonian(hubbard_model)\n",
    "\n",
    "print(hubbard_hamiltonian.one_body)\n",
    "print()\n",
    "print(hubbard_hamiltonian.two_body)\n",
    "print()\n",
    "print(hubbard_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below demonstrates using one of our electronic structure package plugins, OpenFermion-PySCF, to generate a molecular Hamiltonian for a hydrogen molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.66147151365 [] +\n",
      "-1.2178260299951058 [0^ 0] +\n",
      "0.3316650744318082 [0^ 0^ 0 0] +\n",
      "0.09231339177803066 [0^ 0^ 2 2] +\n",
      "0.3316650744318082 [0^ 1^ 1 0] +\n",
      "0.09231339177803066 [0^ 1^ 3 2] +\n",
      "0.09231339177803066 [0^ 2^ 0 2] +\n",
      "0.3267206861819477 [0^ 2^ 2 0] +\n",
      "0.09231339177803066 [0^ 3^ 1 2] +\n",
      "0.3267206861819477 [0^ 3^ 3 0] +\n",
      "0.3316650744318082 [1^ 0^ 0 1] +\n",
      "0.09231339177803066 [1^ 0^ 2 3] +\n",
      "-1.2178260299951058 [1^ 1] +\n",
      "0.3316650744318082 [1^ 1^ 1 1] +\n",
      "0.09231339177803066 [1^ 1^ 3 3] +\n",
      "0.09231339177803066 [1^ 2^ 0 3] +\n",
      "0.3267206861819477 [1^ 2^ 2 1] +\n",
      "0.09231339177803066 [1^ 3^ 1 3] +\n",
      "0.3267206861819477 [1^ 3^ 3 1] +\n",
      "0.32672068618194783 [2^ 0^ 0 2] +\n",
      "0.09231339177803066 [2^ 0^ 2 0] +\n",
      "0.32672068618194783 [2^ 1^ 1 2] +\n",
      "0.09231339177803066 [2^ 1^ 3 0] +\n",
      "-0.5096378744364826 [2^ 2] +\n",
      "0.09231339177803066 [2^ 2^ 0 0] +\n",
      "0.34339576784573445 [2^ 2^ 2 2] +\n",
      "0.09231339177803066 [2^ 3^ 1 0] +\n",
      "0.34339576784573445 [2^ 3^ 3 2] +\n",
      "0.32672068618194783 [3^ 0^ 0 3] +\n",
      "0.09231339177803066 [3^ 0^ 2 1] +\n",
      "0.32672068618194783 [3^ 1^ 1 3] +\n",
      "0.09231339177803066 [3^ 1^ 3 1] +\n",
      "0.09231339177803066 [3^ 2^ 0 1] +\n",
      "0.34339576784573445 [3^ 2^ 2 3] +\n",
      "-0.5096378744364826 [3^ 3] +\n",
      "0.09231339177803066 [3^ 3^ 1 1] +\n",
      "0.34339576784573445 [3^ 3^ 3 3]\n"
     ]
    }
   ],
   "source": [
    "import openfermionpyscf as ofpyscf\n",
    "\n",
    "# Set molecule parameters\n",
    "geometry = [('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.8))]\n",
    "basis = 'sto-3g'\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "\n",
    "# Perform electronic structure calculations and\n",
    "# obtain Hamiltonian as an InteractionOperator\n",
    "hamiltonian = ofpyscf.generate_molecular_hamiltonian(\n",
    "    geometry, basis, multiplicity, charge)\n",
    "\n",
    "# Convert to a FermionOperator\n",
    "hamiltonian_ferm_op = of.get_fermion_operator(hamiltonian)\n",
    "\n",
    "print(hamiltonian_ferm_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "- Create a FermionOperator that is quadratic and Hermitian and use `get_quadratic_hamiltonian` to convert it to a QuadraticHamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "- The BCS mean-field d-wave model of superconductivity has the Hamiltonian\n",
    "$$\n",
    "H = - t \\sum_{\\langle i,j \\rangle} \\sum_\\sigma\n",
    "        (a^\\dagger_{i, \\sigma} a_{j, \\sigma} +\n",
    "         a^\\dagger_{j, \\sigma} a_{i, \\sigma})\n",
    "    - \\sum_{\\langle i,j \\rangle} \\Delta_{ij}\n",
    "      (a^\\dagger_{i, \\uparrow} a^\\dagger_{j, \\downarrow} -\n",
    "       a^\\dagger_{i, \\downarrow} a^\\dagger_{j, \\uparrow} +\n",
    "       a_{j, \\downarrow} a_{i, \\uparrow} -\n",
    "       a_{j, \\uparrow} a_{i, \\downarrow})\n",
    "$$\n",
    "Use the `mean_field_dwave` function to generate an instance of this model with dimensions 10x10 (Tip: Type `of.mean_field_dwave` and then press Shift-Tab. This will reveal the function arguments.)\n",
    "    - Convert the Hamiltonian to a QubitOperator with the JWT. What is the length of the longest Pauli string that appears?\n",
    "    - Convert the Hamiltonian to a QubitOperator with the BKT. What is the length of the longest Pauli string that appears?\n",
    "    - Convert the Hamiltonian to a QuadraticHamiltonian. Get its ground energy using the `ground_energy` method of QuadraticHamiltonian. What would happen if you tried to compute the ground energy by converting to a sparse matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code cell for exercises\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions to exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n",
      "1.0 []\n",
      "\n",
      "True\n",
      "\n",
      "2j [Z0]\n",
      "-2j [X0]\n",
      "\n",
      "0.5 [Z0 Z1 X2] +\n",
      "-0.5j [Z0 Z1 Y2]\n",
      "\n",
      "(0.5+0j) [] +\n",
      "(-0.5+0j) [Z2]\n",
      "\n",
      "(1+0j) [X2 Z3 Z4 Z5 Z6 X7] +\n",
      "(-1.5+0j) [X2 Z3 Z4 Z5 Z6 Y7] +\n",
      "(1.5+0j) [Y2 Z3 Z4 Z5 Z6 X7] +\n",
      "(1+0j) [Y2 Z3 Z4 Z5 Z6 Y7]\n",
      "\n",
      "0\n",
      "0\n",
      "(1+0j) []\n",
      "\n",
      "-3.4185507188738278\n",
      "\n",
      "182\n",
      "15\n",
      "-206.99159932863176\n"
     ]
    }
   ],
   "source": [
    "# Solutions to exercises\n",
    "\n",
    "a2 = of.FermionOperator('2')\n",
    "a2dag = of.FermionOperator('2^')\n",
    "a7 = of.FermionOperator('7')\n",
    "a7dag = of.FermionOperator('7^')\n",
    "\n",
    "print(of.normal_ordered(a2**2))\n",
    "print(of.normal_ordered(a2*a7 + a7*a2))\n",
    "print(of.normal_ordered(a2*a7dag + a7dag*a2))\n",
    "print(of.normal_ordered(a2*a2dag + a2dag*a2))\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "op1 = of.FermionOperator('1^ 2 2^ 4')\n",
    "op2 = of.FermionOperator('1^ 4') - of.FermionOperator('1^ 2^ 2 4')\n",
    "\n",
    "print(of.normal_ordered(op1) == of.normal_ordered(op2))\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "X = of.QubitOperator('X0')\n",
    "Y = of.QubitOperator('Y0')\n",
    "Z = of.QubitOperator('Z0')\n",
    "\n",
    "print(of.commutator(X, Y))\n",
    "print(of.commutator(Z, Y))\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "print(of.jordan_wigner(a2dag))\n",
    "print()\n",
    "print(of.jordan_wigner(a2dag*a2))\n",
    "print()\n",
    "\n",
    "op = (2+3j)*a2dag*a7\n",
    "op += of.hermitian_conjugated(op)\n",
    "print(of.jordan_wigner(op))\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "a2_jw = of.jordan_wigner(a2)\n",
    "a2dag_jw = of.jordan_wigner(a2dag)\n",
    "a7_jw = of.jordan_wigner(a7)\n",
    "a7dag_jw = of.jordan_wigner(a7dag)\n",
    "\n",
    "print(of.anticommutator(a2_jw, a7_jw))\n",
    "print(of.anticommutator(a2_jw, a7dag_jw))\n",
    "print(of.anticommutator(a2_jw, a2dag_jw))\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "hubbard_bk = of.bravyi_kitaev(hubbard_model)\n",
    "hubbard_bk_sparse = of.get_sparse_operator(hubbard_bk)\n",
    "ground_energy_bk, _ = of.get_ground_state(hubbard_bk_sparse)\n",
    "\n",
    "print(ground_energy_bk)\n",
    "print()\n",
    "\n",
    "# ----\n",
    "\n",
    "op = of.FermionOperator('0^ 1') + of.FermionOperator('1^ 0')\n",
    "quad_ham = of.get_quadratic_hamiltonian(op)\n",
    "\n",
    "# ----\n",
    "\n",
    "mean_field_model = of.mean_field_dwave(10, 10, 1.0, 1.0)\n",
    "\n",
    "mean_field_jw = of.jordan_wigner(mean_field_model)\n",
    "print(max(len(term) for term in mean_field_jw.terms))\n",
    "\n",
    "mean_field_bk = of.bravyi_kitaev(mean_field_model)\n",
    "print(max(len(term) for term in mean_field_bk.terms))\n",
    "\n",
    "mean_field_quad_ham = of.get_quadratic_hamiltonian(mean_field_model)\n",
    "print(mean_field_quad_ham.ground_energy())"
   ]
  }
 ],
 "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.5.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}