examples/tutorial_1_basis_change.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 1: Compiling arbitrary single-particle basis rotations in linear depth\n",
    "This is the first of several tutorial explaining key concepts in OpenFermion-Cirq. These tutorials build on one another and should be studied in order. In this tutorial we will discuss the compilation of circuits for implementing arbitrary rotations of the single-particle basis of an electronic structure simulation. As an example, we show how one can use these methods to simulate the evolution of an arbitrary non-interacting fermion model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "### Second quantized fermionic operators\n",
    "\n",
    "In order to represent fermionic systems on a quantum computer one must first discretize space. Usually, one expands the many-body wavefunction in a basis of spin-orbitals $\\varphi_p = \\varphi_p(r)$ which are single-particle basis functions. For reasons of spatial efficiency, all NISQ (and even most error-corrected) algorithms for simulating fermionic systems focus on representing operators in second-quantization. Second-quantized operators are expressed using the fermionic creation and annihilation operators, $a^\\dagger_p$ and $a_p$. The action of $a^\\dagger_p$ is to excite a fermion in spin-orbital $\\varphi_p$ and the action of $a_p$ is to annihilate a fermion from spin-orbital $\\varphi_p$. Specifically, if electron $i$ is represented in a space of spin-orbitals $\\{\\varphi_p(r_i)\\}$ then $a^\\dagger_p$ and $a_p$ are related to Slater determinants through the equivalence,\n",
    "$$\n",
    "\\langle r_0 \\cdots r_{\\eta-1} | a^\\dagger_{0} \\cdots a^\\dagger_{\\eta-1} | \\varnothing\\rangle \\equiv \\sqrt{\\frac{1}{\\eta!}}\n",
    "\\begin{vmatrix}\n",
    "\\varphi_{0}\\left(r_0\\right) & \\varphi_{1}\\left( r_0\\right) & \\cdots & \\varphi_{\\eta-1} \\left( r_0\\right) \\\\\n",
    "\\varphi_{0}\\left(r_1\\right) & \\varphi_{1}\\left( r_1\\right) & \\cdots & \\varphi_{\\eta-1} \\left( r_1\\right) \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots\\\\\n",
    "\\varphi_{0}\\left(r_{\\eta-1}\\right) & \\varphi_{1}\\left(r_{\\eta-1}\\right) & \\cdots & \\varphi_{\\eta-1} \\left(r_{\\eta-1}\\right) \\end{vmatrix}\n",
    "$$\n",
    "where $\\eta$ is the number of electrons in the system, $|\\varnothing \\rangle$ is the Fermi vacuum and $\\varphi_p(r)=\\langle r|\\varphi_p \\rangle$ are the single-particle orbitals that define the basis. By using a basis of Slater determinants, we ensure antisymmetry in the encoded state.\n",
    "\n",
    "### Rotations of the single-particle basis\n",
    "\n",
    "Very often in electronic structure calculations one would like to rotate the single-particle basis. That is, one would like to generate new orbitals that are formed from a linear combination of the old orbitals. Any particle-conserving rotation of the single-particle basis can be expressed as\n",
    "$$\n",
    "\\tilde{\\varphi}_p = \\sum_{q} \\varphi_q u_{pq}\n",
    "\\quad\n",
    "\\tilde{a}^\\dagger_p = \\sum_{q} a^\\dagger_q u_{pq}\n",
    "\\quad\n",
    "\\tilde{a}_p = \\sum_{q} a_q u_{pq}^*\n",
    "$$\n",
    "where $\\tilde{\\varphi}_p$, $\\tilde{a}^\\dagger_p$, and $\\tilde{a}^\\dagger_p$ correspond to spin-orbitals and operators in the rotated basis and $u$ is an $N\\times N$ unitary matrix. From the Thouless theorem, this single-particle rotation\n",
    "is equivalent to applying the $2^N \\times 2^N$ operator\n",
    "$$\n",
    "  U(u) = \\exp\\left(\\sum_{pq} \\left[\\log u \\right]_{pq} \\left(a^\\dagger_p a_q - a^\\dagger_q a_p\\right)\\right) \n",
    "$$\n",
    "where $\\left[\\log u\\right]_{pq}$ is the $(p, q)$ element of the matrix $\\log u$.\n",
    "\n",
    "There are many reasons that one might be interested in performing such basis rotations. For instance, one might be interested in preparing the Hartree-Fock (mean-field) state of a chemical system, by rotating from some initial orbitals (e.g. atomic orbitals or plane waves) into the molecular orbitals of the system. Alternatively, one might be interested in rotating from a basis where certain operators are diagonal (e.g. the kinetic operator is diagonal in the plane wave basis) to a basis where certain other operators are diagonal (e.g. the Coulomb operator is diagonal in the position basis). Thus, it is a very useful thing to be able to apply circuits corresponding to $U(u)$ on a quantum computer in low depth.\n",
    "\n",
    "### Compiling linear depth circuits to rotate the orbital basis\n",
    "\n",
    "OpenFermion-Cirq prominently features routines for implementing the linear depth / linear connectivity basis transformations described in [Phys. Rev. Lett. 120, 110501](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.110501). While we will not discuss this functionality here, we also support routines for compiling the more general form of these transformations which do not conserve particle-number, known as a Bogoliubov transfomation, using routines described in [Phys. Rev. Applied 9, 044036](https://journals.aps.org/prapplied/abstract/10.1103/PhysRevApplied.9.044036). We will not discuss the details of how these methods are implemented here and instead refer readers to those papers. All that one needs in order to compile the circuit $U(u)$ using OpenFermion-Cirq is the $N \\times N$ matrix $u$, which we refer to in documentation as the \"basis_transformation_matrix\". Note that if one intends to apply this matrix to a computational basis state with only $\\eta$ electrons, then one can reduce the number of gates required by instead supplying the $\\eta \\times N$ rectangular matrix that characterizes the rotation of the occupied orbitals only. OpenFermion-Cirq will automatically take advantage of this symmetry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenFermion-Cirq example implementation: exact evolution under tight binding models\n",
    "In this example will show how basis transforms can be used to implement exact evolution under a random Hermitian one-body fermionic operator\n",
    "\\begin{equation}\n",
    "H = \\sum_{pq} T_{pq} a^\\dagger_p a_q.\n",
    "\\end{equation}\n",
    "That is, we will compile a circuit to implement $e^{-i H t}$ for some time $t$. Of course, this is a tractable problem classically but we discuss it here since it is often useful as a subroutine for more complex quantum simulations. To accomplish this evolution, we will use basis transformations. Suppose that $u$ is the basis transformation matrix that diagonalizes $T$. Then, we could implement $e^{-i H t}$ by implementing $U(u)^\\dagger (\\prod_{k} e^{-i \\lambda_k Z_k}) U(u)$ where $\\lambda_k$ are the eigenvalues of $T$. \n",
    "\n",
    "Below, we initialize the T matrix characterizing $H$ and then obtain the eigenvalues $\\lambda_k$ and eigenvectors $u_k$ of $T$. We print out the OpenFermion FermionOperator representation of $T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.5367212624097257+0j) [0^ 0] +\n",
      "(-0.26033703159240107+3.3259173741375454j) [0^ 1] +\n",
      "(1.3433603748462144+1.544987250567917j) [0^ 2] +\n",
      "(-0.26033703159240107-3.3259173741375454j) [1^ 0] +\n",
      "(-2.9143303700812435+0j) [1^ 1] +\n",
      "(-1.52843836446248+1.3527486791390022j) [1^ 2] +\n",
      "(1.3433603748462144-1.544987250567917j) [2^ 0] +\n",
      "(-1.52843836446248-1.3527486791390022j) [2^ 1] +\n",
      "(2.261633626116526+0j) [2^ 2]\n"
     ]
    }
   ],
   "source": [
    "import openfermion\n",
    "import numpy\n",
    "\n",
    "# Set the number of qubits in our example.\n",
    "n_qubits = 3\n",
    "simulation_time = 1.\n",
    "random_seed = 8317\n",
    "\n",
    "# Generate the random one-body operator.\n",
    "T = openfermion.random_hermitian_matrix(n_qubits, seed=random_seed)\n",
    "\n",
    "# Diagonalize T and obtain basis transformation matrix (aka \"u\").\n",
    "eigenvalues, eigenvectors = numpy.linalg.eigh(T)\n",
    "basis_transformation_matrix = eigenvectors.transpose()\n",
    "\n",
    "# Print out familiar OpenFermion \"FermionOperator\" form of H.\n",
    "H = openfermion.FermionOperator()\n",
    "for p in range(n_qubits):\n",
    "    for q in range(n_qubits):\n",
    "        term = ((p, 1), (q, 0))\n",
    "        H += openfermion.FermionOperator(term, T[p, q])\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to make a circuit! First we will use OpenFermion-Cirq to generate the basis transform $U(u)$ from the basis transformation matrix $u$ by calling the Bogoliubov transform function (named as such because this function can also handle non-particle conserving basis transformations). Then, we'll apply local $Z$ rotations to phase by the eigenvalues, then we'll apply the inverse transformation. That will finish the circuit. We're just going to print out the first rotation to keep things easy-to-read, but feel free to play around with the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: ────────────────────────────────────YXXY────────Rz(-0.0π)──────────────\n",
      "                                       │\n",
      "1: ─────────────YXXY────────Z^0.522────#2^-0.656───YXXY────────Rz(0.0π)───\n",
      "                │                                  │\n",
      "2: ───Z^0.762───#2^-0.249───Z^-0.519───────────────#2^-0.479───Rz(-π)─────\n"
     ]
    }
   ],
   "source": [
    "import openfermioncirq\n",
    "import cirq\n",
    "\n",
    "# Initialize the qubit register.\n",
    "qubits = cirq.LineQubit.range(n_qubits)\n",
    "\n",
    "# Start circuit with the inverse basis rotation, print out this step.\n",
    "inverse_basis_rotation = cirq.inverse(openfermioncirq.bogoliubov_transform(qubits, basis_transformation_matrix))\n",
    "circuit = cirq.Circuit(inverse_basis_rotation)\n",
    "print(circuit)\n",
    "\n",
    "# Add diagonal phase rotations to circuit.\n",
    "for k, eigenvalue in enumerate(eigenvalues):\n",
    "    phase = -eigenvalue * simulation_time\n",
    "    circuit.append(cirq.Rz(rads=phase).on(qubits[k]))\n",
    "\n",
    "# Finally, restore basis.\n",
    "basis_rotation = openfermioncirq.bogoliubov_transform(qubits, basis_transformation_matrix)\n",
    "circuit.append(basis_rotation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can check whether our circuit applied to a random initial state with the exact result. Print out the fidelity with the exact result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000000960717732\n"
     ]
    }
   ],
   "source": [
    "# Initialize a random initial state.\n",
    "initial_state = openfermion.haar_random_vector(\n",
    "    2 ** n_qubits, random_seed).astype(numpy.complex64)\n",
    "\n",
    "# Numerically compute the correct circuit output.\n",
    "import scipy\n",
    "hamiltonian_sparse = openfermion.get_sparse_operator(H)\n",
    "exact_state = scipy.sparse.linalg.expm_multiply(\n",
    "    -1j * simulation_time * hamiltonian_sparse, initial_state)\n",
    "\n",
    "# Use Cirq simulator to apply circuit.\n",
    "simulator = cirq.Simulator()\n",
    "result = simulator.simulate(circuit, qubit_order=qubits,\n",
    "                            initial_state=initial_state)\n",
    "simulated_state = result.final_state\n",
    "\n",
    "# Print final fidelity.\n",
    "fidelity = abs(numpy.dot(simulated_state, numpy.conjugate(exact_state)))**2\n",
    "print(fidelity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we see that the circuit correctly effects the intended evolution. We can now use Cirq's compiler to output the circuit using gates native to near-term devices, and then optimize those circuits. We'll output in QASM 2.0 just to demonstrate that functionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "// Generated from Cirq v0.6.0\n",
      "\n",
      "OPENQASM 2.0;\n",
      "include \"qelib1.inc\";\n",
      "\n",
      "\n",
      "// Qubits: [0, 1, 2]\n",
      "qreg q[3];\n",
      "\n",
      "\n",
      "u2(pi*-0.3412628078, pi*0.3412628078) q[2];\n",
      "u2(pi*-0.5794122432, pi*0.5794122432) q[1];\n",
      "u2(pi*0.0897056074, pi*-0.0897056074) q[0];\n",
      "cz q[1],q[2];\n",
      "u3(pi*-0.1242949803, pi*0.6587371922, pi*-0.6587371922) q[2];\n",
      "u3(pi*-0.1242949803, pi*1.4205877568, pi*-1.4205877568) q[1];\n",
      "cz q[1],q[2];\n",
      "u3(pi*-0.3154309963, pi*0.1587371922, pi*-0.1587371922) q[2];\n",
      "u3(pi*-0.8528172267, pi*-0.0794122432, pi*0.0794122432) q[1];\n",
      "cz q[0],q[1];\n",
      "u3(pi*-0.328242091, pi*0.4205877568, pi*-0.4205877568) q[1];\n",
      "u3(pi*-0.328242091, pi*0.0897056074, pi*-0.0897056074) q[0];\n",
      "cz q[0],q[1];\n",
      "u3(pi*-0.9772799497, pi*0.9205877568, pi*-0.9205877568) q[1];\n",
      "u3(pi*-0.0380943193, pi*0.5897056074, pi*-0.5897056074) q[0];\n",
      "cz q[1],q[2];\n",
      "u3(pi*-0.2326621647, pi*-0.3412628078, pi*0.3412628078) q[2];\n",
      "u3(pi*-0.2326621647, pi*0.4205877568, pi*-0.4205877568) q[1];\n",
      "cz q[1],q[2];\n",
      "u3(pi*-0.336801226, pi*1.1587371922, pi*-1.1587371922) q[2];\n",
      "u3(pi*-0.7711683715, pi*0.9205877568, pi*-0.9205877568) q[1];\n",
      "cz q[0],q[1];\n",
      "u3(pi*-0.328242091, pi*0.4205877568, pi*-0.4205877568) q[1];\n",
      "u3(pi*-0.328242091, pi*0.0897056074, pi*-0.0897056074) q[0];\n",
      "cz q[0],q[1];\n",
      "u3(pi*-0.9123845954, pi*-0.0794122432, pi*0.0794122432) q[1];\n",
      "u2(pi*-0.9102943926, pi*0.9102943926) q[0];\n",
      "rz(pi*-0.1553859346) q[0];\n",
      "cz q[1],q[2];\n",
      "u3(pi*-0.1242949803, pi*-0.3412628078, pi*0.3412628078) q[2];\n",
      "u3(pi*-0.1242949803, pi*0.4205877568, pi*-0.4205877568) q[1];\n",
      "cz q[1],q[2];\n",
      "u2(pi*-1.3412628078, pi*1.3412628078) q[2];\n",
      "u2(pi*0.4205877568, pi*-0.4205877568) q[1];\n",
      "rz(pi*1.0794122432) q[1];\n",
      "rz(pi*-0.9205877568) q[2];\n",
      "\n"
     ]
    }
   ],
   "source": [
    "xmon_circuit = cirq.google.optimized_for_xmon(circuit)\n",
    "print(xmon_circuit.to_qasm())"
   ]
  }
 ],
 "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
}