AdiabaticEvolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing out adiabatic evolution on $2 \\times 2$ matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg\n",
    "import scipy.sparse\n",
    "from scipy.optimize import minimize_scalar\n",
    "import matplotlib.pyplot as plt\n",
    "tol = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_A = np.array([[1, 0], [0, -1]])\n",
    "# H_B = np.array([[0, 1], [1, 0]]) # If this is H_B, there isn't a crossing\n",
    "H_B = np.array([[-5, 0], [0, 4]]) # If this is H_B, there is a crossing so ground state goes to excited state\n",
    "H = lambda s: (1-s)*H_A + s*H_B\n",
    "psi_A = np.array([0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0, -1])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(H_A, psi_A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's do 5 steps: H(0), H(0.25), H(0.5), H(0.75), H(1) at t=1 \n",
    "n = 50\n",
    "t = 1\n",
    "res = psi_A\n",
    "for i in range(n): \n",
    "    s = i / (n-1)\n",
    "    res = np.dot(scipy.linalg.expm(-1j * H(s) * t), res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        +0.j        , 0.92175127+0.38778164j])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        +0.j        , 3.68700508+1.55112654j])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(H_B, res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot lowest two eigenvalues for n steps of H(s)\n",
    "x = np.arange(0, 1 + 1e-10, 1/(n-1))\n",
    "y_ground = [] \n",
    "y_excited = [] \n",
    "for i in range(n): \n",
    "    s = i / (n-1)\n",
    "    eig, _ = np.linalg.eigh(H(s)) # Should already be in ascending order \n",
    "    assert len(eig) == 2 \n",
    "    y_ground.append(eig[0])\n",
    "    y_excited.append(eig[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEbCAYAAADTZlM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xUdfb/8ddJQq9SBUIINST0gFIUpYgoSBGQOlEXBCxA/LpV1y2/XXddd1fXoIsUddUMvUlVFLGgqyi9TOgQCCBVOunn98cMLpsFEsgkN5M5z8cjDzN35t45F3Decz/3fs4VVcUYY0zwCXG6AGOMMc6wADDGmCBlAWCMMUHKAsAYY4KUBYAxxgQpCwBjjAlSFgDGGBOkLACMMSZIWQCYYktEVEQuiMifrvbYgXpWiUiqiHx5A+sUWs1Xq+9majaBwwLAFBoR+Y2InM/xk+r7kBt6ldfvF5F78vm2rVT119d5nPM9P/PVdLm+HTfyZiJSRUQW+j60k0Vk+OXnVLUb8PiN78J/ar7e9m+wzsa+/XRfr7581GwCgAWAKTSq+kdVLX/5B7gV2AgsA+Y7W91/GXdFnVE3uO4/gXSgJjACeENEmvmxNn9t/5/Ad36sywQgCwDjCBEpAywFLgCDVDUjx/OJQASwxPdN/Be+5dG+b+mnRWSbiPQt9OKvQUTKAQOB36jqeVX9ElgMxBWl7fuOtk4Dn/ijLhO4LABMoRORksACoCTQV1VTc75GVeOAA0Af3zfxv4pICWAJ8BFQAxgPTBeRG/2WnpsXReSEiHwlIl1uYL0mQKaq7rxi2SbAX0cA+d6+iFQE/gA846eaTACzADCFSkTCgFl4hzDuV9ULN7B6B6A88BdVTVfVVXiPIob5scRfAg2AOsBUvEcgDfO4bnngbI5lZ4AKfqrNH9v/I/CWqqb4qSYTwCwATKERkRDgHSAKuFdVz9zgJmoDB1U1+4plyXg/rG+mnhFXnOz9AEBV16jqOVVNU9V3ga+AXnnc5HmgYo5lFYFzN1Ofv7cvIq2Be4B/+KkeE+DCnC7ABJXJQHvgLlU9kYfX57xZxWGgroiEXBECEcBOboKqTgem56EGyeMmdwJhItJYVXf5lrUCtt1MfQWw/S5AJHBARMB7RBEqIjGqGuunGk0AsSMAUyhE5BXgPuAeVT2Sx9WO4h2OuWwNcBH4hYiU8I3P98E7pOSPGiuLSE8RKS0iYSIyArgL+PCK17wjIu9cbX3fcNYC4A8iUk5E7gD6AYnXec9rbu9mtp/L9qYCDYHWvp/JeK/A6pmX9zfFjwWAKXAi0hz4P7yXfW7LMQ/gqG9o6GpeBJ73XfHzM1VNx/uBfz9wApgEPKyq2/1UagngBeC4b/vjgf45TrrWxTssdC1PAmWAY8BM4AlVvd439Ny2d6Pbv+b2VPWiqn5/+QfvkFKqqh6/gfc3xYjYLSFNcSUiqUAaMFFVf5Pz8U1sryTeq25a5rxsNY/rf4z3RPa3qto9L9u7kZr9Xd+1lpniwwLAGGOClA0BGWNMkLIAMMaYIGUBYIwxQSqg5gFUq1ZNIyMjnS7DGGMCyrp1606oavWcywMqACIjI1m7dq3TZRhjTEARkeSrLbchIGOMCVIWAMYYE6QsAIwxJkhZABhjTJCyADDGmCDleACISKiIbBCRpU7XYowxwcTxAADigSSnizDGmGDjaACISDjQG3jTyTqMMaao2rTpGD//+WccP37R79t2eiLYq8AvuM49TUVkDDAGICIiopDKMsYY56SknGPGjCTcbg9btpwgLCyEu+6qS58+eb09dd44FgAi8gBwTFXX+e7sdFWqOhXvnYxo166d9a42xhRLZ8+mMX/+TtzuJD799ACq0LFjbf75z+4MHhxFtWpl/f6eTh4B3AH0FZFeQGmgooi4VdXlYE3GGFNoMjKyWLFiP263h0WL9pCamkmjRpX53e864XLF0LBh5QJ9f8cCQFWfBZ4F8B0B/Mw+/I0xxZ2q8u233+N2e5g1azsnTlyiatUyjBrVHJcrhvbtayEihVKL0+cAjDEmKOzZcxq324Pb7WH37tOULh1Gv34Ncbli6NkzkhIlQgu9piIRAKr6GfCZw2UYY4xfnTx5iTlzdpCY6OHrrw8jAl27RvDss+0ZOLAJlSqVcrS+IhEAxhhTXKSmZrJ06R7c7iSWL99LRkY2zZtX46WX7mLYsKbUrVvR6RJ/ZAFgjDH5lJ2trF6dgtvtYe7cnZw5k0atWuWYMCGWuLgYWrasXmjj+jfCAsAYY26Sx3MCtzuJ6dM9HDhwjnLlSjBwYGPi4prRtWtdQkOLQrOFa7MAMMaYG3DkyHlmztyO2+1hw4ZjhIYK994byYsv3kW/fg0pV66k0yXmmQWAMcbk4vz5dBYu3IXbncTKlclkZyvt2tXk1Ve7MnRoU2rWLOd0iTfFAsAYY64iMzObTz5JJjHRw8KFu7h4MZN69Sry7LPtcbmiadq0qtMl5psFgDHG+KgqGzYcIzHRw8yZSRw9epHKlUvhcsXgcsVwxx11CAkpeidzb5YFgDEm6CUnn2H6dG/ztaSkU5QsGcoDDzTA5YqhV6/6lCpVPD8qi+deGWNMLk6fTmXevJ0kJnr44osUADp3Dmfy5LY89FATqlQp43CFBc8CwBgTNNLTs/jgg3243R6WLNlDWloWUVFVeOGFOxk+vCn16xds87WixgLAGFOsqSpff30Yt9vD7Nk7OHUqlerVyzB2bCvi4mJo27ZmkZykVRgsAIwxxdKuXT/82Hxt794zlCkTRv/+jYiLi+Gee+o50nytqLEAMMYUG8ePX2T27O243UmsWXMEEejWLYLf/rYjAwY0oUKFwJmkVRgsAIwxAe3SpQwWL96D2+3hww/3k5mZTatW1fnb3+5m2LCm1KlzzTvOBj0LAGNMwMnKyubzz1NITNzG/Pm7OHcunTp1yvPMM21xuWJo0aK60yUGBAsAY0zA2LLlOG63h+nTkzh06DwVKpRk0KAmjBgRTZcuRb/5WlFjAWCMKdIOHTrHzJnbSUz0sHnzccLCQrjvvkhefrkLffs2pEyZEk6XGLAsAIwxRc65c+ksWLATtzuJTz5JRhXat6/Fa691Y8iQplSvXtbpEosFCwBjTJGQkZHFxx8n43Z7eP/93Vy6lEmDBpX4zW864nLF0LjxLU6XWOxYABhjHKOqrF37PW53EjNnJnH8+CWqVCnNI480Iy4uho4dawftJK3CYAFgjCl0+/ad9jVfS2LHjlOUKhVKnz4NcbliuP/++pQsaZO0CoMFgDGmUJw6dYm5c3fidnv48stDANx1Vzg/+1k7Bg1qQuXKpR2uMPhYABhjCkxaWibLlu3F7U5i2bK9pKdnER1dhT//uTPDhzelXr1KTpcY1CwAjDF+lZ2tfPXVIdxuD3Pm7OD06TRq1izLU0+1xuWKoU2bGjauX0RYABhj/GL79pO43UlMn+5h//6zlC0bxoABTXC5ounevR5hYTZJq6ixADDG3LSjRy8wa9Z23G4Pa9ceJSRE6NGjHn/8453079+I8uWt+VpRZgFgjLkhFy9m8P77u3G7PXz00X6yspTY2Jq88koXhg2L5tZbyzldoskjCwBjTK6ysrL59NODJCZuY8GCXZw/n0FERAV+8YvbcbmiiYmp5nSJ5iZYABhjrkpV2bz5OImJHmbMSOLIkQtUqlSKIUOaEhcXQ+fO4YSE2MncQGYBYIz5Lykp55gxI4nERA9bt56gRIkQevVqgMsVzQMPNKR0afvYKC7sb9IYw5kzacyf752k9dlnB1GFjh1rM2nSPQweHEXVqmWcLtEUAAsAY4JURkYWH364H7fbw+LFe0hNzaRRo8r8/vedGDEihoYNKztdoilgFgDGBBFVZc2aI7jdHmbN2sHJk5eoWrUMo0Y1Jy6uGbfffqtN0goiFgDGBIE9e07jdntwuz3s3n2a0qXD6NfP23zt3nsjrflakLIAMKaYOnnyErNnb8ftTuLrrw8jAl27RvDcc+0ZMKAJlSqVcrpE4zDHAkBE6gLvATUBBaaqaoJT9RhTHKSmZrJkyR7cbg/Ll+8jMzOb5s2r8Ze/dGb48Gjq1q3odImmCHHyCCAT+KmqrheRCsA6EflYVT0O1mRMwMnOVlavTiEx0cO8eTs5cyaN2rXL8/TTsbhcMbRqVcPpEk0R5VgAqOoR4Ijv93MikgTUASwAjMkDj+cEiYkepk9P4uDBc5QvX4IBAxoTF9eMrl3rEhpqzdfM9RWJcwAiEgm0AdZc5bkxwBiAiIiIQq3LmKLmyJHzzJq1ncREDxs2HCM0VOjZM5KXXrqLvn0bUq6cNV8zeed4AIhIeWA+8LSqns35vKpOBaYCtGvXTgu5PGMcd/58OgsX7sLtTmLlymSys5V27WqSkNCNIUOiqFnTmq+Zm+NoAIhICbwf/tNVdYGTtRhTlGRmZrNyZTJut4eFC3dx8WIm9epV5Nln2+NyRdO0aVWnSzTFgJNXAQnwFpCkqq84VYcxRYWqsn79UdzuJGbOTOLo0YtUrlwKlyuGuLgYOnWqY83XjF85eQRwBxAHbBGRjb5lz6nqcgdrMqbQJSefYfr0JNxuD0lJpyhZMpTevRsQFxdDr171KVXK8ZFaU0w5eRXQl4B9nTFB6fTpVObN20lioocvvkgBoHPncKZMacugQU2oUsWar5mCZ18tjCkk6elZfPDBPhITt7F06V7S0rKIiqrCCy/cyfDhTalf35qvmcJlAWBMAVJVvv76MG63h9mzd3DqVCo1apTl8cdb4XLF0LZtTWu+ZhxjAWBMAdi164cfm6/t3XuGMmXC6N+/EXFxMfToEUlYmE3SMs6zADDGT44fv/hj87U1a44gAt271+N3v+vEgw82pkIFm6RlihYLAGPy4eLFDJYs2UNioocVK/aTmZlNy5bV+dvf7mbYsKbUqVPB6RKNuSYLAGNuUFZWNp99dhC328P8+bs4dy6dOnXK88wzbXG5YmjRorrTJRqTJxYAxuTRli3HSUz0MGNGEocOnadChZIMGtQElyuGu+8Ot+ZrJuBYABhzHYcOnWPGjCTc7iQ2bz5OWFgI990Xycsvd6FPn4aULVvC6RKNuWkWAMbkcO5cOgsWeCdprVp1AFVo374Wr7/encGDo6hevazTJRrjFxYAxgAZGVl8/LG3+dr77+/m0qVMGjaszG9/25ERI2Jo3PgWp0s0xu8sAEzQUlXWrv0etzuJWbO2c+zYRapUKc2jjzYjLq4ZHTrUsklaplizADBBZ9++077ma0ns2HGKUqVC6dOnIS5XDPffX5+SJUOdLtGYQmEBYILCqVOXmDt3J263hy+/PATA3XeH8/Of38bAgY2pXLm0wxUaU/gsAEyxlZaWybJle3G7k1i2bC/p6VlER1fhz3/uzIgR0UREVHS6RGMcZQFgipXsbOWrrw6RmOhh7twdnD6dRs2aZXnqqda4XDG0aVPDxvWN8QmKAJg9eztr137PuHFtqFevktPlmAKwfftJEhM9TJ+eRHLyWcqWDWPAgCa4XNF0717Pmq8ZcxVBEQBbtpzgH/9YxyuvrKN//0bEx8fSuXO4fRMMcEePXmDWrO0kJnpYt+4oISFCjx71eOGFO+nfvxHly1vzNWOuR1TV6RryrF27drp27dqbWvfgwbNMmrSRqVM3c+pUKq1b1yA+PpahQ5tSunRQ5GCxcOFCOosW7cHt9vDRR/vJylJiY2vickUzbFg0t95azukSjSlyRGSdqrb7n+XBEgCXXbyYwfTpSSQkrGPbtpNUr16GsWNb8cQTraldu7yfKjX+lJWVzapVB3C7PSxYsIvz5zOIiKjAiBExuFzRxMRUc7pEY4o0C4AcVJVPPz1IQsI6lizZQ2hoCIMHRxEfH8vtt9fyy3uYm6eqbNp0HLfb23ztyJELVKpUioceakJcXAx33hlOSIgN4RmTFxYA17Fnz2lef30Db721hXPn0mnfvhbx8bEMGtSEEiVsUlBhOnjwLDNmJJGY6GHbtpOUKBFCr14NcLmieeCBhjZcZ8xN8FsAiEgIUF5Vz/qruLwqqAC47Ny5dN55ZysTJ65n9+7T1K5dniefbM2YMS2tAVgBOnMmjfnzvZO0PvvsIKrQqVNtXK4YBg+OomrVMk6XaExAy1cAiMgM4HEgC/gOqAgkqOrf/F3o9RR0AFyWna188MFeEhLW8/HHyZQqFcqIEdHEx7elZUu72Yc/pKdnsWLFftxuD4sX7yE1NZNGjSrjcsXgcsXQsGFlp0s0ptjIbwBsVNXWIjICiAV+BaxT1Zb+L/XaCisAruTxnOC11zbw3nvbuHgxk7vvDic+vi19+za0G4DcIFVlzZojuN0eZs3awcmTl6hWrQxDhzbF5Yrh9ttvtUtzjSkA+Q2AbUBrYAbwuqp+LiKbVLWV/0u9NicC4LIffkjlzTc38/rrGzhw4ByRkRUZN64No0a1sD4yudi9+wdf8zUPu3efpnTpMPr2bUhcXAw9e0baeRZjClh+A2AC8EtgE9AbiADcqtrZ34Vej5MBcFlmZjaLFu0mIWE9q1enULZsGI880owJE2Jp2rSqo7UVJSdOXGTOnB0kJnr45psjiEDXrhG4XNEMGNCESpVKOV2iMUHD71cBiUiYqmbmu7IbUBQC4EobNhxl4sT1zJixnfT0LHr2jCQ+PpaePesH5SWKqamZLFninaS1fPk+MjOzad68GnFxMQwb1pS6da35mjFOyO8RQE3gz0BtVb1fRGKAjqr6lv9LvbaiFgCXHTt2gSlTNjNp0ka+//4CUVFVGD++DY880qzYtyPIzla++OIgbncSc+fu4OzZdGrXLs/w4d5x/VatajhdojFBL78B8AHwL+DXqtpKRMKADarawv+lXltRDYDL0tOzmDt3BwkJ6/nuu++pVKkUo0Y1Z9y4NtSvX7yuavF4TvzYfO3gwXOUL1+CgQOb4HLF0LVrXTtBbkwRkt8A+E5VbxORDaraxrdso6q2LoBar6moB8Blqso33xxh4sT1zJu3k6ysbPr29Tah69KlbsBe6XLkyHlmztyO2+1hw4ZjhIYKPXtG4nLF0K9fI8qWLeF0icaYq7hWAOR1WuUFEakKqG9jHYAzfqyvWBEROnasTceOtfn7388xadJGpkzZzKJFu2nZsjoTJsQyfHhTypQp+h+Y58+n8/77u0lM9LByZTLZ2cptt91KQkI3hg6NokYNa75mTKDK6xFALPAa0BzYClQHBqnq5oIt778FyhHA1Vy6lMGMGdtJSFjHli0nqFq1DGPGtOTJJ1sTHl7B6fL+S2ZmNitXJuN2e1i4cBcXL2YSGVnxx+ZrdrWTMYEl31cB+cb9owABdqhqhn9LzF0gB8Blqsrnnx8kIWE9ixbtJjQ0hIEDGxMf35YOHWo5NjykqqxffxS3O4mZM5M4evQilSuXYsiQprhc0XTqVCcor2wypjjI7zmAh6+2XFXf80NteVYcAuBK+/ZdbkK3lTNn0rjttluJj4/loYeiKFmycCZH7d9/hhkzvJO0kpJOUaJECA884J2k1atXfUqVsuZrxgS6/AbAa1c8LA10B9ar6iD/lZi74hYAl50/n857721j4sQN7Nhxilq1yvHEE60ZO7ZlgYyx//BDKvPm7SQx0cPq1SkAdO4cjssVzaBBTahSxZqvGVOc+HUimIhUBmap6n35LOo+IAEIBd5U1b9c7/XFNQAuy85WPvpoPwkJ6/jww/2ULBnKsGFNiY+PpU2bmvnadnp6FsuX78Xt9rBkyV7S07OIiqpCXFwMw4c3LXaXqRpj/sPfAVAC2KqqUfkoKBTYCfQAUvB2GR2mqp5rrVPcA+BK27ef5LXXNvDuu9u4cCGDzp3DiY+PpV+/Rnm+wbmq8vXXh0lM9DBnzg5OnUqlRo2yDBvmnaTVtm3NgL0k1RiTd/kdAlqC7xJQIASIAeao6q/yUVBH4Peq2tP3+FkAVX3xWusEUwBcdvp0Km+/vZXXXlvP/v1nqVmzLMOGRRMXF0ObNjWu+gG+c+cp3G7vJK29e89QpkwYDz7YGJcrmh49IvMcIMaY4iG/AXD3FQ8zgWRVTclnQYOA+1T1Md/jOKC9qo7L8boxwBiAiIiItsnJyfl524CVlZXN0qV7effdbSxduoeMjGyio6sQF9fMN6cgjNmzd+B2e/j22+8JCRG6d49gxAhv87UKFYp3SwpjzLUVuVtC5jUArhSMRwBXc+rUpR9P4n755SEAQkKE7GylVavqvuZr0XaTe2MMcJMzgUXkHP8Z+vmvpwBV1fy0dzwE1L3icbhvmcnF2bPpZGZmU6lSKUqVCiUtLYvsbKVu3QqMHNmCRx9tRsWK1m7ZGHN9Th4BhOE9Cdwd7wf/d8BwVd12rXWC+Qhgy5bjTJ+exOLFu0lKOgVAgwaV6N27Ab16NeDMmTQSEtbz9deHqVChJCNHNmf8+Fi7taIxxj9DQCJSA+88AABU9UA+i+oFvIr3MtC3VfVP13t9sAXAoUPnfJO0kti8+ThhYSF06VLX96Ffn8aNb/mfk8DffXeEhIT1zJ69g6ysbHr3bkB8fFu6d4+wK36MCVL5PQncF3gZqA0cA+oBSarazN+FXk8wBMC5c+nMn78Tt9vDqlUHUIUOHWrhcsUweHAU1auXzdN2Dh8+z+TJG5k8eRPHj1+iWbOqTJgQi8sVY107jQky+Q2ATUA3YKWqthGRroBLVUf5v9RrK64BkJGRxccfJ5OY6GHRot1cupRJw4aVcbmiGTEihsaNb7npbaemZjJr1nYSEtazceMxqlQpzejRLXnqqdZ2hy5jgkR+A2CtqrbzBUEbVc0OtpvC+5uqsnbt9z82Xzt+/BJVqpRm6FDvJC1/N4ZTVb788hAJCetYuHA3IvDgg42Jj4/ljjvq2PCQMcVYfu8HcFpEygNfANNF5BhwwZ8FBot9+04zfbp3XH/HjlOUKhVKnz7e5mv33Ve/wJrAiQidO4fTuXM4ycln+Oc/NzJt2mbmzdtJbGxN4uNjGTIkypq/GRNE8noEUA64hHcW8AigEjBdVU8WbHn/LVCPAH74IZU5c7yTtC5ft3/33eG4XDEMGtSEypVL57KFgnHhQjqJiR4mTlxPUtIpatQoy+OPt+KJJ1pz6612oxdjiov8DgE9A8xWVUev0w+kAEhLy2TZsr243UksW+ZtvhYTU9XXfC2aiIiiM/6uqqxcmUxCwnqWLdtLiRIhDBnibULXrt2tTpdnjMmn/AbA74DBwClgNjBXVY/6vcpcFPUAyM5WvvrqEImJHubO3cHp02nUrFmW4cO9vXtat756756iZNeuH3jttfX8619bOX8+g06dahMfH8uAAU2sh5AxAcpf8wBaAkOAgUCKqt7jvxJzV1QDICnp5I/N15KTz1KuXAkGDGiMyxVDt24RAfnBeeZMGu+8s5XXXtvAnj2nCQ+vwFNPtWb06JZUrWr3CzAmkPgrAG4FHgKGAhVUtaX/SsxdUQqAo0cvMHPmdtxuD+vWHSUkROjRox5xcTH079+IcuWKR/O1rKxsli/fR0LCOj755ABlyoThcsUwYUIbmjev7nR5xpg8yO8Q0JN4h4CqA3PxtoK+Zt/+guJ0AFy4kM6iRXtITNzGxx8nk5WltG1bE5crhqFDmxb7E6dbtx5n4sQNJCZ6SE3NpFu3COLjY+nduwGhoYF3lGNMsMhvALyI9yTwxoIoLq+cCICsrGxWrTpAYqKHBQt2ceFCBhERFXC5YhgxIpqYmGqFWk9RcPLkJd58czOvv76RlJRzNGhQiXHj2jByZAsqVbImdMYUNfkeAvLdwasmV8wdyG8voBtVWAGgqmzadBy328OMGUkcOXKBSpVK8dBDTYiLi+HOO8MJCSnaJ3MLQ2ZmNgsX7iIhYT1ffXWI8uVL8OijzRk/vg1NmlRxujxjjE9+jwDGAb8HjgLZvsVa3M4BHDx4lhkzkkhM9LBt20lKlAihV68GxMXF0Lt3A0qXtklS17Ju3fckJKxn1qztZGRk06tXfeLj29KjR70if+WTMcVdfgNgN96btRTqxK+cCiIAzpxJY948b/O1zz8/iCp06lT7x+ZrdsXLjfn++wtMmbKJN97YyNGjF4mOrsL48bE8/HBMsTkxbkygyW8AfAr0UNXMgigur/wVAOnpWXz44T7cbg+LF+8hLS2Lxo1v+bH5mvXQz7+0tEzmzNlBQsJ61q07SuXKpXjssRY89VQbIiMrOV2eMUElvwHwFhAFLAPSLi9X1Vf8WWRu8hMAqso33xzB7fYwe/YOTp68RLVqZRgyJIq4uGbcfvutNlRRAFSVf//7MAkJ61iwYBeq0K9fI+LjY7nrrnD7MzemEOS3GdwB309J309AmTp1E3/963fs2XOa0qXD6NevIS5XDD17RlKiRME0XzNeIsIdd9ThjjvqcPDgWSZN2sjUqZtZuHAXrVpVJz4+lmHDou38ijEOuNGJYGVV9WIB1nNdN3sE8NJLa/joo2RcrmgGDmxi98t12MWLGUyfnkRCwjq2bTtJ9eplGDvW24TObmRvjP/ldwioI/AWUF5VI0SkFTBWVZ/0f6nXdrMBoKo21FAEqSqrVh0gIWE9S5fuITQ0hMGDo5gwIZb27Ws5XZ4xxca1AiCv0zdfBXoCJwFUdRNwl//KK1j24V80iQjdu9dj8eIH2bXrMcaNa8PSpXvo0GE6HTpMZ+bMJDIyspwu05hiK8/z91X1YI5F9n+m8ZuGDSvzj390JSXlcSZO7MbJk5cYPnwZkZHT+NOfvuH4ccdGHo0ptvIaAAdFpBOgIlJCRH4GJBVgXSZIVahQkvHjY9mxYxRLlz5Is2ZVef75L6lbdwojR37Ipk3HnC7RmGIjrwHwOPAUUAc4BLT2PTamQISECL17N+Sjjx5i27ZH+clPmjN79nZat36PLl1msXDhLrKysnPfkDHmmm7oKiCnOd0N1Djrhx9SfU3oNnDgwDkiIysyblwbRo1q4dhtNY0JBPm9CmjiVRafAdaq6iI/1JcnFgAGvE3oFi/eTULCer74IoWyZcN45JFmTJgQS9OmVZ0uz5giJ79XAZXGO+yzy/fTEggHRonIq36r0pg8CAsLYcCAJnz++VA2bHiYIUOa8hGFU4wAABZSSURBVPbbW4mO/hf33TePDz7YS3Z24BzZGuOUvB4BfAPcoapZvsdhwGrgTmCLqsYUaJU+dgRgruXYsQtMnbqZSZM2cuTIBZo0uYXx49vw6KPNKV8+4CavG+NX+T0CuAW4copmOaCKLxDSrr6KMYWnRo1yPP98R/bvH8OMGb2pXLkU48evok6dyTzzzKfs3Xva6RKNKXLyGgB/BTaKyL9E5B1gA/A3ESkHrCyo4oy5USVLhjJsWDRr1rj45psR9O7dgNde20CjRm/Sv//7fPrpAQLpwgdjCtKN3BGsFnC77+F3qnq4wKq6BhsCMjfj0KFzvPHGJqZM2cSJE5do2bI6EybEMnx4U8qUKeF0ecYUuJsaAhKRpr7/xgK1gIO+n1t9y4wp8urUqcALL9zJgQNjeOutngA89tgK6tadynPPrSYl5ZzDFRrjjOseAYjINFUd7bshTE6qqt0KrrT/ZUcAxh9UlS++SCEhYT2LFu1GBAYNasKECbF07FjbekeZYiffN4UvCiwAjL/t33+G11/fwJtvbuHMmTTatatJfHxbBg+OomRJu1eEKR5udgjoF1f8/lCO5/7sv/KMcUZkZCX+/vcupKSM5Z//7M65cxnExS2nXr2p/OEP/+bo0QtOl2hMgcntKqChV/z+bI7n7vNzLcY4pnz5kjz5ZBs8np/w4YcDadOmBr/73b+JiJjKo49+wIYNR50u0Ri/yy0A5Bq/X+2xMQEvJETo2bM+y5cPZPv2kYwe3YJ583YSG5tI584zmTdvB5mZ1oTOFA+5BYBe4/erPc4zEfmbiGwXkc0islBEKt/stowpKFFRVXj99XtISRnLyy934dCh8zz00BIaNpzGX//6LadOXXK6RGPyJbergLKAC3i/7ZcBLt+VQ4DSqnpTF1GLyL3AKlXNFJGXAFT1l7mtZyeBjZOysrJZunQvCQnr+PTTg5QpE0ZcXAwTJsTSrFk1p8sz5ppu6iSwqoaqakVVraCqYb7fLz++6Rk0qvqRqmb6Hn6Dt7GcMUVaaGgI/fo1YtWqIWza9AjDh0fz3nsemjd/hx495rJ06R5rQmcCiuOXgYrIEmC2qrqv8fwYYAxARERE2+Tk5MIsz5jrOnHi4o9N6A4dOk+jRpUZPz6WRx9tRsWKpZwuzxjAgXkAIrISuPUqT/368j0EROTXQDtggOahEBsCMkVVRkYWCxbsIiFhPV9/fZgKFUoycmRzxo+PpWFDO8VlnFXkJoKJyKPAWKC7qubpjt8WACYQfPfdERIS1jNnjveKod69GxAf35bu3SNslrFxRJEKABG5D3gFuFtVj+d1PQsAE0iOHDnP5MmbmDx5E8eOXaRZs6pMmBCLyxVD2bLWhM4UnqIWALuBUsBJ36JvVPXx3NazADCBKC0tk9mzd/Dqq+vYsOEYt9xSmtGjW/DUU22IiKjodHkmCBSpALhZFgAmkKkqX311iISE9SxYsAsRePDBxsTHx3LHHXVseMgUmPzeEcwYk08iwp13hjN3bl/27RvNT3/ajk8+OUDnzrNo187Ne+9tIy0tM/cNGeMnFgDGOCAioiIvvXQ3Bw+OYfLkHqSmZvLIIx9Qr95Ufv/7r/j+e2tCZwqeDQEZUwSoKitXJjNx4nqWLdtLWFgIQ4Y0JT4+lnbtrnY1tTF5Z0NAxhRhIkKPHpEsWTKAnTtH8cQTrVm0aDe33ebmjjtmMGfOdjIyspwu0xQzFgDGFDGNGt1CQkI3UlLG8uqrXTl69CJDhiylQYM3efHFNZw8aU3ojH9YABhTRFWsWIr4+Lbs2DGSxYsfJCrqFp57bjXh4VMYPXoFW7bkeQqNMVdlAWBMERcaGkKfPg1ZuXIwW7c+ysMPxzB9ehItW75L9+5zWLRoN1lZdo8Cc+MsAIwJIM2aVWPKlHs5eHAsf/lLZ3bt+oH+/d+nSZO3+Mc/1nLmTJrTJZoAYgFgTACqWrUMv/xle/buHc3cuX2oXbs8zzzzGeHhkxk3biU7d55yukQTACwAjAlgYWEhDBoUxerVw1i71sWAAY2ZNm0LUVFv06vXfFas2EcgXeptCpfNAzCmmDl69AKTJ2/ijTc2cvToRZo2rcKECbE8/HAM5cqVdLo84wCbB2BMkKhZsxy/+10nkpPH8N5791OuXAmefHIl4eFT+PnPPyM5+YzTJZoiwo4AjCnmVJV///swEyeuZ/78nahC//6NmDAhlrvuCrcmdEHAuoEaYzh48CyTJm1k6tTNnDqVSqtW1YmPj2XYsGhKlw5zujxTQGwIyBhD3boVefHFu0hJGcu0afeSlaWMHLmCiIgp/OY3X3L48HmnSzSFyI4AjAliqsqnnx4kIWEdS5bsITQ0hIceakJ8fFvat6/ldHnGT+wIwBjzP0SEbt0iWLToQXbteozx49uwbNleOnSYTocO05k5M8ma0BVjFgDGGAAaNqzMK690JSXlcV57rRunTqUyfPgyIiOn8ac/fcPx4xedLtH4mQWAMea/VKhQknHjYtm+fSTLlg2gefNqPP/8l9StO4WRIz9k06ZjTpdo/MQCwBhzVSEhQq9eDVixYhAez08YObI5s2dvp3Xr9+jSZRYLFuy0JnQBzgLAGJOr6OiqTJrUg5SUx/nb3+5m//6zDBy4mIYN3+Tvf/+OH35IdbpEcxMsAIwxeXbLLaX52c9uY/fux5g/vy/16lXk5z//nPDwyTz55MckJZ10ukRzAywAjDE3LCwshAEDmvD550PZsOFhhgxpyttvbyUm5l/07DmP5cv3kp0dOJeYByubB2CM8Ytjxy4wdepmJk3ayJEjF2jS5BbGj2/DI480p0IFa0LnJJsHYIwpUDVqlOP55zuyf/8Ypk/vTeXKpRg/fhXh4ZP5v//7lL17TztdosnBAsAY41clS4YyfHg0a9a4+OabEfTu3YDXX99Ao0Zv0q/fQlatOmD3KCgiLACMMQWmfftazJjxAPv3j+a55zrw738fpnv3ObRs+S7Tpm3m0qUMp0sMahYAxpgCV6dOBV544U4OHhzLW2/1JCREGDPmI8LDp/Dss1+QknLO6RKDkp0ENsYUOlXl888PMnHiBhYt2o0IDBzYhPj4WDp2rG33KPAzOwlsjCkyRIQuXSJYsKAfe/Y8xtNPt2XFiv3cccdMbr/djdvtIT3dmtAVNAsAY4yjIiMr8fe/dyElZSyTJt3D+fMZxMUtp169qfy///dvjh694HSJxZYNARljipTsbOXjj/eTkLCeDz7YR8mSoQwdGkV8fFtiY2s6XV5AsiEgY0xACAkRevasz/LlA9m+fSSjR7dg/vxdtG2bSOfOM5k7dweZmdaEzh8sAIwxRVZUVBVef/0eUlLG8vLLXTh06DyDBy+hQYNpvPTSGk6duuR0iQHNAsAYU+RVrlyaZ55px65do3j//f40alSZX/1qNeHhUxg79iO2bTvhdIkBydEAEJGfioiKSDUn6zDGBIbQ0BD69WvEqlVD2Lz5EUaMiOa99zw0b/4OPXrMZenSPdaE7gY4FgAiUhe4FzjgVA3GmMDVokV1pk3rycGDY/jznzuTlHSSPn0W0qTJWyQkrOPs2TSnSyzynDwC+AfwC8Di2hhz06pVK8uzz7Zn377RzJr1ADVrluXppz8lPHwK8fGr2L37B6dLLLIcCQAR6QccUtVNeXjtGBFZKyJrjx8/XgjVGWMCUYkSoQwZ0pSvvhrOt9+OoG/fhrzxxkaaNHmLPn0W8PHH+60JXQ4FNg9ARFYCt17lqV8DzwH3quoZEdkPtFPVXM/i2DwAY8yNOHLkPJMnb2Ly5E0cO3aRmJiqTJgQS1xcDGXLlnC6vEJzrXkAhT4RTERaAJ8AF32LwoHDwO2q+v311rUAMMbcjLS0TGbN2k5Cwno2bDjGLbeUZvToFjz1VBsiIio6XV6BKzIB8D8F2BGAMaaQqCpffXWIhIT1LFiwCxF48MHGTJgQy5131im2TehsJrAxJuiJCHfeGc7cuX3Zt280P/1pOz755AB33TWLtm0TeffdraSlZTpdZqFx/AjgRtgRgDHG3y5cSMftTmLixPV4PCepXr0Mjz/eiieeaE2tWuWdLs8viuwQ0I2wADDGFBRV5ZNPDpCQsI5ly/YSFhbC4MFRxMfHcttttZwuL19sCMgYY65DRLjnnnosWTKAnTtH8cQTrVm8eA+33z6dTp1mMHv2djIyitc9CiwAjDEmh0aNbiEhoRspKWN59dWuHDt2kaFDl1K//jRefHENJ05czH0jAcACwBhjrqFixVLEx7dl585RLFnyINHRVXnuudXUrTuVxx5bwZYtgT051QLAGGNyERIiPPBAQz7++CG2bn2Uhx+OYcaMJFq2fJdu3Wbz/vu7yMoKvHsUWAAYY8wNaNasGlOm3EtKylj+8pfO7N59mgcfXETjxm/xyitrOXMmcJrQWQAYY8xNqFKlDL/8ZXv27h3N3Ll9qFOnPD/96WfUqTOZceNWsnPnKadLzJUFgDHG5ENYWAiDBkWxevUw1q51MXBgY6ZN20JU1Nv06jWfFSv2FdkmdDYPwBhj/Ozo0QtMmbKJN97YxPffX6Bp0yo/NqErX75koddj8wCMMaaQ1KxZjt/+thPJyWNITOxF+fIlePLJlYSHT+FnP/uM/fvPOF0iYEcAxhhT4FSVr78+TELCeubP34kq9O3bkKefbstdd4UXeBM6OwIwxhiHiAidOtVh9uw+7N8/hl/+8nZWrz5Ely6zadPmPd5+ewupqYXfhM4CwBhjClF4eAX+/OfOHDw4hmnT7iUrSxk1agV1607h+ee/5PDh84VWiw0BGWOMg1SVzz47SELCehYv3k1oaAgPPdSE+Pi2tG/vnyZ0NgRkjDFFkIjQtWsE77/fn927H2P8+DYsW7aXDh2m06HDdGbOTCI9vWCa0FkAGGNMEdGgQWVeeaUrKSmP89pr3Th1KpXhw5cRGTmVVasO+P39LACMMaaIqVChJOPGxbJ9+0iWLx9A69Y1aNy4st/fJ8zvWzTGGOMXISHC/fc34P77GxTM9gtkq8YYY4o8CwBjjAlSFgDGGBOkLACMMSZIWQAYY0yQsgAwxpggZQFgjDFBygLAGGOCVEA1gxOR40DyTa5eDTjhx3ICge1zcLB9Dg752ed6qlo958KACoD8EJG1V+uGV5zZPgcH2+fgUBD7bENAxhgTpCwAjDEmSAVTAEx1ugAH2D4HB9vn4OD3fQ6acwDGGGP+WzAdARhjjLmCBYAxxgSpYhcAInKfiOwQkd0i8qurPF9KRGb7nl8jIpGFX6V/5WGfnxERj4hsFpFPRKSeE3X6U277fMXrBoqIikhAXzKYl/0VkcG+v+dtIjKjsGv0tzz8u44QkU9FZIPv33YvJ+r0JxF5W0SOicjWazwvIjLR92eyWURi8/WGqlpsfoBQYA/QACgJbAJicrzmSWCy7/ehwGyn6y6Efe4KlPX9/kQw7LPvdRWAL4BvgHZO113Af8eNgQ3ALb7HNZyuuxD2eSrwhO/3GGC/03X7Yb/vAmKBrdd4vhfwASBAB2BNft6vuB0B3A7sVtW9qpoOzAL65XhNP+Bd3+/zgO4iIoVYo7/lus+q+qmqXvQ9/AYIL+Qa/S0vf88AfwReAlILs7gCkJf9HQ38U1V/AFDVY4Vco7/lZZ8VqOj7vRJwuBDrKxCq+gVw6jov6Qe8p17fAJVFpNbNvl9xC4A6wMErHqf4ll31NaqaCZwBqhZKdQUjL/t8pVF4v0EEslz32XdoXFdVlxVmYQUkL3/HTYAmIvKViHwjIvcVWnUFIy/7/HvAJSIpwHJgfOGU5qgb/f/9uuym8EFERFxAO+Bup2spSCISArwCPOpwKYUpDO8wUBe8R3hfiEgLVT3taFUFaxjwjqq+LCIdgUQRaa6q2U4XFiiK2xHAIaDuFY/Dfcuu+hoRCcN76HiyUKorGHnZZ0TkHuDXQF9VTSuk2gpKbvtcAWgOfCYi+/GOlS4O4BPBefk7TgEWq2qGqu4DduINhECVl30eBcwBUNWvgdJ4G6YVZ3n6/z2vilsAfAc0FpH6IlIS70nexTlesxh4xPf7IGCV+s6uBKhc91lE2gBT8H74B/rYMOSyz6p6RlWrqWqkqkbiPe/RV1XXOlNuvuXl3/X7eL/9IyLV8A4J7S3MIv0sL/t8AOgOICLReAPgeKFWWfgWAw/7rgbqAJxR1SM3u7FiNQSkqpkiMg5YgfcqgrdVdZuI/AFYq6qLgbfwHiruxnuyZahzFedfHvf5b0B5YK7vfPcBVe3rWNH5lMd9LjbyuL8rgHtFxANkAT9X1YA9ss3jPv8UmCYi/4f3hPCjAf5lDhGZiTfIq/nObfwOKAGgqpPxnuvoBewGLgI/ydf7BfiflzHGmJtU3IaAjDHG5JEFgDHGBCkLAGOMCVIWAMYYE6QsAIwxJkhZABhjTJCyADDGmCBlAWCKNF+/9545lj0tIm9cZ53zhVDXBBFJEpHpN7Huv2/w9WNz7q+IbPXNfkVEyojI5yISeo31S4rIF77WJ8b8yALAFHUz+d/Z2kN9y530JNBDVUfc6Iqq2ukGV2kBrL/8QERKA5F4+/0AjAQWqGrWNd4vHfgEGHKjtZrizQLAFHXzgN6+fjD47uBWG1gtIu+LyDrfHbDG5FxRRCKvvLOSiPxMRH5/xWOXiHwrIhtFZMrVvkGL925qW30/T/uWTcZ7o5IPfG0Icq7ziK+uzSLy5VWeP39FfUkiMs23Dx+JSJmr/Bm05IoAwBsIO6/4wB8BLMrl/d/3vc6YH9khoSnSVPWUiHwL3I/3Q24oMEdVVURG+p4vA3wnIvPz2v/GN3wyBLhDVTNEZBLeD8j3rnhNW7y9VtrjvQPTGhH5XFUf9/Xb76qqJ3JstwLwS6C1qqaLSOVcSmkMDFPV0SIyBxgIuHO8phmwQEQu920pDyz1vV9JoIGq7s/l/bcCt+Xlz8YEDzsCMIHgymGgK4d/JojIJrzdPutyY+2PuwNt8QbHRt/jBjlecyewUFUvqOp5YAHQOZftZgFlgJdFpF0e+vHvU9WNvt/X4R3a+ZGI1AWOq2q9K7qbzgK2+F5SDbjyPa76/r6jhXRfQBgDWACYwLAI7607Y/He23idiHQB7gE6qmorvPfDLZ1jvUz++9/4lc8L8K6qtvb9RKnq7/NbqO/Wm82Br4CpIvJkLqtceW+GLP73qLwFsC3Hshhgs+/3S1yxX7m8fykC//aYxo8sAEyR5/v2/SnwNv/59l8J+EFVL4pIU7w3fcnpKFBDRKqKSCnggSue+wQYJCI1AESkiojUy7H+aqC/iJQVkXLAg75l1yQijX1HDLPwDtPkDKUb1RLw5FjWDN8RgO8ewKG+E8PXfH8RqQqcUNWMfNZjihELABMoZgKt+E8AfAiEiUgS8Be8w0D/xfdh9wfgW+BjYPsVz3mA54GPRGSz7/laOdZfD7zjW38N8Kaqbsilzl+LyA4RWQ/UBybd2G7+jxZcEQAiUgVvG/fvr3jNR3iHq673/l2B4nB/ZONHdj8AYwKcb2js/1Q17jqvWQD8SlV3Xus1JvjYEYAxAc53pPLp9SaCAe/bh7/JyY4AjDEmSNkRgDHGBCkLAGOMCVIWAMYYE6QsAIwxJkhZABhjTJCyADDGmCD1/wE9py980jRWEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, y_ground, color='darkblue')\n",
    "plt.plot(x, y_excited, color='darkblue')\n",
    "plt.title(\"$Z$ to [[-5, 0], [0, 4]]\")\n",
    "plt.xlabel('Value of $s$ in $H(s)$')\n",
    "plt.ylabel('Eigenvalues')\n",
    "plt.show()\n",
    "# plt.savefig('adiabatic_crossing.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adiabatic evolution for 2x1 Hubbard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.utils import HubbardSquareLattice \n",
    "\n",
    "x_n = 2 \n",
    "y_n = 1\n",
    "n_dofs = 1\n",
    "periodic = 0\n",
    "spinless = 0 \n",
    "\n",
    "lattice = HubbardSquareLattice(x_n, y_n, n_dofs=n_dofs, periodic=periodic, spinless=spinless)\n",
    "\n",
    "from openfermion.hamiltonians import FermiHubbardModel\n",
    "from openfermion.utils import SpinPairs\n",
    "tunneling = [('neighbor', (0, 0), 1.)] \n",
    "interaction = [('onsite', (0, 0), 2., SpinPairs.DIFF)] \n",
    "potential = None\n",
    "mag_field = 0. \n",
    "particle_hole_sym = False \n",
    "\n",
    "hubbard = FermiHubbardModel(lattice , tunneling_parameters=tunneling, interaction_parameters=interaction, \n",
    "                            potential_parameters=potential, magnetic_field=mag_field, \n",
    "                            particle_hole_symmetry=particle_hole_sym)\n",
    "# 4 qubits so 2^4 = 16 dimensional state space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0 [0^ 2] +\n",
       "-1.0 [1^ 3] +\n",
       "-1.0 [2^ 0] +\n",
       "-1.0 [3^ 1]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hubbard.tunneling_terms()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.transforms import get_sparse_operator \n",
    "from openfermion.utils import get_ground_state\n",
    "\n",
    "H_A = get_sparse_operator(hubbard.tunneling_terms())\n",
    "H_B = get_sparse_operator(hubbard.hamiltonian())\n",
    "H = lambda s: (1-s)*H_A + s*H_B\n",
    "psi_A = get_ground_state(H_A)[1] # We can use this because no degeneracy in 2x1 model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.0000000e+00+0.0000000e+00j,  0.0000000e+00+0.0000000e+00j,\n",
       "        0.0000000e+00+0.0000000e+00j, -4.9999997e-01+1.2246467e-16j,\n",
       "        0.0000000e+00+0.0000000e+00j, -3.0616168e-17+5.6240991e-33j,\n",
       "        4.9999997e-01-6.1232336e-17j,  0.0000000e+00+0.0000000e+00j,\n",
       "        0.0000000e+00+0.0000000e+00j, -4.9999997e-01+6.1232336e-17j,\n",
       "        9.1848507e-17-5.6240995e-33j,  0.0000000e+00+0.0000000e+00j,\n",
       "       -4.9999997e-01+0.0000000e+00j,  0.0000000e+00+0.0000000e+00j,\n",
       "        0.0000000e+00+0.0000000e+00j,  0.0000000e+00+0.0000000e+00j],\n",
       "      dtype=complex64)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Using final_wavefunction and prepare_gaussian_state \n",
    "from openfermion import get_quadratic_hamiltonian\n",
    "from openfermioncirq import prepare_gaussian_state\n",
    "from cirq import Circuit, final_wavefunction, LineQubit\n",
    "psi_A = final_wavefunction(Circuit(prepare_gaussian_state(\n",
    "    LineQubit.range(4), \n",
    "    get_quadratic_hamiltonian(hubbard.tunneling_terms()))))\n",
    "psi_A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.0000000e+00+0.0000000e+00j  0.0000000e+00+0.0000000e+00j\n",
      "  0.0000000e+00+0.0000000e+00j -4.9999997e-01+1.2246467e-16j\n",
      "  0.0000000e+00+0.0000000e+00j -3.0616168e-17+5.6240991e-33j\n",
      "  4.9999997e-01-6.1232336e-17j  0.0000000e+00+0.0000000e+00j\n",
      "  0.0000000e+00+0.0000000e+00j -4.9999997e-01+6.1232336e-17j\n",
      "  9.1848507e-17-5.6240995e-33j  0.0000000e+00+0.0000000e+00j\n",
      " -4.9999997e-01+0.0000000e+00j  0.0000000e+00+0.0000000e+00j\n",
      "  0.0000000e+00+0.0000000e+00j  0.0000000e+00+0.0000000e+00j]\n"
     ]
    }
   ],
   "source": [
    "# Let's do 5 steps: H(0), H(0.25), H(0.5), H(0.75), H(1) at t=1 \n",
    "n = 5\n",
    "t = 1\n",
    "res = psi_A\n",
    "print(res)\n",
    "for i in range(n): \n",
    "    s = i / (n-1)\n",
    "#     res = np.dot(sp.expm(-1j * H(s) * t), res)\n",
    "    e = scipy.linalg.expm(-1j * H(s) * t)\n",
    "    res = e.dot(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j,\n",
       "        0.00000000e+00+0.00000000e+00j,  6.59063564e-02-3.53071442e-01j,\n",
       "        0.00000000e+00+0.00000000e+00j, -3.06161680e-17+5.62409909e-33j,\n",
       "       -3.58082448e-02+6.08041626e-01j,  0.00000000e+00+0.00000000e+00j,\n",
       "        0.00000000e+00+0.00000000e+00j,  3.58082448e-02-6.08041626e-01j,\n",
       "        9.18485073e-17-5.62409946e-33j,  0.00000000e+00+0.00000000e+00j,\n",
       "        6.59063564e-02-3.53071442e-01j,  0.00000000e+00+0.00000000e+00j,\n",
       "        0.00000000e+00+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.utils import inner_product \n",
    "\n",
    "def overlap(a, b):\n",
    "    \"\"\"Calculates the overlap between vectors a and b. This metric is also known as fidelity. \"\"\"\n",
    "    inner = inner_product(a, b)\n",
    "    return (np.conjugate(inner) * inner).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9964795950032181"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overlap(res, get_ground_state(H_B)[1]) # With 5 steps, we get 0.99648 overlap. Increases if we increase steps. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing out adiabatic evolution for 2x2 Hubbard. This tunneling term has degnerate ground states so we have to be more specific which one we want to start with. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_n = 2\n",
    "lattice = HubbardSquareLattice(x_n, y_n, n_dofs=n_dofs, periodic=periodic, spinless=spinless)\n",
    "hubbard = FermiHubbardModel(lattice , tunneling_parameters=tunneling, interaction_parameters=interaction, \n",
    "                            potential_parameters=potential, magnetic_field=mag_field, \n",
    "                            particle_hole_symmetry=particle_hole_sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Perturbed eigenvector v_per[:, 0] has overlap 0.9770315592604035.\n",
      "Perturbed eigenvector v_per[:, 1] has overlap 1.016132559310905e-20.\n",
      "Perturbed eigenvector v_per[:, 2] has overlap 8.801403183749669e-26.\n",
      "Perturbed eigenvector v_per[:, 3] has overlap 3.619538211946138e-26.\n",
      "Perturbed eigenvector v_per[:, 4] has overlap 2.0798430394317092e-26.\n",
      "7.365761758597068e-12 0\n",
      "0.029993297631743365 1\n",
      "0.14076420572868092 4\n",
      "0.1785725708397301 34\n",
      "Tunneling eigenvector v_tun[:, 34] had the maximum overlap of 0.1785725708397301 with the best perturbed state.\n"
     ]
    }
   ],
   "source": [
    "from openfermion.transforms import get_sparse_operator\n",
    "\n",
    "# Lowest eigenvalue/vector of Hubbard Hamiltonian\n",
    "# The column v_hub[:, i] is the eigenvector corresponding to the eigenvalue w_hub[i]\n",
    "hub_sparse = get_sparse_operator(hubbard.hamiltonian())\n",
    "w_hub, v_hub = scipy.sparse.linalg.eigsh(hub_sparse, k=1, which='SA') \n",
    "\n",
    "# From looking at the spectrum,  I noticed the ground state was degenerate. \n",
    "# There were 16 eigenvectors with the lowest eigenvalue. \n",
    "# Actually, the initial state might not even be a ground state. Let's try \n",
    "# as many of the eigenvectors as we can. \n",
    "tun_sparse = get_sparse_operator(hubbard.tunneling_terms())\n",
    "# k can be at most n-2 where n is the dimension of the matrix\n",
    "w_tun, v_tun = scipy.sparse.linalg.eigsh(tun_sparse, k=np.shape(tun_sparse)[0]-2, which='SA')\n",
    "\n",
    "# Do a perturbation and find state with most overlap \n",
    "s = 1e-4\n",
    "int_sparse = get_sparse_operator(hubbard.interaction_terms())\n",
    "perturbed_sparse = tun_sparse + s * int_sparse\n",
    "\n",
    "w_per,v_per = scipy.sparse.linalg.eigsh(perturbed_sparse, k=5, which='SA') \n",
    "\n",
    "for i in range(len(w_per)):\n",
    "    fid = overlap(v_hub[:, 0], v_per[:, i])\n",
    "    print('Perturbed eigenvector v_per[:, {}] has overlap {}.'.format(i, fid))\n",
    "    \n",
    "# WOW, there's a clear winner here v_per[:, 0]\n",
    "per_state_most_overlap = v_per[:, 0]\n",
    "\n",
    "max_overlap = 0\n",
    "index_max_overlap = 0\n",
    "\n",
    "for i in range(len(w_tun)):\n",
    "    fid = overlap(per_state_most_overlap, v_tun[:, i])\n",
    "    if fid > max_overlap: \n",
    "        print(fid, i)\n",
    "        max_overlap = fid \n",
    "        index_max_overlap = i \n",
    "print(\"Tunneling eigenvector v_tun[:, {}] had the maximum overlap of {} with the best perturbed state.\".format(\n",
    "    index_max_overlap, max_overlap))\n",
    "initial_state = v_tun[:, index_max_overlap]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_A = get_sparse_operator(hubbard.tunneling_terms())\n",
    "H_B = get_sparse_operator(hubbard.hamiltonian())\n",
    "H = lambda s: (1-s)*H_A + s*H_B\n",
    "psi_A = initial_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "254"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_tun.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.13819128511606324 4\n",
      "0.1299495385671853 29\n",
      "0.1753085802122591 34\n"
     ]
    }
   ],
   "source": [
    "# Let's do 5 steps: H(0), H(0.25), H(0.5), H(0.75), H(1) at t=1 \n",
    "n = 5\n",
    "t = 3\n",
    "for state in range(40):\n",
    "#     res = psi_A\n",
    "    res = v_tun[:, state]\n",
    "    for i in range(n): \n",
    "        s = i / (n-1)\n",
    "    #     res = np.dot(sp.expm(-1j * H(s) * t), res)\n",
    "        e = scipy.linalg.expm(-1j * H(s) * t)\n",
    "        res = e.dot(res)\n",
    "    x = overlap(res, v_hub[:, 0])\n",
    "    if x > 0.1: \n",
    "        print(x, state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I noticed that every time we run the cell that determines the tunneling eigenvectors, we get different overlaps. I should have noticed this earlier! Since we have a degenerate ground state for the tunneling term, the eigenvectors span some subspace, so we can pick any basis of that subspace as the eigenvectors. This means that letting the algorithm pick gives us a random one. \n",
    "\n",
    "Once we do that, we should try to find a *linear combination* in that ground state subspace that's closest to the true ground state eigenvector of the Hubbard Hamiltonian. We can use Least-Squares to do this!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16 should be 16 for 2x2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/warren/github/qc-mentorship/env/lib/python3.6/site-packages/ipykernel_launcher.py:11: FutureWarning: arrays to stack must be passed as a \"sequence\" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(256, 16)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min_eig = -4 \n",
    "indices_ground = []\n",
    "for i in range(len(w_tun)):\n",
    "    if w_tun[i] - min_eig < tol: \n",
    "        indices_ground.append(i)\n",
    "print(len(indices_ground), 'should be 16 for 2x2')\n",
    "ground_state_array = []\n",
    "for i in indices_ground:\n",
    "    ground_state_array.append(v_tun[:, i])\n",
    "    \n",
    "mat = np.hstack(row.reshape((-1, 1)) for row in ground_state_array)\n",
    "mat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(256,)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "projection = mat @ (np.linalg.inv(mat.T @ mat) @ mat.T @ v_hub[:,0])\n",
    "projection.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9545796915314164"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overlap(projection, v_hub[:, 0]) \n",
    "# THERE WE GO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi_A = projection\n",
    "\n",
    "# Let's do 5 steps: H(0), H(0.25), H(0.5), H(0.75), H(1) at t=1 \n",
    "n = 5\n",
    "t = 5\n",
    "res = psi_A\n",
    "for i in range(n): \n",
    "    s = i / (n-1)\n",
    "#     res = np.dot(sp.expm(-1j * H(s) * t), res)\n",
    "    e = scipy.linalg.expm(-1j * H(s) * t)\n",
    "    res = e.dot(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.974689145161343"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overlap(res, v_hub[:, 0]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9373361730974824"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overlap(res, psi_A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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
}