docs/tutorials/circuits_3_arbitrary_basis_trotter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rz159NaSyQNW"
   },
   "source": [
    "##### Copyright 2020 The OpenFermion Developers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "id": "XwmPOS8rySzo"
   },
   "outputs": [],
   "source": [
    "#@title Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9f8d4dd70ff9"
   },
   "source": [
    "# Circuits 3: Low rank, arbitrary basis molecular simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E7CNq-HzyXBY"
   },
   "source": [
    "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://quantumai.google/openfermion/tutorials/circuits_3_arbitrary_basis_trotter\"><img src=\"https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png\" />View on QuantumAI</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/quantumlib/OpenFermion/blob/master/docs/tutorials/circuits_3_arbitrary_basis_trotter.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://github.com/quantumlib/OpenFermion/blob/master/docs/tutorials/circuits_3_arbitrary_basis_trotter.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/github_logo_1x.png\" />View source on GitHub</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a href=\"https://storage.googleapis.com/tensorflow_docs/OpenFermion/docs/tutorials/circuits_3_arbitrary_basis_trotter.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/download_icon_1x.png\" />Download notebook</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wJqrHiNcyrMH"
   },
   "source": [
    "## Setup\n",
    "\n",
    "Install the OpenFermion package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5cbe6b680387"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import openfermion\n",
    "except ImportError:\n",
    "    !pip install git+https://github.com/quantumlib/OpenFermion.git@master#egg=openfermion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "90c6d9693f39"
   },
   "source": [
    "## Low rank decomposition of the Coulomb operator\n",
    "\n",
    "The algorithm discussed in this tutorial is described in [arXiv:1808.02625](https://arxiv.org/abs/1808.02625).\n",
    "\n",
    "In [Circuits 1](./circuits_1_basis_change.ipynb) we discussed methods for very compiling single-particle basis transformations of fermionic operators in $O(N)$ depth on a linearly connected architecture. We looked at the particular example of simulating a free fermion model by using Bogoliubov transformations to diagonalize the model.\n",
    "\n",
    "In [Circuits 2](./circuits_2_diagonal_coulomb_trotter.ipynb) we discussed methods for compiling Trotter steps of electronic structure Hamiltonian in $O(N)$ depth on a linearly connected architecture when expressed in a basis diagonalizing the Coulomb operator so that\n",
    "$$\n",
    "H = \\sum_{pq} T_{pq} a^\\dagger_p a_q + \\sum_{pq} V_{pq} a^\\dagger_p a_p a^\\dagger_q a_q.\n",
    "$$\n",
    "\n",
    "Here we will discuss how both of those techniques can be combined, along with some insights from electronic structure, in order to simulate arbitrary basis molecular Hamiltonians taking the form\n",
    "$$\n",
    "H = \\sum_{pq} T_{pq} a^\\dagger_p a_q + \\sum_{pqrs} V_{pqrs} a^\\dagger_p a_q a^\\dagger_r a_s\n",
    "$$\n",
    "in depth scaling only as $O(N^2)$ on a linear array of qubits. First, we note that the one-body part of the above expression is easy to simulate using the techniques introduced in [Circuits 1](./circuits_1_basis_change.ipynb). Thus, the real challenge is to simulate the two-body part of the operator.\n",
    "\n",
    "We begin with the observation that the rank-4 tensor $V$, with the values $V_{pqrs}$ representing the coefficient of $a^\\dagger_p a_q a^\\dagger_r a_s$ can be flattened into an $N^2 \\times N^2$ array by making $p,q$ one index and $r,s$ the other. This is the electronic repulsion integral (ERI) matrix in chemist notation. We will refer to the ERI matrix as $W$. By diagonalizing $W$, one obtains $W g_\\ell = \\lambda_\\ell g_\\ell$ where the eigenvector $g_\\ell$ is a vector of dimension $N^2$. If we reshape $g_\\ell$ into an $N \\times N$ vector, we realize that\n",
    "$$\n",
    "\\sum_{pqrs} V_{pqrs} a^\\dagger_p a_q a^\\dagger_r a_s = \\sum_{\\ell=0}^{L-1} \\lambda_\\ell \\left(\\sum_{pq} \\left[g_{\\ell}\\right]_{pq} a^\\dagger_p a_q\\right)^2.\n",
    "$$\n",
    "This is related to the concept of density fitting in electronic structure, which is often accomplished using a Cholesky decomposition. It is fairly well known in the quantum chemistry community that the ERI matrix is positive semi-definite and despite having linear dimension $N^2$, has rank of only $L = O(N)$. Thus, the eigenvalues $\\lambda_\\ell$ are positive and there are only $O(N)$ of them.\n",
    "\n",
    "Next, we diagonalize the one-body operators inside of the square so that\n",
    "$$\n",
    "R_\\ell \\left(\\sum_{pq} \\left[g_\\ell\\right]_{pq} a^\\dagger_p a_q\\right) R_\\ell^\\dagger = \\sum_{p} f_{\\ell p} a^\\dagger_p a_p\n",
    "$$\n",
    "where the $R_\\ell$ represent single-particle basis transformations of the sort we compiled in [Circuits 1](./circuits_1_basis_change.ipynb). Then,\n",
    "$$\n",
    "\\sum_{\\ell=0}^{L-1} \\lambda_\\ell \\left(\\sum_{pq} \\left[g_{\\ell}\\right]_{pq} a^\\dagger_p a_q\\right)^2 =\n",
    "\\sum_{\\ell=0}^{L-1} \\lambda_\\ell \\left(R_\\ell \\left(\\sum_{p} f_{\\ell p} a^\\dagger_p a_p\\right) R_\\ell^\\dagger\\right)^2 = \\sum_{\\ell=0}^{L-1} \\lambda_\\ell \\left(R_\\ell \\left(\\sum_{p} f_{\\ell p} a^\\dagger_p a_p\\right) R_\\ell^\\dagger R_\\ell \\left(\\sum_{p} f_{\\ell p} a^\\dagger_p a_p\\right) R_\\ell^\\dagger\\right)\n",
    "= \\sum_{\\ell=0}^{L-1} \\lambda_\\ell R_\\ell \\left(\\sum_{pq} f_{\\ell p} f_{\\ell q} a^\\dagger_p a_p a^\\dagger_q a_q\\right) R_\\ell^\\dagger.\n",
    "$$\n",
    "We now see that we can simulate a Trotter step under the arbitrary basis two-body operator as\n",
    "$$\n",
    "\\prod_{\\ell=0}^{L-1} R_\\ell \\exp\\left(-i\\sum_{pq} f_{\\ell p} f_{\\ell q} a^\\dagger_p a_p a^\\dagger_q a_q\\right) R_\\ell^\\dagger\n",
    "$$\n",
    "where we note that the operator in the exponential take the form of a diagonal Coulomb operator. Since we can implement the $R_\\ell$ circuits in $O(N)$ depth (see [Circuits 1](./circuits_1_basis_change.ipynb)) and we can implement Trotter steps under diagonal Coulomb operators in $O(N)$ layers of gates (see [Circuits 2](./circuits_2_diagonal_coulomb_trotter.ipynb)) we see that we can implement Trotter steps under arbitrary basis electronic structure Hamiltionians in $O(L N) = O(N^2)$ depth, and all on a linearly connected device. This is a big improvement over the usual way of doing things, which would lead to no less than $O(N^5)$ depth! In fact, it is also possible to do better by truncating rank on the second diagonalization but we have not implemented that (details will be discussed in aforementioned paper-in-preparation).\n",
    "\n",
    "Note that these techniques are also applicable to realizing evolution under other two-body operators, such as the generator of unitary coupled cluster. Note that one can create variational algorithms where a variational parameter specifies the rank at which to truncate the $\\lambda_\\ell$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0800a1366c50"
   },
   "source": [
    "## Example implementation: Trotter steps of LiH in molecular orbital basis\n",
    "\n",
    "We will now use these techniques to implement Trotter steps for an actual molecule. We will focus on LiH at equilibrium geometry, since integrals for that system are provided with every OpenFermion installation. However, by installing [OpenFermion-PySCF](https://github.com/quantumlib/OpenFermion-PySCF) or [OpenFermion-Psi4](https://github.com/quantumlib/OpenFermion-Psi4) one can use these techniques for any molecule at any geometry. We will generate LiH in an active space consisting of 4 qubits. First, we obtain the Hamiltonian as an InteractionOperator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9b51e47b58f7"
   },
   "outputs": [],
   "source": [
    "import openfermion\n",
    "\n",
    "# Set Hamiltonian parameters for LiH simulation in active space.\n",
    "diatomic_bond_length = 1.45\n",
    "geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., diatomic_bond_length))]\n",
    "basis = 'sto-3g'\n",
    "multiplicity = 1\n",
    "active_space_start = 1\n",
    "active_space_stop = 3\n",
    "\n",
    "# Generate and populate instance of MolecularData.\n",
    "molecule = openfermion.MolecularData(geometry, basis, multiplicity, description=\"1.45\")\n",
    "molecule.load()\n",
    "\n",
    "# Get the Hamiltonian in an active space.\n",
    "molecular_hamiltonian = molecule.get_molecular_hamiltonian(\n",
    "    occupied_indices=range(active_space_start),\n",
    "    active_indices=range(active_space_start, active_space_stop))\n",
    "print(openfermion.get_fermion_operator(molecular_hamiltonian))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9265fc9c3dae"
   },
   "source": [
    "We see from the above output that this is a fairly complex Hamiltonian already. Next we will use the `simulate_trotter` function from [Circuits 1](./circuits_1_basis_change.ipynb), but this time using a different type of Trotter step associated with these low rank techniques. To keep this circuit very short for pedagogical purposes we will force a truncation of the eigenvalues $\\lambda_\\ell$ at a predetermined value of `final_rank`. While we also support a canned `LOW_RANK` option for the Trotter steps, in order to pass this value of `final_rank` we will instantiate a custom Trotter algorithm type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d530047fea8c"
   },
   "outputs": [],
   "source": [
    "import cirq\n",
    "import openfermion\n",
    "from openfermion.circuits import trotter\n",
    "\n",
    "# Trotter step parameters.\n",
    "time = 1.\n",
    "final_rank = 2\n",
    "\n",
    "# Initialize circuit qubits in a line.\n",
    "n_qubits = openfermion.count_qubits(molecular_hamiltonian)\n",
    "qubits = cirq.LineQubit.range(n_qubits)\n",
    "\n",
    "# Compile the low rank Trotter step using OpenFermion.\n",
    "custom_algorithm = trotter.LowRankTrotterAlgorithm(final_rank=final_rank)\n",
    "circuit = cirq.Circuit(\n",
    "    trotter.simulate_trotter(\n",
    "            qubits, molecular_hamiltonian,\n",
    "            time=time, omit_final_swaps=True,\n",
    "            algorithm=custom_algorithm),\n",
    "    strategy=cirq.InsertStrategy.EARLIEST)\n",
    "\n",
    "# Print circuit.\n",
    "cirq.drop_negligible_operations(circuit)\n",
    "print(circuit.to_text_diagram(transpose=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "40f8df9e4f45"
   },
   "source": [
    "We were able to print out the circuit this way but forcing `final_rank` of 2 is not very accurate. In the cell below, we compile the Trotter step with full rank so $L = N^2$ and depth is actually $O(N^3)$ and repeat the Trotter step multiple times to show that it actually converges to the correct result. Since we are not forcing the rank truncation we can use the built-in `LOW_RANK` Trotter step type. Note that the rank of the Coulomb operators is asymptotically $O(N)$ but for very small molecules in small basis sets only a few eigenvalues can be truncated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e2f406181547"
   },
   "outputs": [],
   "source": [
    "# Initialize a random initial state.\n",
    "import numpy\n",
    "random_seed = 8317\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(molecular_hamiltonian)\n",
    "exact_state = scipy.sparse.linalg.expm_multiply(\n",
    "    -1j * time * hamiltonian_sparse, initial_state)\n",
    "\n",
    "# Trotter step parameters.\n",
    "n_steps = 3\n",
    "\n",
    "# Compile the low rank Trotter step using OpenFermion.\n",
    "qubits = cirq.LineQubit.range(n_qubits)\n",
    "circuit = cirq.Circuit(\n",
    "    trotter.simulate_trotter(\n",
    "            qubits, molecular_hamiltonian,\n",
    "            time=time, n_steps=n_steps,\n",
    "            algorithm=trotter.LOW_RANK),\n",
    "    strategy=cirq.InsertStrategy.EARLIEST)\n",
    "\n",
    "# Use Cirq simulator to apply circuit.\n",
    "simulator = cirq.Simulator()\n",
    "result = simulator.simulate(circuit, qubit_order=qubits, initial_state=initial_state)\n",
    "simulated_state = result.final_state_vector\n",
    "\n",
    "# Print final fidelity.\n",
    "fidelity = abs(numpy.dot(simulated_state, numpy.conjugate(exact_state))) ** 2\n",
    "print('Fidelity with exact result is {}.\\n'.format(fidelity))\n",
    "\n",
    "# Print circuit.\n",
    "cirq.drop_negligible_operations(circuit)\n",
    "print(circuit.to_text_diagram(transpose=True))"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "circuits_3_arbitrary_basis_trotter.ipynb",
   "toc_visible": true
  },
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}