CIRQ_Openfermion.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CIRQ_Openfermion.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "CKF-xHLLLojw"
      },
      "source": [
        "!pip install cirq\n",
        "!pip install xlrd\n",
        "!pip install pysf\n",
        "!pip install openfermion\n",
        "!pip install openfermionpyscf\n",
        "import numpy as np"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dMZgI_A5_DJX"
      },
      "source": [
        "import cirq\n",
        "import scipy\n",
        "import pyscf\n",
        "import openfermion\n",
        "import openfermionpyscf \n",
        "import itertools\n",
        "\n",
        "from pyscf import gto, scf, fci\n",
        "from openfermion.transforms import jordan_wigner, get_fermion_operator\n",
        "from openfermionpyscf import run_pyscf\n",
        "from openfermion.utils import count_qubits\n",
        "from openfermion.linalg import jw_hartree_fock_state\n",
        "from openfermion.circuits import simulate_trotter\n",
        "from openfermion import (\n",
        "    get_sparse_operator, get_ground_state, FermionOperator,\n",
        "    jw_get_ground_state_at_particle_number, MolecularData,\n",
        "    expectation, uccsd_convert_amplitude_format,\n",
        "    get_interaction_operator, QubitOperator, eigenspectrum,\n",
        "    InteractionOperator\n",
        ")"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uXr004duTvwT"
      },
      "source": [
        "## Function Definitions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-BbDAfpQpORJ"
      },
      "source": [
        "chemicalAccuracy = 1.5936*10**-3\n",
        "\n",
        "# Define necessary Pauli operators (two-dimensional) as matrices\n",
        "pauliX = np.array([[0,1],\n",
        "                 [1,0]],\n",
        "                dtype = complex)\n",
        "pauliZ = np.array([[1,0],\n",
        "                 [0,-1]],\n",
        "                dtype = complex)\n",
        "pauliY = np.array([[0,-1j],\n",
        "                 [1j,0]],\n",
        "                dtype = complex)\n",
        "\n",
        "def stringToMatrix(pauliString):\n",
        "  '''\n",
        "  Converts a Pauli string to its matrix form.\n",
        "\n",
        "  Arguments:\n",
        "    pauliString (str): the Pauli string (e.g. \"IXYIZ\")\n",
        "\n",
        "  Returns:\n",
        "    matrix (np.ndarray): the corresponding matrix, in the computational basis\n",
        "\n",
        "  '''\n",
        "\n",
        "  matrix = np.array([1])\n",
        "\n",
        "  # Iteratively construct the matrix, going through each single qubit Pauli term\n",
        "  for pauli in pauliString:\n",
        "      if pauli == \"I\":\n",
        "        matrix = np.kron(matrix,np.identity(2))\n",
        "      elif pauli == \"X\":\n",
        "        matrix = np.kron(matrix,pauliX)\n",
        "      elif pauli == \"Y\":\n",
        "        matrix = np.kron(matrix,pauliY)\n",
        "      elif pauli == \"Z\":\n",
        "        matrix = np.kron(matrix,pauliZ)\n",
        "\n",
        "  return matrix\n",
        "  \n",
        "def fromVectortoKet(stateVector):\n",
        "  '''\n",
        "  Transforms a vector representing a basis state to the corresponding ket.\n",
        "\n",
        "  Arguments:\n",
        "    stateVector (np.ndarray): basis vector in the 2^n-dimensional Hilbert space\n",
        "\n",
        "  Returns:\n",
        "    ket (list): a list of length n representing the corresponding ket \n",
        "  '''\n",
        "\n",
        "  dim = len(stateVector)\n",
        "  ket = []\n",
        "\n",
        "  while dim>1:\n",
        "    if any (stateVector[i] for i in range(int(dim/2))):\n",
        "      # Ket is of the form |0>|...>. \n",
        "\n",
        "      #Fix |0> as the msq.\n",
        "      ket.append(0)\n",
        "\n",
        "      # Get the vector representing the state of the remaining qubits.\n",
        "      stateVector = stateVector[:int(dim/2)]\n",
        "\n",
        "    else:\n",
        "      # Ket is of the form |1>|...>. \n",
        "      \n",
        "      #Fix |0> as the msq.\n",
        "      ket.append(1)\n",
        "\n",
        "      # Get the vector representing the state of the remaining qubits.\n",
        "      stateVector = stateVector[int(dim//2):]\n",
        "\n",
        "    dim = dim/2\n",
        "\n",
        "  return ket\n",
        "\n",
        "def fromKettoVector(ket):\n",
        "  '''\n",
        "  Transforms a ket representing a basis state to the corresponding state vector.\n",
        "\n",
        "  Arguments:\n",
        "    ket (list): a list of length n representing the ket \n",
        "\n",
        "  Returns:\n",
        "    stateVector (np.ndarray): the corresponding basis vector in the \n",
        "      2^n dimensional Hilbert space\n",
        "  '''\n",
        "  stateVector = [1]\n",
        "\n",
        "  # Iterate through the ket, calculating the tensor product of the qubit states\n",
        "  for i in ket:\n",
        "    qubitVector = [not i,i]\n",
        "    stateVector = np.kron(stateVector,qubitVector)\n",
        "\n",
        "  return stateVector\n",
        "\n",
        "def calculateOverlap(stateCoordinates1,stateCoordinates2):\n",
        "    '''\n",
        "    Calculates the overlap between two states, given their coordinates.\n",
        "\n",
        "    Arguments:\n",
        "      stateCoordinates1 (np.ndarray): the coordinates of one of the states in \n",
        "        some orthonormal basis,\n",
        "      stateCoordinates2 (np.ndarray): the coordinates of the other state, in \n",
        "        the same basis\n",
        "\n",
        "    Returns: \n",
        "      overlap (float): the overlap between two states (absolute value of the \n",
        "        inner product).\n",
        "    '''\n",
        "\n",
        "    bra = np.conj(stateCoordinates1)\n",
        "    ket = stateCoordinates2\n",
        "    overlap = np.abs(np.dot(bra,ket))\n",
        "    \n",
        "    return overlap\n",
        "\n",
        "def findSubStrings(mainString,hamiltonian,checked = []):\n",
        "    '''\n",
        "    Finds and groups all the strings in a Hamiltonian that only differ from \n",
        "    mainString by identity operators.\n",
        "\n",
        "    Arguments:\n",
        "      mainString (str): a Pauli string (e.g. \"XZ)\n",
        "      hamiltonian (dict): a Hamiltonian (with Pauli strings as keys and their \n",
        "        coefficients as values)\n",
        "      checked (list): a list of the strings in the Hamiltonian that have already\n",
        "        been inserted in another group\n",
        "\n",
        "    Returns: \n",
        "      groupedOperators (dict): a dictionary whose keys are boolean strings \n",
        "        representing substrings of the mainString (e.g. if mainString = \"XZ\", \n",
        "        \"IZ\" would be represented as \"01\"). It includes all the strings in the \n",
        "        hamiltonian that can be written in this form (because they only differ \n",
        "        from mainString by identities), except for those that were in checked\n",
        "        (because they are already part of another group of strings).\n",
        "      checked (list):  the same list passed as an argument, with extra values\n",
        "        (the strings that were grouped in this function call).\n",
        "    '''\n",
        "    \n",
        "    groupedOperators = {}\n",
        "    \n",
        "    # Go through the keys in the dictionary representing the Hamiltonian that \n",
        "    #haven't been grouped yet, and find those that only differ from mainString \n",
        "    #by identities\n",
        "    for pauliString in hamiltonian:\n",
        "        \n",
        "        if pauliString not in checked:\n",
        "            # The string hasn't been grouped yet\n",
        "            \n",
        "            if(all((op1 == op2 or op2 == \"I\") \\\n",
        "                   for op1,op2 in zip(mainString,pauliString))):\n",
        "                # The string only differs from mainString by identities\n",
        "                \n",
        "                # Represent the string as a substring of the main one\n",
        "                booleanString = \"\".join([str(int(op1 == op2)) for op1,op2 in \\\n",
        "                                       zip(mainString,pauliString)])\n",
        "                    \n",
        "                # Add the boolean string representing this string as a key to \n",
        "                #the dictionary of grouped operators, and associate its \n",
        "                #coefficient as its value\n",
        "                groupedOperators[booleanString] = hamiltonian[pauliString]\n",
        "                \n",
        "                # Mark the string as grouped, so that it's not added to any \n",
        "                #other group\n",
        "                checked.append(pauliString)\n",
        "                \n",
        "    return (groupedOperators,checked)\n",
        "\n",
        "def groupHamiltonian(hamiltonian):\n",
        "    '''\n",
        "    Organizes a Hamiltonian into groups where strings only differ from \n",
        "    identities, so that the expectation values of all the strings in each \n",
        "    group can be calculated from the same measurement array.\n",
        "\n",
        "    Arguments: \n",
        "      hamiltonian (dict): a dictionary representing a Hamiltonian, with Pauli \n",
        "        strings as keys and their coefficients as values.\n",
        "\n",
        "    Returns: \n",
        "      groupedHamiltonian (dict): a dictionary of subhamiltonians, each of \n",
        "        which includes Pauli strings that only differ from each other by \n",
        "        identities. \n",
        "        The keys of groupedHamiltonian are the main strings of each group: the \n",
        "        ones with least identity terms. The value associated to a main string is \n",
        "        a dictionary, whose keys are boolean strings representing substrings of \n",
        "        the respective main string (with 1 where the Pauli is the same, and 0\n",
        "        where it's identity instead). The values are their coefficients.\n",
        "    '''\n",
        "    groupedHamiltonian = {}\n",
        "    checked = []\n",
        "    \n",
        "    # Go through the hamiltonian, starting by the terms that have less\n",
        "    #identity operators\n",
        "    for mainString in \\\n",
        "        sorted(hamiltonian,key = lambda pauliString: pauliString.count(\"I\")):\n",
        "            \n",
        "        # Call findSubStrings to find all the strings in the dictionary that \n",
        "        #only differ from mainString by identities, and organize them as a \n",
        "        #dictionary (groupedOperators)\n",
        "        groupedOperators,checked = findSubStrings(mainString,hamiltonian,checked)\n",
        "        \n",
        "        # Use the dictionary as a value for the mainString key in the \n",
        "        #groupedHamiltonian dictionary\n",
        "        groupedHamiltonian[mainString] = groupedOperators\n",
        "        \n",
        "        # If all the strings have been grouped, exit the for cycle\n",
        "        if(len(checked) == len(hamiltonian.keys())):\n",
        "           break\n",
        "       \n",
        "    return groupedHamiltonian \n",
        "    \n",
        "def convertHamiltonian(openfermionHamiltonian):\n",
        "  '''\n",
        "  Formats a qubit Hamiltonian obtained from openfermion, so that it's a suitable\n",
        "  argument for functions such as measureExpectationEstimation.\n",
        "\n",
        "  Arguments:\n",
        "    openfermionHamiltonian (openfermion.qubitOperator): the Hamiltonian.\n",
        "\n",
        "  Returns:\n",
        "    formattedHamiltonian (dict): the Hamiltonian as a dictionary with Pauli\n",
        "      strings (eg 'YXZI') as keys and their coefficients as values.\n",
        "  '''\n",
        "\n",
        "  formattedHamiltonian = {}\n",
        "  qubitNumber = count_qubits(openfermionHamiltonian)\n",
        "\n",
        "  # Iterate through the terms in the Hamiltonian\n",
        "  for term in openfermionHamiltonian.get_operators():\n",
        "\n",
        "    operators = []\n",
        "    coefficient = list(term.terms.values())[0]\n",
        "    pauliString = list(term.terms.keys())[0]\n",
        "    previousQubit = -1\n",
        "\n",
        "    for (qubit,operator) in pauliString:\n",
        "\n",
        "      # If there are qubits in which no operations are performed, add identities \n",
        "      #as necessary, to make sure that the length of the string will match the \n",
        "      #number of qubits\n",
        "      identities = (qubit-previousQubit-1)\n",
        "      if identities>0: \n",
        "        operators.append('I'*identities)\n",
        "\n",
        "      operators.append(operator)\n",
        "      previousQubit = qubit\n",
        "    \n",
        "    # Add final identity operators if the string still doesn't have the \n",
        "    #correct length (because no operations are performed in the last qubits)\n",
        "    operators.append('I'*(qubitNumber-previousQubit-1))\n",
        "\n",
        "    formattedHamiltonian[\"\".join(operators)] = coefficient\n",
        "\n",
        "  return formattedHamiltonian\n",
        "\n",
        "def hamiltonianToMatrix(hamiltonian):\n",
        "    '''\n",
        "    Convert a Hamiltonian (from OpenFermion) to matrix form.\n",
        "    \n",
        "    Arguments:\n",
        "      hamiltonian (openfermion.InteractionOperator): the Hamiltonian to be\n",
        "        transformed.\n",
        "\n",
        "    Returns:\n",
        "      matrix (np.ndarray): the Hamiltonian, as a matrix in the computational \n",
        "        basis\n",
        "    \n",
        "    ''' \n",
        "    \n",
        "    qubitNumber = hamiltonian.n_qubits\n",
        "    \n",
        "    hamiltonian = jordan_wigner(hamiltonian)\n",
        "    formattedHamiltonian = convertHamiltonian(hamiltonian)\n",
        "    groupedHamiltonian = groupHamiltonian(formattedHamiltonian)\n",
        "\n",
        "    matrix = np.zeros((2**qubitNumber,2**qubitNumber),dtype = complex)\n",
        "\n",
        "    # Iterate through the strings in the Hamiltonian, adding the respective \n",
        "    #contribution to the matrix\n",
        "    for string in groupedHamiltonian:\n",
        "\n",
        "      for substring in groupedHamiltonian[string]:\n",
        "        pauli = (\"\".join(\"I\"*(not int(b)) + a*int(b) \\\n",
        "                         for (a,b) in zip(string,substring)))\n",
        "        \n",
        "        matrix += stringToMatrix(pauli) * groupedHamiltonian[string][substring]\n",
        "\n",
        "    return matrix\n",
        "\n",
        "def groundStatesFromDiagonalization(listOfHamiltonians):\n",
        "    '''\n",
        "    Obtains the ground states and energies corresponding to a list of \n",
        "    Hamiltonians, by performing exact diagonalization.\n",
        "\n",
        "    Arguments:\n",
        "      listOfHamiltonians (list): the list of hamiltonians as dictionaries.\n",
        "\n",
        "    Returns: \n",
        "      exactGroundEnergies (list): the list of ground energies corresponding to\n",
        "        the list of Hamiltonians\n",
        "      exactGroundStates (list): the list of the ground states corresponding to \n",
        "        the list of hamiltonians\n",
        "    diagonalization.\n",
        "    '''\n",
        "    \n",
        "    exactGroundEnergies = []\n",
        "    exactGroundStates = []\n",
        "    \n",
        "    # Obtain the number of qubits and calculate the dimension.\n",
        "    n = len(list(listOfHamiltonians[0].keys())[0])\n",
        "    dim = n**2\n",
        "\n",
        "    # Find the ground state and energy of each Hamiltonian in the list.\n",
        "    for hamiltonian in listOfHamiltonians:\n",
        "        \n",
        "        # Get the matrix form of the Hamiltonian, in the computational basis.\n",
        "        hamiltonianMatrix = np.zeros((dim,dim),dtype = complex)\n",
        "\n",
        "        for pauliString in hamiltonian:\n",
        "            coefficient = hamiltonian[pauliString]\n",
        "            hamiltonianMatrix += coefficient * stringToMatrix(pauliString)\n",
        "            \n",
        "        w,v = np.linalg.eigh(hamiltonianMatrix)\n",
        "        \n",
        "        # Add the minimum eigenvalue of the Hamiltonian to the list of ground\n",
        "        #energies\n",
        "        exactGroundEnergies.append((np.amin(w)))\n",
        "\n",
        "        # Add the corresponging eigenstate to the list of ground states\n",
        "        exactGroundStates.append(v[np.argmin(w)])\n",
        "        \n",
        "    return (exactGroundEnergies,exactGroundStates)\n",
        "\n",
        "def unitaryToRotationGates(unitaryMatrix,qubit,keepPhase = False):\n",
        "    ''' \n",
        "    Decomposes a single qubit unitary matrix into rotation gates.\n",
        "\n",
        "    Arguments: \n",
        "      unitaryMatrix (np.ndarray): the matrix to be decomposed\n",
        "      qubit (cirq.LineQubit): the qubit the matrix acts on\n",
        "      keepPhase (bool): a flag indicating whether the global phase should be \n",
        "        kept\n",
        "\n",
        "    Returns: \n",
        "      (list): a list of single qubit rotation gates, equivalent to the action of \n",
        "        unitaryMatrix (up to a global phase, if keepPhase = False)\n",
        "    ''' \n",
        "    \n",
        "    # Compute the angles of each rotation from the entries of the matrix\n",
        "    rz1angle = -np.angle(unitaryMatrix[0][0])+np.angle(-unitaryMatrix[0][1])\n",
        "    ryangle = 2*np.arccos(np.abs(unitaryMatrix[0][0]))\n",
        "    rz2angle = -np.angle(unitaryMatrix[0][0])+np.angle(unitaryMatrix[1][0])\n",
        "    \n",
        "    if keepPhase:\n",
        "        gp = np.angle(unitaryMatrix[0][0])+(rz1angle+rz2angle)/2\n",
        "    \n",
        "    # Create the necessary rotation gates\n",
        "    r1 = cirq.rz(rz1angle)\n",
        "    r2 = cirq.ry(ryangle)\n",
        "    r3 = cirq.rz(rz2angle)\n",
        "    \n",
        "    yield [r1(qubit),r2(qubit),r3(qubit)]\n",
        "\n",
        "    if keepPhase:\n",
        "        # The global phase should be kept; add the correction.\n",
        "        yield cirq.MatrixGate(np.identity(2)*np.exp(1j*gp)).on(qubit)\n",
        "\n",
        "def measureString(pauliString,repetitions,statePreparationGates,qubits):\n",
        "    ''' \n",
        "    Measures the expectation value of a Pauli string using the CIRQ simulator \n",
        "    (simulating sampling).\n",
        "\n",
        "    Arguments:\n",
        "      pauliString (str): the Pauli string to be measured.\n",
        "      repetitions (int): the number of circuit repetitions to be used.\n",
        "      statePreparationGates (list): the list of CIRQ gates that prepare the \n",
        "        state in which to obtain the expectation value.\n",
        "\n",
        "    Returns: \n",
        "      expecationValue (float): the expectation value of pauliString, with \n",
        "        sampling noise.\n",
        "    ''' \n",
        "\n",
        "    # Initialize circuit.\n",
        "    circuit = cirq.Circuit(statePreparationGates)\n",
        "\n",
        "    # Optimize circuit.\n",
        "    cirq.optimizers.EjectZ().optimize_circuit(circuit)\n",
        "    cirq.optimizers.DropNegligible().optimize_circuit(circuit)\n",
        "    \n",
        "    # Append necessary rotations and measurements for each qubit.\n",
        "    for i in range(len(qubits)):\n",
        "        op = pauliString[i]\n",
        "        \n",
        "        # Rotate qubit i to the X basis if that's the desired measurement.\n",
        "        if (op == \"X\"):\n",
        "            circuit.append(cirq.H(qubits[i]))\n",
        "        \n",
        "        # Rotate qubit i to the Y basis if that's the desired measurement.\n",
        "        if (op == \"Y\"):\n",
        "            # Apply adjoint of the Clifford S gate\n",
        "            circuit.append(cirq.ops.ZPowGate(exponent = -1/2).on(qubits[i]))\n",
        "\n",
        "            # Finish rotating to Y basis, proceeding as for X\n",
        "            circuit.append(cirq.H(qubits[i]))\n",
        "            \n",
        "        # Measure qubit i in the computational basis, unless operator is I.\n",
        "        if (op != \"I\"):\n",
        "            circuit.append(cirq.measure(qubits[i],key = str(i)))\n",
        "            \n",
        "    # Sample the desired number of repetitions from the circuit, unless\n",
        "    #there are no measurements (identity term).\n",
        "    if (pauliString != \"I\"*len(qubits)):\n",
        "        s = cirq.Simulator()\n",
        "        results = s.run(circuit,repetitions = repetitions)\n",
        "    #print(circuit)\n",
        "\n",
        "    # Calculate the expectation value of the Pauli string by averaging over  \n",
        "    #all the repetitions.\n",
        "    \n",
        "    total = 0\n",
        "    \n",
        "    for j in range(repetitions):\n",
        "        meas = 1\n",
        "        for i in range(len(qubits)):\n",
        "            if (pauliString[i] != \"I\"):\n",
        "                meas = meas*(1-2*results.data[str(i)][j])\n",
        "        total += meas\n",
        "        \n",
        "    expectationValue = total/repetitions\n",
        "    \n",
        "    return expectationValue\n",
        "\n",
        "def measureExpectation(mainString,subHamiltonian,repetitions,\\\n",
        "                       statePreparationGates,qubits):\n",
        "    ''' \n",
        "    Measures the expectation value of a subHamiltonian using the CIRQ simulator \n",
        "    (simulating sampling). By construction, all the expectation values of the \n",
        "    strings in subHamiltonian can be obtained from the same measurement array.\n",
        "\n",
        "    Arguments: \n",
        "      mainString (str): the main Pauli string. This is the string in the group\n",
        "        with the least identity terms. It defines the circuit that will be used.\n",
        "      subHamiltonian (dict): a dictionary whose keys are boolean strings\n",
        "        representing substrings of the main one, and whose values are the \n",
        "        respective coefficients.\n",
        "      repetitions (int): the number of repetitions to be performed, the \n",
        "      statePreparationGates (list): the list of CIRQ gates that prepare (from \n",
        "        |0..0>) the state in which to obtain the expectation value.\n",
        "      qubits (list): list of cirq.LineQubit to apply the gates on\n",
        "\n",
        "    Returns:\n",
        "      totalExpectationValue (float): the total expectation value of \n",
        "        subHamiltonian, with sampling noise.\n",
        "    ''' \n",
        "    \n",
        "    # Initialize circuit.\n",
        "    circuit = cirq.Circuit()\n",
        "    \n",
        "    # Append to the circuit the gates that prepare the state corresponding to\n",
        "    #the received parameters.\n",
        "    circuit.append(statePreparationGates)\n",
        "    cirq.optimizers.EjectZ().optimize_circuit(circuit)\n",
        "    cirq.optimizers.DropNegligible().optimize_circuit(circuit)\n",
        "    \n",
        "    # Append necessary rotations and measurements for each qubit.\n",
        "    for i in range(len(qubits)):\n",
        "        op = mainString[i]\n",
        "        \n",
        "        # Rotate qubit i to the X basis if that's the desired measurement.\n",
        "        if (op == \"X\"):\n",
        "            circuit.append(cirq.H(qubits[i]))\n",
        "            \n",
        "        # Rotate qubit i to the X basis if that's the desired measurement.\n",
        "        if (op == \"Y\"):\n",
        "            # Apply adjoint of the Clifford S gate\n",
        "            circuit.append(cirq.ops.ZPowGate(exponent = -1/2).on(qubits[i]))\n",
        "\n",
        "            # Apply the Hadamard gate\n",
        "            circuit.append(cirq.H(qubits[i]))\n",
        "            \n",
        "        #Measure qubit i in the computational basis, unless operator is I.\n",
        "        if (op != \"I\"):\n",
        "            circuit.append(cirq.measure(qubits[i],key = str(i)))\n",
        "            \n",
        "    # Sample the desired number of repetitions from the circuit, unless\n",
        "    #there are no measurements (identity term).\n",
        "    if (mainString != \"I\"*len(qubits)):\n",
        "        s = cirq.Simulator()\n",
        "        results = s.run(circuit,repetitions = repetitions)\n",
        "\n",
        "    # For each substring, initialize the sum of all measurements as zero\n",
        "    total = {}\n",
        "    for subString in subHamiltonian:\n",
        "        total[subString] = 0\n",
        "    \n",
        "    # Calculate the expectation value of each Pauli string by averaging over  \n",
        "    #all the repetitions\n",
        "    for j in range(repetitions):\n",
        "        meas = {}\n",
        "        \n",
        "        # Initialize the measurement in repetition j for all substrings\n",
        "        for subString in subHamiltonian:\n",
        "            meas[subString] = 1\n",
        "        \n",
        "        # Go through the measurements on all the qubits\n",
        "        for i in range(len(qubits)):\n",
        "            \n",
        "            if (mainString[i] != \"I\"):\n",
        "                # There's a measurement associated with this qubit\n",
        "                \n",
        "                # Use this single qubit measurement for the calculation of the\n",
        "                #measurement of each full substring in this repetition. If the\n",
        "                #substring has a \"0\" in the position corresponding to this\n",
        "                #qubit, the operator associated is I, and the measurement\n",
        "                #is ignored (raised to the power of 0)\n",
        "                for subString in subHamiltonian:\n",
        "                    meas[subString] = meas[subString]*((1-2*results.data[str(i)][j])\\\n",
        "                                                     **int(subString[i]))\n",
        "                        \n",
        "        # Add this measurement to the total, for each string\n",
        "        for subString in subHamiltonian:\n",
        "            total[subString]+=meas[subString]\n",
        "        \n",
        "    totalExpectationValue = 0\n",
        "    \n",
        "    # Calculate the expectation value of the subHamiltonian, by multiplying\n",
        "    #the expectation value of each substring by the respective coefficient\n",
        "    for subString in subHamiltonian:\n",
        "        \n",
        "        # Get the expectation value of this substring by taking the average\n",
        "        #over all the repetitions\n",
        "        expectationValue = total[subString]/repetitions\n",
        "        \n",
        "        # Add this value to the total expectation value, weighed by its \n",
        "        #coefficient\n",
        "        totalExpectationValue+=expectationValue*subHamiltonian[subString]\n",
        "    \n",
        "    return(totalExpectationValue)\n",
        "\n",
        "count = 0\n",
        "trackOptimization = False\n",
        "\n",
        "def stateEnergy(statePreparationGates,hamiltonian,repetitions,qubits):\n",
        "    ''' \n",
        "    Returns: the experimental energy expectation in a given state.\n",
        "    Arguments: a list of 6 parameters defining the state in which to obtain the \n",
        "    expectation value, a dictionary with the coefficients of each Pauli string\n",
        "    term in the Hamiltonian, the number of repetitions to be used to calculate \n",
        "    the expectation value of each string.\n",
        "    ''' \n",
        "    \n",
        "    # Print the percentage of the maximum number of function evaluations that \n",
        "    #has been used so far in the classical optimization, if it's a multiple of\n",
        "    #10% (to inform on the stage of the optimization process).\n",
        "    global count, trackOptimization\n",
        "    if(trackOptimization):\n",
        "        if (count%(300/10) == 0):\n",
        "            print(round(count/(300/100)),\"%\",sep = '')\n",
        "        count = count+1\n",
        "\n",
        "    groupedHamiltonian = groupHamiltonian(hamiltonian)\n",
        "\n",
        "    experimentalEnergyExpectation = 0\n",
        "    \n",
        "    # Obtain experimental expectation value for each necessary Pauli string by\n",
        "    #calling the measureExpectation function, and perform the necessary weighed\n",
        "    #sum to obtain the energy expectation value.\n",
        "    for mainString in groupedHamiltonian:\n",
        "         \n",
        "        expectationValue = measureExpectation(mainString,\\\n",
        "                                            groupedHamiltonian[mainString],\\\n",
        "                                                repetitions,statePreparationGates,qubits)\n",
        "        experimentalEnergyExpectation+=expectationValue\n",
        "\n",
        "    return experimentalEnergyExpectation\n",
        "\n",
        "def exactStateEnergy(stateCoordinates,hamiltonian):\n",
        "    ''' \n",
        "    Calculates the exact energy in a specific state.\n",
        "\n",
        "    Arguments:\n",
        "      stateCoordinates (np.ndarray): the state in which to obtain the \n",
        "        expectation value.\n",
        "      hamiltonian (dict): the Hamiltonian of the system.\n",
        "    \n",
        "    Returns:\n",
        "      exactEnergy (float): the energy expecation value in the state.\n",
        "    ''' \n",
        "\n",
        "    exactEnergy = 0\n",
        "    \n",
        "    # Obtain the theoretical expectation value for each Pauli string in the\n",
        "    #Hamiltonian by matrix multiplication, and perform the necessary weighed\n",
        "    #sum to obtain the energy expectation value.\n",
        "    for pauliString in hamiltonian:\n",
        "        \n",
        "        ket = np.array(stateCoordinates,dtype = complex)\n",
        "        bra = np.conj(ket)\n",
        "        \n",
        "        ket = np.matmul(stringToMatrix(pauliString),ket)\n",
        "        expectationValue = np.real(np.dot(bra,ket))\n",
        "        \n",
        "        exactEnergy+=\\\n",
        "            hamiltonian[pauliString]*expectationValue\n",
        "            \n",
        "    return exactEnergy\n",
        "    \n",
        "def exactStateEnergySparse(stateVector,sparseHamiltonian):\n",
        "    ''' \n",
        "    Calculates the exact energy in a specific state, using sparse matrices.\n",
        "\n",
        "    Arguments:\n",
        "      stateVector (Union[np.ndarray, scipy.sparse.csc_matrix): the state in \n",
        "        which to obtain the expectation value.\n",
        "      sparseHamiltonian (scipy.sparse.csc_matrix): the Hamiltonian of the system.\n",
        "    \n",
        "    Returns:\n",
        "      energy (float): the energy expecation value in the state.\n",
        "    ''' \n",
        "\n",
        "    if not isinstance(stateVector,scipy.sparse.csc_matrix):\n",
        "      ket = scipy.sparse.csc_matrix(stateVector,dtype=complex).transpose()\n",
        "    else:\n",
        "      ket = stateVector\n",
        "      \n",
        "    bra = ket.transpose().conj()\n",
        "\n",
        "    energy = (bra * sparseHamiltonian * ket)[0,0].real\n",
        "    \n",
        "    return energy\n",
        "    \n",
        "def energyFromExpOperator(operator,referenceState,sparseHamiltonian):\n",
        "  '''\n",
        "  Calculates the energy in the state obtained by applying e^operator to the\n",
        "  reference state, using sparse matrices.\n",
        "\n",
        "  Arguments:\n",
        "    operator (openfermion.FermionOperator): the operator to be exponentiated\n",
        "      and applied to the reference state.\n",
        "    referenceState (np.ndarray): the reference state.\n",
        "    sparseHamiltonian (scipy.sparse.csc_matrix): the Hamiltonian of the system,\n",
        "      in matrix form.\n",
        "\n",
        "  Returns:\n",
        "    energy (float): the energy.\n",
        "  '''\n",
        "  # Obtain the dimension and, from it, calculate the number of qubits\n",
        "  dimension, _ = sparseHamiltonian.shape\n",
        "  qubitNumber = int(np.log(dimension)/np.log(2)) \n",
        "  \n",
        "  # Obtain the sparse matrix representing the operator, and act on the state\n",
        "  #with its exponentiated version\n",
        "  sparseOperator = get_sparse_operator(operator,qubitNumber)\n",
        "  ket = scipy.sparse.csc_matrix(referenceState,dtype=complex).transpose()\n",
        "  ket = scipy.sparse.linalg.expm_multiply(sparseOperator,ket)\n",
        "\n",
        "  # Calculate the energy expectation value\n",
        "  bra = ket.transpose().conj()\n",
        "  energy = (bra * sparseHamiltonian * ket)[0,0].real\n",
        "\n",
        "  return energy"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tdSEQd6RUKt7"
      },
      "source": [
        "## Molecule Definitions "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f4MfKM6OMC6s"
      },
      "source": [
        "# H2\n",
        "geometry = [['H',[0,0,0]],['H',[0,0,0.74]]]\n",
        "basis = 'sto-3g'\n",
        "multiplicity = 1\n",
        "charge = 0\n",
        "h2molecule = MolecularData(geometry,basis,multiplicity,charge)\n",
        "h2molecule = run_pyscf(h2molecule,run_fci = True,run_ccsd = True)\n",
        "\n",
        "# HeH+\n",
        "r = 1 # interatomic distance in angstrom\n",
        "geometry = [['He',[0,0,0]],['H',[0,0,r]]]\n",
        "basis = 'sto-3g'\n",
        "multiplicity = 1\n",
        "charge = +1\n",
        "helonium = MolecularData(geometry,basis,multiplicity,charge)\n",
        "helonium = run_pyscf(helonium,run_fci = True,run_ccsd = True)\n",
        "\n",
        "# Get the PySCF molecule as well (it's not necessary as there is the OF-PySCF\n",
        "#plugin, just for confirmation)\n",
        "helonium_pysf = gto.M(atom = geometry, charge = charge, basis = basis,verbose = 0)\n",
        "\n",
        "# HeH+\n",
        "r = 1 # interatomic distance in angstrom\n",
        "geometry = [['He',[0,0,0]],['H',[0,0,r]]]\n",
        "basis = 'sto-3g'\n",
        "multiplicity = 1\n",
        "charge = +1\n",
        "helonium = MolecularData(geometry,basis,multiplicity,charge)\n",
        "helonium = run_pyscf(helonium,run_fci = True,run_ccsd = True)\n",
        "\n",
        "# LiH\n",
        "bondLength = 1.45 # interatomic distance in angstrom\n",
        "geometry = [['Li',[0,0,0]],['H',[0,0,bondLength]]]\n",
        "basis = 'sto-3g'\n",
        "multiplicity = 1\n",
        "charge = 0\n",
        "liH = MolecularData(geometry,basis,multiplicity,charge)\n",
        "liH = run_pyscf(liH,run_fci = True,run_ccsd = True)\n",
        "\n",
        "# Alternative to using run_pyscf: load from OpenFermion (data for this \n",
        "#particular molecule at this particule interatomic distance is available in \n",
        "#a file that comes with OF)\n",
        "liHOF = MolecularData(geometry,basis,multiplicity,charge,description = '1.45')\n",
        "liHOF.load()"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mTsz50AcUelJ"
      },
      "source": [
        "## Choosing the molecule"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gBO3PXeH5hYs"
      },
      "source": [
        "molecule = helonium"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vCyKDV6u6C8w"
      },
      "source": [
        "qubitNumber = molecule.n_qubits\n",
        "electronNumber = molecule.n_electrons\n",
        "orbitalNumber = molecule.n_orbitals"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NYy8gsGFy4LA"
      },
      "source": [
        "## Getting some data on the molecule"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o5Q9-f5c5urH"
      },
      "source": [
        "### Hamiltonian"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OLy_HUj53YnZ"
      },
      "source": [
        "The Hamiltonian of the molecule (in terms of fermionic raising/creation and lowering/annihilation operators) can be obtained using **get_molecular_hamiltonian**, a method of the **MolecularData** class. \n",
        "\n",
        "When no arguments are supplied, the method considers all orbitals active. Since the *STO-3G* minimal basis was chosen in setting the molecule, *all orbitals* means all the ones considered in a minimal basis: those that belong to a (fully or partially) occupied shell. \n",
        "\n",
        ">*Example*: In helonium, for instance, this is equivalent to setting >*active_indices = [0,1]*. Indices are of *spatial* orbitals. \n",
        "\n",
        "Orbital that belong to the same shell have the same principal quantum number. So if we have $_1H: 1s^1$ or $_2He: 1s^2$ we use one orbital ($1s$); if we have $_3Li:1s^2 1s^1$ we use 5 ($1s$, $2s$, $2p_x$, $2p_y$, $2p_z$).\n",
        "\n",
        "*STO* in *STO-3G* stands for *Slater type orbitals*, which means that we're using as molecular orbitals (wavefunctions of a many-particle system) anti-symmetrized products of atomic orbitals (wavefunctions of the individual particles). *3G* means thate we're approximating each Slater determinant by three Gaussian orbitals.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Fmbtkcf53ZLi"
      },
      "source": [
        "hamiltonian = molecule.get_molecular_hamiltonian() "
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XGYOqRd0ybGb"
      },
      "source": [
        "The Hamiltonian is returned as an instance of **InteractionOperator**. This is a special class for storing fermionic operators that:\n",
        "\n",
        "- Consist only of one-body and two-body terms\n",
        "- Consists only of terms that preserve particle number and spin\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4UfTa8px0Xlz",
        "outputId": "1ef18954-decc-4d8c-ddee-14cc46816665"
      },
      "source": [
        "print(type(hamiltonian))"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "<class 'openfermion.ops.representations.interaction_operator.InteractionOperator'>\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5eMkgaCl9LpU"
      },
      "source": [
        "The data structure consists of:\n",
        "- A constant term (attribute **constant**)\n",
        "- An NxN array of floats (matrix) representing the one body terms (attribute **one_body_tensor**)\n",
        "- An NxNxNxN array of floats representing the two body terms (attribute **two_body_tensor**)\n",
        "\n",
        "Where N is the number of qubits necessary to represent the molecule, that corresponds to the number of fermionic modes, or spin orbitals, that we choose to consider active.\n",
        "\n",
        "As an example, if term $h_{p, q, r, s} a^\\dagger_p a^\\dagger_q a_r a_s$ occurs in *hamiltonian*, then *hamiltonian.two_body_tensor* [p,q,r,s] = $h_{p,q,r,s}$.\n",
        "\n",
        "Since each term is particle number conserving, it must have as many creation operators as annihilation operators. So it's enough to decide a specific order for them to make the tensors non-ambiguous. The convention is *normal ordering* them. This corresponts to putting all creation operators to the left."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vZaqkVpM9NH5",
        "outputId": "8c79dc16-43a6-4711-fc73-b890f276213f"
      },
      "source": [
        "print(\"Constant:\\n\",hamiltonian.constant)\n",
        "print(\"Shape of One Body Tensor:\\n\",hamiltonian.one_body_tensor.shape)\n",
        "print(\"Shape of Two Body Tensor:\\n\",hamiltonian.two_body_tensor.shape)\n",
        "print(\"Hamiltonian:\")\n",
        "print(hamiltonian)"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Constant:\n",
            " 1.05835442184\n",
            "Shape of One Body Tensor:\n",
            " (4, 4)\n",
            "Shape of Two Body Tensor:\n",
            " (4, 4, 4, 4)\n",
            "Hamiltonian:\n",
            "() 1.05835442184\n",
            "((0, 1), (0, 0)) -2.4305217022828143\n",
            "((0, 1), (2, 0)) 0.17523764451269103\n",
            "((1, 1), (1, 0)) -2.4305217022828143\n",
            "((1, 1), (3, 0)) 0.17523764451269103\n",
            "((2, 1), (0, 0)) 0.1752376445126909\n",
            "((2, 1), (2, 0)) -1.3235920553067362\n",
            "((3, 1), (1, 0)) 0.1752376445126909\n",
            "((3, 1), (3, 0)) -1.3235920553067362\n",
            "((0, 1), (0, 1), (0, 0), (0, 0)) 0.4748839522005155\n",
            "((0, 1), (0, 1), (0, 0), (2, 0)) -0.08761816148293655\n",
            "((0, 1), (0, 1), (2, 0), (0, 0)) -0.08761816148293655\n",
            "((0, 1), (0, 1), (2, 0), (2, 0)) 0.05857335548569318\n",
            "((0, 1), (1, 1), (1, 0), (0, 0)) 0.4748839522005155\n",
            "((0, 1), (1, 1), (1, 0), (2, 0)) -0.08761816148293655\n",
            "((0, 1), (1, 1), (3, 0), (0, 0)) -0.08761816148293655\n",
            "((0, 1), (1, 1), (3, 0), (2, 0)) 0.05857335548569318\n",
            "((0, 1), (2, 1), (0, 0), (0, 0)) -0.08761816148293655\n",
            "((0, 1), (2, 1), (0, 0), (2, 0)) 0.05857335548569318\n",
            "((0, 1), (2, 1), (2, 0), (0, 0)) 0.28505319223565884\n",
            "((0, 1), (2, 1), (2, 0), (2, 0)) 0.031643292748480115\n",
            "((0, 1), (3, 1), (1, 0), (0, 0)) -0.08761816148293655\n",
            "((0, 1), (3, 1), (1, 0), (2, 0)) 0.05857335548569318\n",
            "((0, 1), (3, 1), (3, 0), (0, 0)) 0.28505319223565884\n",
            "((0, 1), (3, 1), (3, 0), (2, 0)) 0.031643292748480115\n",
            "((1, 1), (0, 1), (0, 0), (1, 0)) 0.4748839522005155\n",
            "((1, 1), (0, 1), (0, 0), (3, 0)) -0.08761816148293655\n",
            "((1, 1), (0, 1), (2, 0), (1, 0)) -0.08761816148293655\n",
            "((1, 1), (0, 1), (2, 0), (3, 0)) 0.05857335548569318\n",
            "((1, 1), (1, 1), (1, 0), (1, 0)) 0.4748839522005155\n",
            "((1, 1), (1, 1), (1, 0), (3, 0)) -0.08761816148293655\n",
            "((1, 1), (1, 1), (3, 0), (1, 0)) -0.08761816148293655\n",
            "((1, 1), (1, 1), (3, 0), (3, 0)) 0.05857335548569318\n",
            "((1, 1), (2, 1), (0, 0), (1, 0)) -0.08761816148293655\n",
            "((1, 1), (2, 1), (0, 0), (3, 0)) 0.05857335548569318\n",
            "((1, 1), (2, 1), (2, 0), (1, 0)) 0.28505319223565884\n",
            "((1, 1), (2, 1), (2, 0), (3, 0)) 0.031643292748480115\n",
            "((1, 1), (3, 1), (1, 0), (1, 0)) -0.08761816148293655\n",
            "((1, 1), (3, 1), (1, 0), (3, 0)) 0.05857335548569318\n",
            "((1, 1), (3, 1), (3, 0), (1, 0)) 0.28505319223565884\n",
            "((1, 1), (3, 1), (3, 0), (3, 0)) 0.031643292748480115\n",
            "((2, 1), (0, 1), (0, 0), (0, 0)) -0.08761816148293655\n",
            "((2, 1), (0, 1), (0, 0), (2, 0)) 0.28505319223565884\n",
            "((2, 1), (0, 1), (2, 0), (0, 0)) 0.05857335548569318\n",
            "((2, 1), (0, 1), (2, 0), (2, 0)) 0.03164329274848008\n",
            "((2, 1), (1, 1), (1, 0), (0, 0)) -0.08761816148293655\n",
            "((2, 1), (1, 1), (1, 0), (2, 0)) 0.28505319223565884\n",
            "((2, 1), (1, 1), (3, 0), (0, 0)) 0.05857335548569318\n",
            "((2, 1), (1, 1), (3, 0), (2, 0)) 0.03164329274848008\n",
            "((2, 1), (2, 1), (0, 0), (0, 0)) 0.05857335548569318\n",
            "((2, 1), (2, 1), (0, 0), (2, 0)) 0.03164329274848008\n",
            "((2, 1), (2, 1), (2, 0), (0, 0)) 0.031643292748480115\n",
            "((2, 1), (2, 1), (2, 0), (2, 0)) 0.37311534589993395\n",
            "((2, 1), (3, 1), (1, 0), (0, 0)) 0.05857335548569318\n",
            "((2, 1), (3, 1), (1, 0), (2, 0)) 0.03164329274848008\n",
            "((2, 1), (3, 1), (3, 0), (0, 0)) 0.031643292748480115\n",
            "((2, 1), (3, 1), (3, 0), (2, 0)) 0.37311534589993395\n",
            "((3, 1), (0, 1), (0, 0), (1, 0)) -0.08761816148293655\n",
            "((3, 1), (0, 1), (0, 0), (3, 0)) 0.28505319223565884\n",
            "((3, 1), (0, 1), (2, 0), (1, 0)) 0.05857335548569318\n",
            "((3, 1), (0, 1), (2, 0), (3, 0)) 0.03164329274848008\n",
            "((3, 1), (1, 1), (1, 0), (1, 0)) -0.08761816148293655\n",
            "((3, 1), (1, 1), (1, 0), (3, 0)) 0.28505319223565884\n",
            "((3, 1), (1, 1), (3, 0), (1, 0)) 0.05857335548569318\n",
            "((3, 1), (1, 1), (3, 0), (3, 0)) 0.03164329274848008\n",
            "((3, 1), (2, 1), (0, 0), (1, 0)) 0.05857335548569318\n",
            "((3, 1), (2, 1), (0, 0), (3, 0)) 0.03164329274848008\n",
            "((3, 1), (2, 1), (2, 0), (1, 0)) 0.031643292748480115\n",
            "((3, 1), (2, 1), (2, 0), (3, 0)) 0.37311534589993395\n",
            "((3, 1), (3, 1), (1, 0), (1, 0)) 0.05857335548569318\n",
            "((3, 1), (3, 1), (1, 0), (3, 0)) 0.03164329274848008\n",
            "((3, 1), (3, 1), (3, 0), (1, 0)) 0.031643292748480115\n",
            "((3, 1), (3, 1), (3, 0), (3, 0)) 0.37311534589993395\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jtxr-yXu6z8_"
      },
      "source": [
        "It is necessary to represent the problem in a way that allows solving it with a quantum computer. Fermionic states will have to be mapped into qubit states, and fermionic operators will have to be mapped into qubit operators.\n",
        "\n",
        "For example, it is necessary to transform the Hamiltonian into one that acts on the qubits that represent the fermionic state.\n",
        "\n",
        "This is achieved by performing the **Jordan Wigner** transformation, which maps a **FermionOperator** into a **QubitOperator**. Since the function only accepts FermionOperators, the Hamiltonian must be transformed into one first. This is always possible, because instances of **InteractionOperator** are fermion operators with special characteristics.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VkmZSzAh-ZES"
      },
      "source": [
        "qubitHamiltonian = jordan_wigner(get_fermion_operator(hamiltonian))\n",
        "\n",
        "# Eliminate terms that are negligeable up to a certain tolerance\n",
        "qubitHamiltonian.compress()"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SsyCFU3a-bOf"
      },
      "source": [
        "If we pass the qubit Hamiltonian as an argument to **count_qubits**, it will return the minimum number of qubits it acts on.  Under the Jordan Wigner transformation, each qubit represents a specific fermionic mode (spin orbital). \n",
        "\n",
        "The number will be the same as **molecule.n_qubits**: the number of qubits necessary to represent a state of the molecule.\n",
        "\n",
        ">*Example*: \n",
        "In $HeH^+$, there are 4 spin orbitals: two $1s$ spatial orbitals from each atom, each with two associated spin orbitals ($\\alpha$ and $\\beta$ spins). \n",
        ">\n",
        ">In order to represent a state of this molecule, we need four qubits: the state of each one will represent the occupation number of the corresponding spin orbital (0 or 1, given the Pauli exclusion principle).\n",
        ">\n",
        ">A qubit operator acting on the state has to live in the same dimension or in a lower one. As for the Hamiltonian, it necessarily involves all of the particles representing the molecule: if they didn't affect the energy, they wouldn't be in the molecule."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xZDj6zEU50IP",
        "outputId": "60f435eb-40b5-43d4-b277-f547b9c7d988"
      },
      "source": [
        "# Compute the minimum number of qubits on which operator acts; it's at most \n",
        "#molecule.n_qubits\n",
        "print(\"Qubit Number:\",count_qubits(qubitHamiltonian))\n",
        "\n",
        "# Get the number of spatial orbitals of the molecule. \n",
        "# Since for each spacial orbital there are two spin orbitals (alpha, 'up', or\n",
        "#beta, 'down', spins), this is half the number of qubits.\n",
        "print(\"Orbital Number:\",molecule.n_orbitals)\n",
        "\n",
        "# Get the energies associated to the orbitals.\n",
        "# Lower index -> lower energy.\n",
        "print(\"Orbital Energies:\",molecule.orbital_energies)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Qubit Number: 4\n",
            "Orbital Number: 2\n",
            "Orbital Energies: [-1.48075592 -0.30052523]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jGKUgCuL51eB"
      },
      "source": [
        "### Hartree Fock State"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "RYmZZzxN2XS7",
        "outputId": "e8417c9c-6348-43e6-f062-84f5fbbb51b7"
      },
      "source": [
        "# Get the Hartree Fock state, assuming the Jordan Wigner mapping.\n",
        "# Second argument corresponds to SPIN orbitals, hence the factor of two.\n",
        "hfState = jw_hartree_fock_state(molecule.n_electrons,molecule.n_orbitals*2) \n",
        "\n",
        "print(\"Hartree Fock Vector:\",hfState)\n",
        "\n",
        "# The vector corresponds to a computational basis state with the qubits up to\n",
        "#molecule.n_electrons in state |1>, and all the other ones in state |0> \n",
        "#(|11...00>).\n",
        "# Corresponds to a single Slater Determinant, which includes the lowest energy\n",
        "#atomic orbitals. \n",
        "# By acting on this state with excitation operators, the UCCSD ansatz will \n",
        "#create a superposition of Slater Determinants that may correspond to an \n",
        "#entangled state.\n",
        "print(\"Hartree Fock Ket:\",fromVectortoKet(hfState))\n",
        "\n",
        "# So the Hartree Fock state preparation circuit is implemented on CIRQ as:\n",
        "qubits = cirq.LineQubit.range(molecule.n_qubits)\n",
        "\n",
        "hfGates = [cirq.X(qubits[i]) for i in range(molecule.n_electrons)]\n",
        "\n",
        "# Which is equivalent to:\n",
        "hfGates = [cirq.X(qubits[i]) for i in range(molecule.n_qubits) \\\n",
        "  if fromVectortoKet(hfState)[i]] "
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Hartree Fock Vector: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]\n",
            "Hartree Fock Ket: [1, 1, 0, 0]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hkG7GXU85_65"
      },
      "source": [
        "### Ground Energy"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nCss7V3oqU0B",
        "outputId": "a4b3f19a-a27c-4997-ad14-ff68287049fe"
      },
      "source": [
        "# Get the sparse matrix representing the Hamiltonian in the qubit space\n",
        "sparseHamiltonian = get_sparse_operator(hamiltonian)\n",
        "\n",
        "# Reshape the qubit Hamiltonian into a dictionary associating to Pauli strings\n",
        "#(e.g. 'XXYY') their coefficients (as suitable for use in \n",
        "#stateEnergy)\n",
        "formattedHamiltonian = convertHamiltonian(qubitHamiltonian)\n",
        "\n",
        "# Group strings that only differ by identity terms\n",
        "groupedHamiltonian = groupHamiltonian(formattedHamiltonian)\n",
        "\n",
        "# One can get the Hamiltonian's matrix form by adding the matrices corresponding\n",
        "#to each of the Pauli strings, as \n",
        "#hamiltonianMatrix = hamiltonianToMatrix(hamiltonian)\n",
        "#However, it's really inefficient; it's better to use the sparse matrix from\n",
        "#OpenFermion\n",
        "\n",
        "qubits = cirq.LineQubit.range(molecule.n_qubits)\n",
        "shots = 1000\n",
        "\n",
        "# Measure HF energy by getting the statistics on the CIRQ simulator\n",
        "myenergy = stateEnergy(hfGates,formattedHamiltonian,shots,qubits)\n",
        "\n",
        "# Calculate exact HF energy using matrix algebra\n",
        "calculatedEnergy = exactStateEnergy(hfState,formattedHamiltonian) \n",
        "\n",
        "# For greater efficiency, it's better to use sparse matrices:\n",
        "calculatedEnergySparse = exactStateEnergySparse(hfState,sparseHamiltonian)\n",
        "\n",
        "tol = 10**-10\n",
        "assert(calculatedEnergy-calculatedEnergySparse < tol)\n",
        "\n",
        "# Get the energy from PySCF\n",
        "hf_pyscf = scf.RHF(helonium_pysf)\n",
        "hfEnergy_pyscf = hf_pyscf.kernel()\n",
        "\n",
        "print(\"HF\\nMy measured energy: {}\".format(myenergy))\n",
        "\n",
        "print(\"My calculated energy: {}\".format(calculatedEnergy))\n",
        "\n",
        "print(\"Openfermion energy (attribute): {}\".format(molecule.hf_energy))\n",
        "\n",
        "print(\"Openfermion energy (calculated): {}\".format(np.real\\\n",
        "      (openfermion.expectation(get_sparse_operator(qubitHamiltonian),hfState))))\n",
        "\n",
        "print(\"PySCF energy: {}\".format(hfEnergy_pyscf))\n",
        "\n",
        "error = molecule.hf_energy - molecule.fci_energy\n",
        "print(\"\\nError:\",error)\n",
        "print(\"(in % of chemical accuracy: {:.3f}%)\\n\".format(error/chemicalAccuracy*100))"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HF\n",
            "My measured energy: -2.853916883515915\n",
            "My calculated energy: -2.852921078324599\n",
            "Openfermion energy (attribute): -2.852921078324597\n",
            "Openfermion energy (calculated): -2.852921078324599\n",
            "PySCF energy: -2.852921078324597\n",
            "\n",
            "Error: 0.0072840442552060125\n",
            "(in % of chemical accuracy: 457.081%)\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HEq5GTF0FHt0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6cd58094-76b4-4afe-a8e8-224199ef764b"
      },
      "source": [
        "# Print the ground energy obtained from the classical variational full \n",
        "#configuration interaction approach (exact within the STO-3G basis set)\n",
        "print(\"FCI Energy (OpenFermion): {}\".format(molecule.fci_energy))\n",
        "\n",
        "# Get the same energy from PySCF\n",
        "cisolver = fci.FCI(helonium_pysf, hf_pyscf.mo_coeff)\n",
        "fciEnergy_pysf = cisolver.kernel()[0]\n",
        "print(\"FCI Energy (pyscf):\",fciEnergy_pysf)\n",
        "\n",
        "# Print the ground energy obtained from the classical variational coupled \n",
        "#cluster singles and doubles approach\n",
        "print(\"CCSD Energy: {}\\n\".format(molecule.ccsd_energy))\n",
        "\n",
        "# Get the ground state and energy of the Hamiltonian using OpenFermion's\n",
        "#get_ground_state function, which resorts scipy.sparse.linalg.eigsh\n",
        "# It's significantly lower than FCI.\n",
        "energyFromSparseOperator, stateFromSparseOperator = \\\n",
        "  get_ground_state(sparseHamiltonian)\n",
        "\n",
        "print(\"Energy from sparse operator: {}\".format(energyFromSparseOperator))\n",
        "\n",
        "calculatedEnergy = exactStateEnergySparse(stateFromSparseOperator,sparseHamiltonian) \n",
        "\n",
        "print(\"Same state, my calculated energy: {}\".format(calculatedEnergy))\n",
        "\n",
        "# Diagonalize Hamiltonian to get the exact ground state and energy, using my\n",
        "#own function. Same energy.\n",
        "formattedHamiltonian = convertHamiltonian(qubitHamiltonian)\n",
        "\n",
        "(exactGroundEnergies,exactGroundStates) = groundStatesFromDiagonalization\\\n",
        "  ([formattedHamiltonian])\n",
        "\n",
        "print(\"My energy from exact diagonalization: {}\".format\\\n",
        "      (np.real(exactGroundEnergies[0])))\n",
        "\n",
        "exactState = exactGroundStates[0]\n",
        "\n",
        "# Calculate overlap of the exact ground state obtained from exact \n",
        "#diagonalization with the hartree-fock state. \n",
        "print(\"Overlap with HF state:\",calculateOverlap(exactState,hfState))\n",
        "\n",
        "# Get eigenvalues of the Hamiltonian. The FCI energy is there, but there's a \n",
        "#lower value. \n",
        "# There shouldn't be? The Slater determinant with the lowest energy isn't the\n",
        "#one corresponding to the Hartree Fock solution, rather one one the wrong \n",
        "#number of electrons.\n",
        "print(\"Eigenspectrum of the Hamiltonian:\\n\",eigenspectrum(qubitHamiltonian))\n",
        "\n",
        "# The lower eigenvalue can't correspond to a physical state of the \n",
        "#system. The ground states obtained from get_ground_state or from my \n",
        "#diagonalization, that supposedly have lower energy than the FCI state, have \n",
        "#nonzero coefficients for e.g. |0001>. The real system should have always\n",
        "#2 electrons. FCI operators are all particle number conserving, they can't\n",
        "#reach a state with the wrong number of electrons starting from the Hartree Fock\n",
        "#state, and that's a good thing. To get the ground state from the matrix \n",
        "#of the Hamiltonian, it's necessary to first restrict it to live in the \n",
        "#space spanned by computational basis states that have the correct number of \n",
        "#electrons. See  openfermion.linalg.sparse_tools.jw_number_restrict_operator\n",
        "\n",
        "# Get all possible combinations of qubits to hold the correct number of \n",
        "#electrons. This is done by selecting electronNumber qubit indices, with no \n",
        "#repetition. 1 qubit in state |1> corresponds to 1 electron, because under \n",
        "#the Jordan Wigner mapping it represents the occupation of a specific spin\n",
        "#orbital. For HeH+ we have two electrons, so we get all the different pairs\n",
        "#of indices.\n",
        "occupations = itertools.combinations(range(qubitNumber), electronNumber)\n",
        "\n",
        "# Transform the selected qubit labels given by 'occupations' into indices.\n",
        "# Example: the Hartree Fock state |1100> will have index 2^0 + 2^1 = 3.\n",
        "select_indices = [sum([2**n for n in occupation]) for occupation in occupations]\n",
        "\n",
        "# Construct a restricted version of the Hamiltonian, by selecting only those \n",
        "#lines and columns that correspond to computational basis states that have the \n",
        "#correct number of particles.\n",
        "restrictedHamiltonian = sparseHamiltonian[np.ix_(select_indices, select_indices)]\n",
        "\n",
        "energyFromSparseOperator, stateFromSparseOperator = \\\n",
        "  get_ground_state(restrictedHamiltonian)\n",
        "\n",
        "print(\"\\nEnergy from sparse operator, restricted to the correct particle number: \",\\\n",
        "      (energyFromSparseOperator))\n",
        "\n",
        "# The length of the resulting ground state is not the dimension of the Hilbert\n",
        "#space, but the dimension of the subspace spanned by the computational basis\n",
        "#states that correspond to Slater Determinants with the correct number of \n",
        "#electrons.\n",
        "# For HeH+ there are only 6 coefficients (full Hilbert space is of dimension 16),\n",
        "#corresponding to the 6 computational basis states that two and only two |1> \n",
        "#states).\n",
        "print(\"Length of restricted vector:\",len(stateFromSparseOperator))\n",
        "\n",
        "# None of this should be necessary and this only happens for HeH+. Also happens\n",
        "#using Qiskit and Pyscf. Something to do with the representation of the energy\n",
        "#in Pyscf?"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "FCI Energy (OpenFermion): -2.860205122579803\n",
            "FCI Energy (pyscf): -2.860205122579803\n",
            "CCSD Energy: -2.860205139888532\n",
            "\n",
            "Energy from sparse operator: -3.1578592138000037\n",
            "Same state, my calculated energy: -3.1578592138000006\n",
            "My energy from exact diagonalization: -3.1578592138000006\n",
            "Overlap with HF state: 0.0\n",
            "Eigenspectrum of the Hamiltonian:\n",
            " [-3.15785921 -3.15785921 -2.86020512 -2.70774238 -2.2456425  -2.2456425\n",
            " -2.24279966 -2.24279966 -2.24279966 -2.0978611  -1.39924664 -1.39924664\n",
            " -0.74596009 -0.23815828 -0.23815828  1.05835442]\n",
            "\n",
            "Energy from sparse operator, restricted to the correct particle number:  -2.8602051225798033\n",
            "Length of restricted vector: 6\n"
          ]
        }
      ]
    }
  ]
}