HEA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cirq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Molecular Encoding - SuppInfo 2\n",
    "\n",
    "Step 1 is to map our molecules onto qubits. We want to find the ground states of 3 molecules: \n",
    "1. $H_2$ has 2 $1s$ orbitals, so 4 spin-orbitals. \n",
    "2. $LiH$ has $1s, 2s, 2p_x$ (we assume zero-filling for $2p_y, 2p_z$) for $Li$ and $1s$ for $H$, so 8 spin-orbitals. \n",
    "3. Same as above, but now we have 2 $H$ atoms, so 10 spin-orbitals. \n",
    "\n",
    "In the STO-3G (fitting 3 gaussians to the Slater atomic orbitals), we can write second quantized Hamiltonian as:\n",
    "$$H = H_1 + H_2 = \\sum\n",
    "\n",
    "For $LiH$ and $BeH_2$ we consider **perfect filling** in the inner $1s$ orbitals. \n",
    "\n",
    "**Okay I don't understand any of this - I'll just use the $H_2$ representation given and learn this later.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.ops import QubitOperator\n",
    "h2_hub = (0.011280 * QubitOperator('Z0 Z1') + \n",
    "          0.397936 * QubitOperator('Z0') + \n",
    "          0.397936 * QubitOperator('Z1') + \n",
    "          0.180931 * QubitOperator('X0 X1')) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "θ_0_0_1 was not in param_set\n",
      "θ_1_0_1 was not in param_set\n",
      "The HEA ansatz: \n",
      "0           1\n",
      "│           │\n",
      "@───────────X\n",
      "│           │\n",
      "X───────────@\n",
      "│           │\n",
      "Rx(θ_0_0_2) Rx(θ_1_0_2)\n",
      "│           │\n",
      "Rz(θ_0_0_3) Rz(θ_1_0_3)\n",
      "│           │\n",
      "@───────────X\n",
      "│           │\n",
      "X───────────@\n",
      "│           │\n",
      "Rz(θ_0_1_1) Rz(θ_1_1_1)\n",
      "│           │\n",
      "Rx(θ_0_1_2) Rx(θ_1_1_2)\n",
      "│           │\n",
      "Rz(θ_0_1_3) Rz(θ_1_1_3)\n",
      "│           │\n"
     ]
    }
   ],
   "source": [
    "# Define objective function as Hamiltonian averaging \n",
    "from openfermioncirq import HamiltonianObjective\n",
    "obj = HamiltonianObjective(h2_hub)\n",
    "\n",
    "# Define ansatz \n",
    "from openfermioncirq import VariationalAnsatz \n",
    "from openfermioncirq.variational.letter_with_subscripts import LetterWithSubscripts\n",
    "class HEA(VariationalAnsatz):\n",
    "    def __init__(self, n, d=1):\n",
    "        \"\"\"\n",
    "        n: number of qubits\n",
    "        d: number of entangling + rotation layers \n",
    "        \"\"\"\n",
    "        self.d = d \n",
    "        self.n = n\n",
    "        super().__init__(None) \n",
    "    \n",
    "    def params(self): \n",
    "        \"\"\"The parameters of the ansatz. N(3d + 2) parameters \"\"\"\n",
    "        for q in range(self.n): # Loop over qubits\n",
    "            for i in range(self.d + 1): # Loop over 0,1,...,d\n",
    "                for j in range(1, 4): # Loop over 1,2,3 for ZXZ \n",
    "                    # 2 starting X then Z rotations on each qubit\n",
    "                    if i == 0 and j == 1: continue # No first Z rotation\n",
    "                    \n",
    "                    yield LetterWithSubscripts('θ', q, i, j) \n",
    "    \n",
    "    def param_bounds(self): \n",
    "        bounds = []\n",
    "        for param in self.params():\n",
    "            bounds.append((-2*np.pi, 2*np.pi))\n",
    "        return bounds\n",
    "    \n",
    "    def _generate_qubits(self):\n",
    "        \"\"\"Produce qubits that can be used by the ansatz\"\"\"\n",
    "        return cirq.LineQubit.range(self.n)\n",
    "    \n",
    "    def operations(self, qubits): \n",
    "        \"\"\"Produce the operations of the ansatz circuit\"\"\"\n",
    "        param_set = set(self.params()) \n",
    "        \n",
    "        for i in range(self.d + 1): \n",
    "            for q in range(self.n):\n",
    "                if q != self.n - 1: \n",
    "                    yield cirq.ops.CNOT.on(self.qubits[q], self.qubits[q+1])\n",
    "                    yield cirq.ops.CNOT.on(self.qubits[q+1], self.qubits[q])\n",
    "                for j in range(1,4):\n",
    "                    sym = LetterWithSubscripts('θ', q, i, j)\n",
    "                    if sym not in param_set: \n",
    "                        print('{} was not in param_set'.format(sym))\n",
    "                        continue \n",
    "                        \n",
    "                    if j == 1 or j == 3: \n",
    "                        yield cirq.rz(sym).on(self.qubits[q])\n",
    "                    elif j == 2: \n",
    "                        yield cirq.rx(sym).on(self.qubits[q])\n",
    "                                                \n",
    "                    else: print('j = ', j, '. But Should be 1,2,or 3.')\n",
    "                        \n",
    "# Initialize ansatz \n",
    "ansatz = HEA(2, 1)\n",
    "print('The HEA ansatz: ')\n",
    "print(ansatz.circuit.to_text_diagram(transpose=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created a study with 2 qubits and 10 parameters\n",
      "(4, 4)\n",
      "Value of objective with default initial params is -0.011280000000000012\n"
     ]
    }
   ],
   "source": [
    "# Create a variational study \n",
    "from openfermioncirq import VariationalStudy \n",
    "study = VariationalStudy('HEA', ansatz, obj, initial_state=np.array([0,1,0,0]))\n",
    "\n",
    "print('Created a study with {} qubits and {} parameters'.format(\n",
    "        len(study.ansatz.qubits), study.num_params))\n",
    "print(obj._hamiltonian_linear_op.shape)\n",
    "\n",
    "print('Value of objective with default initial params is {}'.format(\n",
    "        study.value_of(ansatz.default_initial_params())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.8048889791046269\n"
     ]
    }
   ],
   "source": [
    "from openfermioncirq.optimization import ScipyOptimizationAlgorithm, OptimizationParams\n",
    "\n",
    "# Optimize\n",
    "algorithm = ScipyOptimizationAlgorithm(\n",
    "    kwargs={'method': 'COBYLA'},\n",
    "#     options={'maxiter': 100},\n",
    "    uses_bounds=False)\n",
    "optimization_params = OptimizationParams(\n",
    "    algorithm=algorithm)\n",
    "result = study.optimize(optimization_params)\n",
    "print(result.optimal_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}