12_01_quantum_spin_chain.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Spins\n",
    "\n",
    "We are going consider a magnetic insulator, \n",
    "a Mott insulator. In this problem we have a valence electron per\n",
    "atom, with very localized wave-functions. The Coulomb repulsion between\n",
    "electrons on the same orbital is so strong, that electrons are bound to\n",
    "their host atom, and cannot move. For this reason, charge disappears\n",
    "from the equation, and the only remaining degree of freedom is the spin\n",
    "of the electrons. The corresponding local state can therefore be either\n",
    "$|\\uparrow\\rangle$ or $|\\downarrow\\rangle$. The only interaction taking\n",
    "place is a process that “flips” two anti-aligned neighboring spins\n",
    "$|\\uparrow\\rangle|\\downarrow\\rangle \\rightarrow |\\downarrow\\rangle|\\uparrow\\rangle$.\n",
    "\n",
    "Let us now consider a collection of spins residing on site of a\n",
    "–one-dimensional for simplicity– lattice. An arbitrary state of\n",
    "$N$-spins can be described by using the $S^z$ projection\n",
    "($\\uparrow,\\downarrow$) of each spin as: $|s_1,s_2,..., s_N\\rangle$. As\n",
    "we can easily see, there are $2^N$ of such configurations.\n",
    "\n",
    "Let us now consider a collection of spins residing on site of a\n",
    "–one-dimensional for simplicity– lattice. An arbitrary state of\n",
    "$N$-spins can be described by using the $S^z$ projection\n",
    "($\\uparrow,\\downarrow$) of each spin as: $|s_1,s_2,..., s_N\\rangle$. As\n",
    "we can easily see, there are $2^N$ of such configurations.\n",
    "\n",
    "We shall describe the interactions between neighboring spins using the\n",
    "so-called Heisenberg Hamiltonian:\n",
    "\n",
    "$$\\hat{H}=\\sum_{i=1}^{N-1} \\hat{\\mathbf{S}}_i \\cdot \\hat{\\mathbf{S}}_{i+1}\n",
    "$$\n",
    "\n",
    "where $\\hat{\\mathbf{S}}_i = (\\hat{S}^x,\\hat{S}^y,\\hat{S}^z)$ is the spin\n",
    "operator acting on the spin on site $i$.\n",
    "\n",
    "Since we are concerned about spins one-half, $S=1/2$, all these\n",
    "operators have a $2\\times2$ matrix representation, related to the\n",
    "well-known Pauli matrices:\n",
    "\n",
    "$$S^z = \\left(\n",
    "\\begin{array}{cc}\n",
    "1/2 & 0 \\\\\n",
    "0 & -1/2\n",
    "\\end{array}\n",
    "\\right),\n",
    "S^x = \\left(\n",
    "\\begin{array}{cc}\n",
    "0 & 1/2 \\\\\n",
    "1/2 & 0\n",
    "\\end{array}\n",
    "\\right),\n",
    "S^y = \\left(\n",
    "\\begin{array}{cc}\n",
    "0 & -i/2 \\\\\n",
    "i/2 & 0\n",
    "\\end{array}\n",
    "\\right),\n",
    "$$\n",
    "\n",
    "These matrices act on two-dimensional vectors defined by the basis\n",
    "states $|\\uparrow\\rangle$ and $|\\downarrow\\rangle$. It is useful to\n",
    "introduce the identities:\n",
    "$$\\hat{S}^\\pm = \\left(\\hat{S}^x \\pm i \\hat{S}^y\\right),$$ where $S^+$\n",
    "and $S^-$ are the spin raising and lowering operators. It is intuitively\n",
    "easy to see why by looking at how they act on the basis states:\n",
    "$\\hat{S}^+|\\downarrow\\rangle = |\\uparrow\\rangle$ and\n",
    "$\\hat{S}^-|\\uparrow\\rangle = |\\downarrow\\rangle$. Their corresponding\n",
    "$2\\times2$ matrix representations are: $$S^+ = \\left(\n",
    "\\begin{array}{cc}\n",
    "0 & 1 \\\\\n",
    "0 & 0\n",
    "\\end{array}\n",
    "\\right),\n",
    "S^- = \\left(\n",
    "\\begin{array}{cc}\n",
    "0 & 0 \\\\\n",
    "1 & 0\n",
    "\\end{array}\n",
    "\\right),\n",
    "$$\n",
    "\n",
    "We can now re-write the Hamiltonian (\\[heisenberg\\]) as:\n",
    "$$\\hat{H}=\\sum_{i=1}^{N-1}\n",
    "\\hat{S}^z_i \\hat{S}^z_{i+1} +\n",
    "\\frac{1}{2}\\left[\n",
    "\\hat{S}^+_i \\hat{S}^-_{i+1} +\n",
    "\\hat{S}^-_i \\hat{S}^+_{i+1}\n",
    "\\right]\n",
    "$$ The first term in this expression is diagonal and\n",
    "does not flip spins. This is the so-called Ising term. The second term\n",
    "is off-diagonal, and involves lowering and raising operators on\n",
    "neighboring spins, and is responsible for flipping anti-aligned spins.\n",
    "This is the “$XY$” part of the Hamiltonian.\n",
    "\n",
    "The Heisenberg spin chain is a paradigmatic model in condensed matter.\n",
    "Not only it is attractive due to its relative simplicity, but can also\n",
    "describe real materials that can be studied experimentally. The\n",
    "Heisenberg chain is also a prototypical integrable system, that can be\n",
    "solved exactly by the Bethe Ansatz, and can be studied using\n",
    "bosonization techniques and conformal field theory.\n",
    "\n",
    "The Heisenberg spin chain is a paradigmatic model in condensed matter.\n",
    "Not only it is attractive due to its relative simplicity, but can also\n",
    "describe real materials that can be studied experimentally. The\n",
    "Heisenberg chain is also a prototypical integrable system, that can be\n",
    "solved exactly by the Bethe Ansatz, and can be studied using\n",
    "bosonization techniques and conformal field theory.\n",
    "\n",
    "In these lectures, we will be interested in obtaining its ground state\n",
    "properties of this model by numerically solving the time-independent\n",
    "Schrödinger equation: $$\\hat{H}|\\Psi\\rangle = E|\\Psi\\rangle,$$ where $H$\n",
    "is the Hamiltonian of the problem, $|\\Psi\\rangle$ its eigenstates, with\n",
    "the corresponding eigenvalues, or energies $E$.\n",
    "\n",
    "Exact diagonalization\n",
    "=====================\n",
    "\n",
    "![blocks](figures/fig_blocks.png)\n",
    "#### Pictorial representation of the Hamiltonian building recursion explained in the text. At each step, the block size is increased by adding a spin at a time.\n",
    "\n",
    "In this section we introduce a technique that will allow us to calculate\n",
    "the ground state, and even excited states of small Heisenberg chains.\n",
    "Exact Diagonalization (ED) is a conceptually simple technique which\n",
    "basically consists of diagonalizing the Hamiltonian matrix by brute\n",
    "force. Same as for the spin operators, the Hamiltonian also has a\n",
    "corresponding matrix representation. In principle, if we are able to\n",
    "compute all the matrix elements, we can use a linear algebra package to\n",
    "diagonalize it and obtain all the eigenvalues and eigenvectors @lapack.\n",
    "\n",
    "In these lectures we are going to follow a quite unconventional\n",
    "procedure to describe how this technique works. It is important to point out that this\n",
    "is a quite inefficient and impractical way to diagonalize the\n",
    "Hamiltonian, and more sophisticated techniques are generally used in\n",
    "practice.\n",
    "\n",
    "Two-spin problem\n",
    "----------------\n",
    "\n",
    "The Hilbert space for the two-spin problem consists of four possible\n",
    "configurations of two spins\n",
    "$$\\left\\{ |\\uparrow\\uparrow\\rangle,|\\uparrow\\downarrow\\rangle,|\\downarrow\\uparrow\\rangle,|\\downarrow\\downarrow\\rangle \\right\\}$$\n",
    "\n",
    "The problem is described by the Hamiltonian: $$\\hat{H}=\n",
    "\\hat{S}^z_1 \\hat{S}^z_2 +\n",
    "\\frac{1}{2}\\left[\n",
    "\\hat{S}^+_1 \\hat{S}^-_2 +\n",
    "\\hat{S}^-_1 \\hat{S}^+_2\n",
    "\\right]$$\n",
    "\n",
    "The corresponding matrix will have dimensions $4 \\times 4$. In order to\n",
    "compute this matrix we shall use some simple matrix algebra to first\n",
    "obtain the single-site operators in the expanded Hilbert space. This is\n",
    "done by following the following simple scheme: And operator $O_1$ acting\n",
    "on the left spin, will have the following $4 \\times 4$ matrix form:\n",
    "$$\\tilde{O}_1 = O_1 \\otimes {1}_2\n",
    "$$ Similarly, for an operator $O_2$ acting on the right\n",
    "spin: $$\\tilde{O}_2 = {1}_2 \\otimes O_2\n",
    "$$ where we introduced the $n \\times n$ identity matrix\n",
    "${1}_n$. The product of two operators acting on different sites\n",
    "can be obtained as: $$\\tilde{O}_{12} = O_1 \\otimes O_2$$ It is easy to\n",
    "see that the Hamiltonian matrix will be given by: $$H_{12}=\n",
    "S^z \\otimes S^z +\n",
    "\\frac{1}{2}\\left[\n",
    "S^+ \\otimes S^- +\n",
    "S^- \\otimes S^+\n",
    "\\right]$$ where we used the single spin ($2 \\times 2$) matrices $S^z$\n",
    "and $S^\\pm$. We leave as an exercise for the reader to show that the\n",
    "final form of the matrix is:\n",
    "\n",
    "$$H_{12} = \\left(\n",
    "\\begin{array}{cccc}\n",
    "1/4 & 0 & 0 & 0 \\\\\n",
    "0 & -1/4 & 1/2 & 0 \\\\\n",
    "0 & 1/2 & -1/4 & 0 \\\\\n",
    "0 & 0 & 0 & 1/4 \\\\\n",
    "\\end{array}\n",
    "\\right),\n",
    "$$\n",
    "\n",
    "Obtaining the eigenvalues and eigenvectors is also a straightforward\n",
    "exercise: two of them are already given, and the entire problem now\n",
    "reduces to diagonalizing a two by two matrix. We therefore obtain the\n",
    "well known result: The ground state\n",
    "$|s\\rangle = 1/\\sqrt{2}\\left[ |\\uparrow\\downarrow\\rangle - |\\downarrow\\uparrow\\rangle \\right]$,\n",
    "has energy $E_s=-3/4$, and the other three eigenstates\n",
    "$\\left\\{|\\uparrow\\uparrow\\rangle,|\\downarrow\\downarrow\\rangle,1/\\sqrt{2}\\left[ |\\uparrow\\downarrow\\rangle + |\\downarrow\\uparrow\\rangle \\right] \\right\\}$\n",
    "form a multiplet with energy $E_t=1/4$.\n",
    "\n",
    "Many spins\n",
    "----------\n",
    "\n",
    "Imagine now that we add a third spin to the right of our two spins. We\n",
    "can use the previous result to obtain the new $8 \\times 8$ Hamiltonian\n",
    "matrix as: $$H_{3}=\n",
    "H_{2} \\otimes {1}_2 +\n",
    "\\tilde{S}^z_2 \\otimes S^z +\n",
    "\\frac{1}{2}\\left[\n",
    "\\tilde{S}^+_2 \\otimes S^- +\n",
    "\\tilde{S}^-_2 \\otimes S^+\n",
    "\\right]$$ Here we used the single spin $S^z_1$, $S^\\pm_1$, and the\n",
    "\\`tilde\\` matrices defined in Eqs.(\\[tildeL\\]) and (\\[tildeR\\]):\n",
    "$$\\tilde{S}^z_2 = {1}_2 \\otimes S^z,$$ and\n",
    "$$\\tilde{S}^\\pm_2 = {1}_2 \\otimes S^\\pm,$$\n",
    "\n",
    "It is easy to see that this leads to a recursion scheme to construct the\n",
    "$2^i \\times 2^i$ Hamiltonian matrix the $i^\\mathrm{th}$ step as:\n",
    "\n",
    "$$H_{i}=\n",
    "H_{i-1} \\otimes {1}_2 +\n",
    "\\tilde{S}^z_{i-1} \\otimes S^z +\n",
    "\\frac{1}{2}\\left[\n",
    "\\tilde{S}^+_{i-1} \\otimes S^- +\n",
    "\\tilde{S}^-_{i-1} \\otimes S^+\n",
    "\\right]$$\n",
    "\n",
    "with $$\\tilde{S}^z_{i-1} = {1}_{2^{i-2}} \\otimes S^z,$$ and\n",
    "$$\\tilde{S}^\\pm_{i-1} = {1}_{2^{i-2}} \\otimes S^\\pm,$$\n",
    "\n",
    "This recursion algorithm can be visualized as a left ‘block’, to which\n",
    "we add new ‘sites’ or spins to the right, one at a time, as shown in\n",
    "Fig.\\[fig:block\\].The block has a ‘block Hamiltonian’, $H_L$, that is\n",
    "iteratively built by connecting to the new spins through the\n",
    "corresponding interaction terms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ADDING SITE  1  DIML=  1\n",
      "ADDING SITE  2  DIML=  2\n",
      "ADDING SITE  3  DIML=  4\n",
      "ADDING SITE  4  DIML=  8\n",
      "ADDING SITE  5  DIML=  16\n",
      "ADDING SITE  6  DIML=  32\n",
      "ADDING SITE  7  DIML=  64\n",
      "[-3.3749326  -2.98224049 -2.98224049 -2.98224049 -2.50372907 -2.50372907\n",
      " -2.50372907 -2.33380386 -2.10812936 -2.10812936 -2.10812936 -2.0567576\n",
      " -2.0567576  -2.0567576  -1.91720711 -1.84883362 -1.84883362 -1.84883362\n",
      " -1.8316135  -1.8316135  -1.8316135  -1.8316135  -1.8316135  -1.67470637\n",
      " -1.67470637 -1.67470637 -1.64520463 -1.46566636 -1.46566636 -1.46566636\n",
      " -1.46566636 -1.46566636 -1.42441526 -1.42441526 -1.42441526 -1.39501087\n",
      " -1.22827398 -1.22827398 -1.22827398 -1.22563298 -1.22563298 -1.22563298\n",
      " -1.22563298 -1.22563298 -1.12098189 -1.03810412 -1.03810412 -1.03810412\n",
      " -1.03810412 -1.03810412 -0.99669162 -0.99669162 -0.99669162 -0.98276266\n",
      " -0.98276266 -0.98276266 -0.80242001 -0.80242001 -0.80242001 -0.80242001\n",
      " -0.80242001 -0.72918536 -0.72918536 -0.72918536 -0.72745015 -0.72745015\n",
      " -0.72745015 -0.6984272  -0.62176807 -0.62176807 -0.62176807 -0.62176807\n",
      " -0.62176807 -0.62128849 -0.61664185 -0.61664185 -0.61664185 -0.48599396\n",
      " -0.48599396 -0.48599396 -0.4510669  -0.4510669  -0.4510669  -0.4510669\n",
      " -0.4510669  -0.44720335 -0.44720335 -0.44720335 -0.39170281 -0.39170281\n",
      " -0.39170281 -0.39170281 -0.39170281 -0.37901811 -0.31675346 -0.31675346\n",
      " -0.31675346 -0.30093324 -0.30093324 -0.30093324 -0.30093324 -0.30093324\n",
      " -0.2170348  -0.2170348  -0.2170348  -0.17387953 -0.17387953 -0.17387953\n",
      " -0.17387953 -0.17387953 -0.17387953 -0.17387953 -0.12025775 -0.12025775\n",
      " -0.12025775 -0.07856615 -0.07856615 -0.07856615 -0.07856615 -0.07856615\n",
      " -0.04799743 -0.04799743 -0.04799743 -0.04799743 -0.04799743 -0.0155351\n",
      "  0.04289322  0.04289322  0.04289322  0.04289322  0.04289322  0.04289322\n",
      "  0.04289322  0.12769422  0.12769422  0.12769422  0.16805495  0.16805495\n",
      "  0.16805495  0.25423037  0.25423037  0.25423037  0.25423037  0.25423037\n",
      "  0.3208237   0.3208237   0.3208237   0.34784113  0.35531584  0.35531584\n",
      "  0.35531584  0.35531584  0.35531584  0.36731657  0.36731657  0.36731657\n",
      "  0.36731657  0.36731657  0.36731657  0.36731657  0.41720711  0.53535109\n",
      "  0.53535109  0.53535109  0.59949643  0.59949643  0.59949643  0.64679072\n",
      "  0.64679072  0.64679072  0.64679072  0.64679072  0.75        0.75\n",
      "  0.75        0.75        0.75        0.75        0.75        0.777125\n",
      "  0.777125    0.777125    0.81450906  0.81450906  0.81450906  0.81450906\n",
      "  0.81450906  0.85384723  0.95091982  0.95091982  0.95091982  0.97628634\n",
      "  0.97628634  0.97628634  0.97628634  0.97628634  1.039326    1.039326\n",
      "  1.039326    1.039326    1.039326    1.13268343  1.13268343  1.13268343\n",
      "  1.13268343  1.13268343  1.13268343  1.13268343  1.17046791  1.17046791\n",
      "  1.17046791  1.18631931  1.18631931  1.18631931  1.18631931  1.18631931\n",
      "  1.32342873  1.32342873  1.32342873  1.38251439  1.39334097  1.39334097\n",
      "  1.39334097  1.39334097  1.39334097  1.45710678  1.45710678  1.45710678\n",
      "  1.45710678  1.45710678  1.45710678  1.45710678  1.49369887  1.49369887\n",
      "  1.49369887  1.58935298  1.58935298  1.58935298  1.58935298  1.58935298\n",
      "  1.67387953  1.67387953  1.67387953  1.67387953  1.67387953  1.67387953\n",
      "  1.67387953  1.75        1.75        1.75        1.75        1.75\n",
      "  1.75        1.75        1.75        1.75      ]\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "\n",
    "# PARAMETERS\n",
    "nsites = 8\n",
    "\n",
    "#Single site operators\n",
    "sz0 = np.zeros(shape=(2,2)) # single site Sz\n",
    "splus0 = np.zeros(shape=(2,2)) # single site S+\n",
    "sz0[0,0]         = -0.5\n",
    "sz0[1,1]         =  0.5\n",
    "splus0[1,0]      =  1.0\n",
    "\n",
    "term_szsz = np.zeros(shape=(4,4)) #auxiliary matrix to store Sz.Sz\n",
    "term_szsz = np.kron(sz0,sz0)\n",
    "\n",
    "term_spsm = np.zeros(shape=(4,4)) #auxiliary matrix to store 1/2 S+.S-\n",
    "term_spsm = np.kron(splus0,np.transpose(splus0))*0.5\n",
    "term_spsm += np.transpose(term_spsm)\n",
    "\n",
    "h12 = term_szsz+term_spsm\n",
    "\n",
    "H = np.zeros(shape=(2,2))\n",
    "for i in range(1,nsites):\n",
    "\n",
    "    diml = 2**(i-1)        # 2^i\n",
    "    dim = diml*2\n",
    "    \n",
    "    print (\"ADDING SITE \",i,\" DIML= \",diml)\n",
    "\n",
    "    Ileft = np.eye(diml)\n",
    "    Iright = np.eye(2)\n",
    "\n",
    "    # We add the term for the interaction S_i.S_{i+1}\n",
    "\n",
    "    aux = np.zeros(shape=(dim,dim))\n",
    "    aux = np.kron(H,Iright)\n",
    "    H = aux\n",
    "\n",
    "    H = H + np.kron(Ileft,h12)\n",
    "\n",
    "w, v = np.linalg.eigh(H) #Diagonalize the matrix\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1 -0.13440094427539903\n",
      "0.2 -0.2745106220291552\n",
      "0.30000000000000004 -0.4192870911355526\n",
      "0.4 -0.5675111660216381\n",
      "0.5 -0.7178307681913986\n",
      "0.6 -0.8688168893699839\n",
      "0.7000000000000001 -1.0190260819336356\n",
      "0.8 -1.167063391675081\n",
      "0.9 -1.3116396792650657\n",
      "1.0 -1.4516183124034265\n",
      "1.1 -1.5860479724220509\n",
      "1.2000000000000002 -1.714180378979959\n",
      "1.3000000000000003 -1.8354736576118653\n",
      "1.4000000000000001 -1.949583532313432\n",
      "1.5000000000000002 -2.056345366485309\n",
      "1.6 -2.155750314045565\n",
      "1.7000000000000002 -2.2479186067208534\n",
      "1.8000000000000003 -2.333072470791156\n",
      "1.9000000000000001 -2.4115105081323973\n",
      "2.0 -2.483584724611145\n",
      "2.1 -2.549680826427394\n",
      "2.2 -2.610201969519424\n",
      "2.3000000000000003 -2.6655558432015374\n",
      "2.4000000000000004 -2.716144781467402\n",
      "2.5000000000000004 -2.7623584995787276\n",
      "2.6 -2.8045690239636674\n",
      "2.7 -2.8431273972140207\n",
      "2.8000000000000003 -2.878361778912821\n",
      "2.9000000000000004 -2.910576613895644\n",
      "3.0000000000000004 -2.940052593457619\n",
      "3.1 -2.9670471865212606\n",
      "3.2 -2.9917955639500233\n",
      "3.3000000000000003 -3.0145117788059745\n",
      "3.4000000000000004 -3.035390098252337\n",
      "3.5000000000000004 -3.0546064094468397\n",
      "3.6 -3.0723196429088726\n",
      "3.7 -3.0886731733195503\n",
      "3.8000000000000003 -3.103796170357368\n",
      "3.9000000000000004 -3.1178048817330026\n",
      "4.0 -3.1308038377109333\n",
      "4.1 -3.142886971630536\n",
      "4.2 -3.154138654702368\n",
      "4.3 -3.16463464600466\n",
      "4.3999999999999995 -3.174442960414687\n",
      "4.5 -3.1836246583926866\n",
      "4.6 -3.1922345622559467\n",
      "4.7 -3.2003219039634057\n",
      "4.8 -3.2079309095711137\n",
      "4.9 -3.215101325488383\n",
      "5.0 -3.221868891516565\n",
      "5.1 -3.2282657654275697\n",
      "5.2 -3.2343209035664526\n",
      "5.3 -3.2400604016634493\n",
      "5.4 -3.2455077997309707\n",
      "5.5 -3.250684354611485\n",
      "5.6 -3.25560928344012\n",
      "5.7 -3.2602999809967206\n",
      "5.8 -3.2647722136486683\n",
      "5.9 -3.269040292330319\n",
      "6.0 -3.273117226767761\n",
      "6.1 -3.2770148629393656\n",
      "6.2 -3.28074400556245\n",
      "6.3 -3.2843145272139194\n",
      "6.4 -3.287735465526885\n",
      "6.5 -3.2910151097550813\n",
      "6.6 -3.294161077861073\n",
      "6.7 -3.2971803851619477\n",
      "6.8 -3.3000795054560577\n",
      "6.9 -3.3028644254554345\n",
      "7.0 -3.3055406932597697\n",
      "7.1 -3.3081134615283663\n",
      "7.2 -3.31058752593522\n",
      "7.3 -3.312967359428785\n",
      "7.4 -3.3152571427610007\n",
      "7.5 -3.3174607916994545\n",
      "7.6 -3.319581981290975\n",
      "7.7 -3.3216241675047438\n",
      "7.8 -3.3235906065466003\n",
      "7.9 -3.325484372104261\n",
      "8.0 -3.3273083707543263\n",
      "8.1 -3.329065355736452\n",
      "8.2 -3.330757939277211\n",
      "8.3 -3.3323886036259625\n",
      "8.4 -3.3339597109468744\n",
      "8.5 -3.335473512195249\n",
      "8.6 -3.3369321550919087\n",
      "8.7 -3.338337691296708\n",
      "8.8 -3.339692082870839\n",
      "8.9 -3.3409972081075736\n",
      "9.0 -3.342254866802048\n",
      "9.1 -3.3434667850227884\n",
      "9.2 -3.344634619440467\n",
      "9.3 -3.3457599612632443\n",
      "9.4 -3.3468443398222694\n",
      "9.5 -3.347889225846049\n",
      "9.6 -3.3488960344579564\n",
      "9.700000000000001 -3.349866127927187\n",
      "9.8 -3.3508008182000277\n",
      "9.9 -3.351701369235202\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHipJREFUeJzt3XuUVeWZ5/HvU3fqyp2iChAQRK6iIio6icbLYEysTtLpmFmdmFuznGkn6V69VsY0nTFjJjN2Mqs7mdbphM61Rzt2xk6QKN41MaZbRRTkDiUgVEFxp6iiqMup88wf51RxKA5UwbnsU2f/Pmux6pxdb539nMT1/vZ+33fvbe6OiIiET0HQBYiISDAUACIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkEpLAJjZUjPbZmaNZnZ/kt+Xmtk/x3//hplNTcd+RUTk4hWl+gFmVgg8AtwGNAFrzGyVu29OaPZF4Ji7zzCzu4G/Bj41yOfqCjURkYvg7jaUduk4A1gMNLr7TnfvBh4HGga0aQB+Fn/9BHCLmQ2pQBERyYyUzwCAemBvwvsm4NpztXH3iJm1AmOAw4N9uG5VISIyNBd6XJ1zk8BmtszM3gq6DhGRfJeOAGgGJie8nxTflrSNmRUBNcCRZB/m7ivcfVEa6hIRkfNIRwCsAWaa2TQzKwHuBlYNaLMKuCf++g+Bl11jOyIigUp5DiA+pn8f8BxQCPzY3TeZ2YPAW+6+CvgR8H/NrBE4SiwkREQkQJarB+J9y0BztT4RkVzTNwmczWWgIiIyDCkARERCKh3XAYiISBLRqHOyO0J7V4T2zggnOk+/buvsob0rQltn7F9715nvV3zmasZXl2W0PgWAiMgA7k5Hd29Ch9yT0HFHaBvYicfbtSd06u2dEdq7I1zsNOaxjh4FgIjIULk7XZEoJzp7Yh3wGUfY8Q47/vr00XhPwu9Pd+rRNK0/KS8ppLK0iKqyIirLiqnqe11aRGVZEVXxbZXxbVVlsX+TR49ITwHnoQAQkZzQFek9feTcFTmrEz/v0XhXT//fRtLUc5cVF1BZWkx12Zmdc2VpcX8nnawT79teVVpMRWkhRYW5O9WqABCRlPT0RhM6454zOvG2zp6E4ZIknXjCcEl3bzQt9ZQUFZw+yu7vuAd20MVUlhXFOvfShDYJnXpxDnfc6aIAEAmpSG+Uk129sSPtroQOu/PMjnmwo/GuSHo67qIC6++0q+Id9Jkd+ZlH3lVlxWcMmfR13KVFhWmpJwwUACLDTG98ZcnpDronNp49YJy7bUAnPvBo/FRPb1rqKSywhOGRvg65OGF4JNaR93fafe/LzuzES4sKLvhulpIaBYBIlkSjTkdPb//qkTOHRnrOmKxsTxxOGdCRt3dF0lKPGVSWFlE9oLM+qxM/4yi7OOEoPdamrFgd93ClABAZou5IlKZjHf1H28k68f7VJUmOxtu7Ln5J4ECJHXPfOHd1wtH0GdvKknfi5SWF6rhDTgEgMkA06uw91sHWlja2t7Sx9UDs567DJ1NeYVJeUpgwXl181iTk2ePeZ09OVpYUUVCgjltSpwCQ0HJ3DrV1se1AG9taYv+2H2hj+4H2pOPjZjB59AhGlZecc0ng+TrxipLcXhIo4aMAkFA40dnD9pY2tsWP5rfGO/tjHT1J20+oLmVWbTWzJlRy2YQqLq+tZsb4SkaUaIWJ5A8FgOSVzp5e3jvUzvYDbf1DONta2tjX2pm0fVVZEZfXVsU7+djPWbVVjCwvyXLlItmnAJBhqTfq7DnawbaWE2xraWfbgRNsa2lj95EOepOM05cUFTBzfCWzaquYFe/kZ9VWUVtdpolQCS0FgOQ0d+fAib5x+lhnv/1AGzsOttHZc/YFSAUG08dVnO7kJ1RxWW0VU8dUUKiJU5EzKAAkZ7R29LD94JlDN9sOtNF6Kvk4/cSastOdfLzDnzG+krJijdOLDIUCQLKus6eXxoPt/R183wqclhPJx+lrRhSfNXRz2YQqakYUZ7lykfyiAJCMOtTWxZrdR89YZrn7yMmkt9otKy5g5vgzh24ur61ifFWpxulFMkABIGnX3hXhuY0trFzXzO8bD5/V2RcWGDPi4/SXJRzVTxldrnF6kSxSAEhadEeivLr9ECvXNfPilgP9E7RFBcYNM8Ywt66GWbWVzJpQzfRxFRqnF8kBKQWAmY0G/hmYCuwG/sjdjyVp1wtsiL/d4+53pbJfyQ3RqLN2zzFWvtPM0xv2czzhoqprpo6iYWE9d86fyKgKrakXyUXmKdydysy+DRx194fM7H5glLv/lyTt2t298gI/2yG2DFByy/YDbax8p5kn1+2j+fip/u2XTaikYWE9DQvrmDSqPMAKRcKpb67M3Yc0lppqAGwDbnL3/WY2EfiNu89K0k4BMMztbz3FqnX7WLluH1v2n+jfPrGmjLuuqKNhYT2zJ1ZpslYkQNkOgOPuPjL+2oBjfe8HtIsA64AI8JC7rxzCZysAAtba0cMzG/ezcl0zb+w62n8r4+qyIu5cMJG7rqjn2mmjdWdKkRxxoQEw6ByAmb0I1Cb51fLEN+7ufZ12Epe4e7OZTQdeNrMN7v7eOfa3DFg2WF2SGZ09vbyy9SAr1zXzytZD/c9pLSkq4NbZ42lYWM9Ns8bpsXsieSArQ0AD/uanwFPu/sQg7XQGkCW9Uef1nUdY+U4zz25soS3+xCkzWHLpGBoW1rN0Xi3VZbrwSiSXpf0MYBCrgHuAh+I/n0xS0Cigw927zGwscAPw7RT3K2mwsbmVle808+t393HgRFf/9nn11fzBwno+ekUdE6rLAqxQRDIp1TOAMcAvgCnA+8SWgR41s0XAve7+JTNbAvwAiAIFwHfd/UdD+GydAWTI3qMdfPOpzTy/+UD/timjy2lYWEfDwjpmjK8KsDoRuVhZnQTOJAVA+nVFevmHV3fy8CuNdPZEqSgp5BNXT6JhYT1XTRmpFTwiw1y2h4BkmPjNtoN8Y9Umdh/pAOCjV9Sx/MOzqa3REI9IWCkA8lzTsdhwz3ObYsM9M8ZX8mDDXJZcOjbgykQkaAqAPNUV6eWHv9vF3728g86eKOUlhfzZrTP53JJplBTpweQiogDIS7/dfohvrNrErsMnAfjIgoksv3M2E2tGBFyZiOQSBUAeaT5+iv/+1Gae2dgCwKXjKniwYR43zNBwj4icTQGQB/qGex5+uZFTPb2UlxTylVtm8vkbNNwjIuemABjmXo0P9+yMD/fcuWAif6XhHhEZAgXAMNUdibL8Vxv4f2ubAJg+roIH75rHjTM13CMiQ6MAGIZOdkW499G1/G7HYcqKC/jKLZfxxRs13CMiF0YBMMwcae/iCz9dw/qmVsZUlPCTz1/Dgkln3YFbRGRQCoBhZO/RDj774zfZdfgkk0eP4B+/cC3TxlYEXZaIDFMKgGFi874T3POTNznU1sXsidX87AvXML5Kt3EQkYunABgGXt95hD/52Vu0dUW4bvpoVnx2ke7NLyIpUwDkuGc37ufLj6+jOxLlw/Nr+dtPLdTTuEQkLRQAOezR19/n609uxB0+e/0lPPDRuRTq+bsikiYKgBzk7nz3xR1876UdAPzFbZdx34dm6H79IpJWCoAc4+58/cmNPPr6HgoM/sfH5nP34ilBlyUieUgBkGP+4Xc7efT1PZQWFfB3n76S2+fWBl2SiOQpXTqaQ37feJiHntkKwPfuVucvIpmlAMgRTcc6uO+f3ibqcN/NM1g6T52/iGSWAiAHdPb0cu+jaznW0cMHLxvHn992WdAliUgIKAAC5u785a82sLH5BFNGl/O/775SSz1FJCtSCgAz+6SZbTKzqJktOk+7pWa2zcwazez+VPaZb/7x397nl283M6K4kB985mpqynWFr4hkR6pnABuBjwOvnquBmRUCjwB3AHOAT5vZnBT3mxfe3HWUbz61GYC//sMFzJ5YHXBFIhImKS0DdfctwGAXKC0GGt19Z7zt40ADsDmVfQ93La2d/KfH3iYSdb504zTuuqIu6JJEJGSyMQdQD+xNeN8U35aUmS0zs7cyXlWAuiK9/MfH1nK4vYvrp4/h/jsuD7okEQmhQc8AzOxFINmaxOXu/mS6C3L3FcAKM/N0f3aueOSV93hnz3Hqasp4+D9cSVGh5uJFJPsGDQB3vzXFfTQDkxPeT4pvC6Xdh0/y/d++B8DffmohYypLA65IRMIqG4eea4CZZjbNzEqAu4FVWdhvznF3vvHrTXRHonz8qnqunT4m6JJEJMRSXQb6MTNrAq4Hnjaz5+Lb68xsNYC7R4D7gOeALcAv3H1TamUPT89vPsBvth2iqqyIr90xO+hyRCTkzD03h9r75gBytb4L1dEd4ba/eZXm46f4b3fN5Z4lU4MuSUTyTN+KTHcf0tWkmn3MkkdeaaT5+Cnm1lXzx9ddEnQ5IiIKgGzYeaidFa/uBODBhnm61YOI5AQFQIa5Ow+s2kRPr/OpRZO5+pJRQZckIgIoADLumY0t/G7HYWpGFPPVpbOCLkdEpJ8CIIO6I1G+9fQWAL66dJbW/ItITlEAZNDKd5ppPn6KGeMrufsaPddXRHKLAiBDeqPO38ev+P3Tmy/VxK+I5BwFQIY8vWE/uw6fZMrocj66QHf6FJHcowDIgGjU+T+vNAJw7wcv1c3eRCQnqWfKgJe2HmRrSxu11WV84upz3vlaRCRQCoA0c3cejh/9/8kHplNaVBhwRSIiySkA0uz3jUdYv/c4oytK+PTiyYP/gYhIQBQAafbwKzsA+OKN0ygvSemJmyIiGaUASKMNTa28vvMoVWVFfOZ63fBNRHKbAiCN/unNPQB88urJVJcVB1yNiMj5KQDS5GRXhFXrYk+61Ni/iAwHCoA0+fX6fZzs7uWaqaOYOaEq6HJERAalAEiTn8eHf3TPHxEZLhQAabBpXyvrm1qpLivizgUTgy5HRGRIFABp8PibewH4+FWTKCvWhV8iMjwoAFLU0R1h5Tuxyd+7NfkrIsNISgFgZp80s01mFjWzRedpt9vMNpjZOjN7K5V95pqn391PW1eEK6eM5PLa6qDLEREZslQvVd0IfBz4wRDa3uzuh1PcX85ZGV/6+alFOvoXkeElpQBw9y0AZuF82Mmhti7+7b0jFBcad8zT5K+IDC/ZmgNw4HkzW2tmy7K0z4x7dlMLUYd/N3McNeW68ldEhpdBzwDM7EWgNsmvlrv7k0Pcz43u3mxm44EXzGyru796jv0tA4ZFSDy1fh8Ad87X0b+IDD+DBoC735rqTty9Of7zoJn9ClgMJA0Ad18BrDAzT3W/mXTwRCdv7j5KSWEBt82dEHQ5IiIXLONDQGZWYWZVfa+B24lNHg9rqzfsxx0+OGucbvwmIsNSqstAP2ZmTcD1wNNm9lx8e52ZrY43mwC8ZmbrgTeBp9392VT2mwue3rAfgI/oyl8RGabMPTdHWvqGgHKxvv2tp7j+f75MaVEBa79+G5WlevCLiASvb0Wmuw9paaauBL4Iz286AMDNs8ar8xeRYUsBcBFe2noQgNvmaPJXRIYvBcAFOtkV4fX3jmAGN80aF3Q5IiIXTQFwgV5rPEx3b5QrJ49kTGVp0OWIiFw0BcAFenlLbPjnltka/hGR4U0BcAGiUeeVbbEAuHnW+ICrERFJjQLgAmzad4KDbV1MrClj9kQ991dEhjcFwAV4aWts+eeHLh8f2jugikj+UABcgN9uPwRo+EdE8oMCYIhOdPawfu9xigqM6y4dE3Q5IiIpUwAM0Rs7jxJ1uHLKSF39KyJ5QQEwRL9vjD3NcsmlYwOuREQkPRQAQ/RaPABunKkAEJH8oAAYgpbWThoPtlNRUsjCySODLkdEJC0UAEPwr+/Fjv6vnT6G4kL9TyYi+UG92RD8vvEIAEu0+kdE8ogCYAje2BULgOumKwBEJH8oAAax7/gpmo6doqqsiNkTq4MuR0QkbRQAg1iz+ygAiy4ZRWGBbv8gIvlDATCIN3bFAmDxNA3/iEh+UQAM4s3+ABgVcCUiIumlADiPI+1dNB5sp7SogPn1Wv8vIvklpQAws++Y2VYze9fMfmVmSXtJM1tqZtvMrNHM7k9ln9m0ZvcxIHb/n5IiZaWI5JdUe7UXgHnuvgDYDnxtYAMzKwQeAe4A5gCfNrM5Ke43K97eEwuAa6aODrgSEZH0SykA3P15d4/E374OTErSbDHQ6O473b0beBxoSGW/2fL2+7EAuGqKxv9FJP+kc1zjC8AzSbbXA3sT3jfFt+W07kiUd5tbgdgQkIhIvhn0xvZm9iJQm+RXy939yXib5UAEeCzVgsxsGbAs1c9J1eb9J+iORJk+roKR5SVBlyMiknaDBoC733q+35vZ54CPALe4uydp0gxMTng/Kb7tXPtbAawws2SflTV9wz9Xa/hHRPJUqquAlgJfBe5y945zNFsDzDSzaWZWAtwNrEplv9mwNj4BfNUlCgARyU+pzgE8DFQBL5jZOjP7PoCZ1ZnZaoD4JPF9wHPAFuAX7r4pxf1m3Lo9xwGN/4tI/krp4bbuPuMc2/cBH054vxpYncq+sulwexfNx09RXlLIzPFVQZcjIpIRuropiXebYkf/8+trdAM4EclbCoAk1u2NLf/U4x9FJJ8pAJJYvzd2BrBgkgJARPKXAmAAd+8fArpick3A1YiIZI4CYICmY6c41tHDmIoS6keOCLocEZGMUQAM8G5TbPx//qQazDQBLCL5SwEwwIb4/X/m12v4R0TymwJggI3xAJinABCRPKcASODuOgMQkdBQACRoOnaK1lM9jK0sYWJNWdDliIhklAIgQd/wz9w6TQCLSP5TACTYtO8EAPPqqwOuREQk8xQACTbuO30GICKS7xQACfrOAObW6QxARPKfAiDuYFsnh9q6qCotYvKo8qDLERHJOAVAXN/R/+y6agp0C2gRCQEFQNyW/bEAmDNRwz8iEg4KgLjN8TOAORr/F5GQUADEbdYZgIiEjAIAONXdy+7DJykqMGZOqAy6HBGRrFAAAFtbThB1mDG+ktKiwqDLERHJCgUAsLWlDYDZGv4RkRApSuWPzew7wEeBbuA94PPufjxJu91AG9ALRNx9USr7Tbe+FUCzJ1YFXImISPakegbwAjDP3RcA24Gvnaftze6+MNc6fzgdAJfX6gxARMIjpQBw9+fdPRJ/+zowKfWSssvd2bpfQ0AiEj7pnAP4AvDMOX7nwPNmttbMlqVxnylrOnaKtq4IYytLGVdVGnQ5IiJZM+gcgJm9CNQm+dVyd38y3mY5EAEeO8fH3OjuzWY2HnjBzLa6+6vn2N8yIGshcXoCWOP/IhIugwaAu996vt+b2eeAjwC3uLuf4zOa4z8PmtmvgMVA0gBw9xXACjNL+lnptrV//F8BICLhktIQkJktBb4K3OXuHedoU2FmVX2vgduBjansN536zgBmaQJYREIm1TmAh4EqYsM668zs+wBmVmdmq+NtJgCvmdl64E3gaXd/NsX9ps2WFi0BFZFwSuk6AHefcY7t+4APx1/vBK5IZT+Z0tkTuwVEYYExY7xuASEi4RLqK4F3HGgn6jB9bIVuASEioRPqAOgb/pmlCWARCaFQB4AuABORMAt3AGgCWERCLLQB4O5aAioioRbaADjU1sXRk91UlRVRV1MWdDkiIlkX2gDovwVEbTVmFnA1IiLZF+IAiN8CQuP/IhJS4Q2A/X3j/woAEQmn0AbAFj0GUkRCLpQB0NMbpfFg/Axggs4ARCScQhkA7x1qp6fXmTqmnIrSlG6HJCIybIUyAPQMYBGR0AaAxv9FREIaALoFhIhI6ALA3dm8ry8AdAYgIuEVugA41NbFkZPdVJcVMWnUiKDLEREJTOgCYFP86H9OnW4BISLhFsIAaAU0/CMiEsIAiJ0BzK2rCbgSEZFghTYA5tXrDEBEwi1UAdDa0cOeox2UFhUwY1xl0OWIiAQq5QAws2+a2btmts7MnjezunO0u8fMdsT/3ZPqfi/GxoTx/6LCUGWfiMhZ0tELfsfdF7j7QuAp4L8ObGBmo4EHgGuBxcADZjYqDfu+IBuaYwEwv17j/yIiKQeAu59IeFsBeJJm/x54wd2Puvsx4AVgaar7vlDvNh0HYMEkBYCISFpuhWlm3wI+C7QCNydpUg/sTXjfFN+W7LOWAcvSUddA6/fGzgAWTBqZiY8XERlWhnQGYGYvmtnGJP8aANx9ubtPBh4D7kulIHdf4e6LUvmMZA63d9F8/BTlJYXMGK8JYBGRIZ0BuPutQ/y8x4DVxMb7EzUDNyW8nwT8ZoifmRbr9sSGf+bX11BYoCuARUTSsQpoZsLbBmBrkmbPAbeb2aj45O/t8W1Z887eYwAsnKLhHxERSM8cwENmNguIAu8D9wKY2SLgXnf/krsfNbNvAmvif/Ogux9Nw76H7O33Y2cAV03J+uIjEZGcZO7JFu0Ez8wcYrdvTlVPb5QF33ieUz29vPVXtzK2sjTlzxQRyTV9N7h09yGNc4fiaqhN+05wqqeX6WMr1PmLiMSFIgDW7IqNNl0zdXTAlYiI5I5QBMDrO48AsHiaAkBEpE/eB0CkN8qb8TOA6y8dE3A1IiK5I+8DYH1TK21dEaaNraBupB4BKSLSJ+8D4Hc7DgFw44yxAVciIpJb8j4AfrMtFgAfuGxcwJWIiOSWvA6Aw+1drG86TklhATfM0Pi/iEiivA6Al7YcwB2WzBhDeUlabnwqIpI38joAnt3YAsDtc2oDrkREJPfkbQAc7+jmtcbDFBjcPndC0OWIiOScvA2A1Rta6Ol1llw6Vrd/EBFJIm8D4Im1sQeQfezKpA8eExEJvbwMgG0tbby95zgVJYUsnafxfxGRZPIyAH782i4APn7VJCpKtfpHRCSZvAuAnYfa+Ze3mygw+PwNU4MuR0QkZ+Xd4fGRk93U1pSx5NIxTB+nh7+LiJxLXj4RrLOnl65IlJoRxWmvS0QkV13oE8HyMgBERMJIj4QUEZEhUQCIiIRUSpPAZvZNoAGIAgeBz7n7viTteoEN8bd73P2uVPYrIiKpS2kOwMyq3f1E/PWXgTnufm+Sdu3ufkFLcjQHICJyYbI6B9DX+cdVAOqtRUSGiZSvAzCzbwGfBVqBm8/RrMzM3gIiwEPuvjLV/YqISGoGHQIysxeBZDfUWe7uTya0+xpQ5u4PJPmMendvNrPpwMvALe7+3jn2twxYBlw99K8hIiJ9sn4dgJlNAVa7+7xB2v0UeMrdnxiknYaTREQuwlADINVVQDPdfUf8bQOwNUmbUUCHu3eZ2VjgBuDbg332UL/AgH295e6LLvTvhrMwfmcI5/cO43eGcH7vbH3nVOcAHjKzWcSWgb4P3AtgZouAe939S8Bs4AdmFiU26fyQu29Ocb8iIpKilALA3T9xju1vAV+Kv/5XYH4q+xERkfTLtyuBVwRdQADC+J0hnN87jN8Zwvm9s/Kdc/ZmcCIikln5dgYgIiJDlBcBYGZLzWybmTWa2f1B15MNZjbZzF4xs81mtsnMvhJ0TdliZoVm9o6ZPRV0LdliZiPN7Akz22pmW8zs+qBryjQz+/P4f9sbzeznZlYWdE2ZYGY/NrODZrYxYdtoM3vBzHbEf47KxL6HfQCYWSHwCHAHMAf4tJnNCbaqrIgAf+Huc4DrgD8NyfcG+AqwJegisux7wLPufjlwBXn+/c2sHvgysCh+bVEhcHewVWXMT4GlA7bdD7zk7jOBl+Lv027YBwCwGGh0953u3g08TuyahLzm7vvd/e346zZiHUJ9sFVlnplNAu4Efhh0LdliZjXAB4AfAbh7t7sfD7aqrCgCRphZEVAOnHWn4Xzg7q8CRwdsbgB+Fn/9M+APMrHvfAiAemBvwvsmQtARJjKzqcCVwBvBVpIV3wW+Suzak7CYBhwCfhIf+vqhmVUEXVQmuXsz8L+APcB+oNXdnw+2qqya4O77469bgAmZ2Ek+BEComVkl8C/Anw24O2veMbOPAAfdfW3QtWRZEXAV8PfufiVwkgwNCeSK+Jh3A7HwqwMqzOyPg60qGB5bqpmR5Zr5EADNwOSE95Pi2/KemRUT6/wfc/dfBl1PFtwA3GVmu4kN9X3IzB4NtqSsaAKa3L3vDO8JYoGQz24Fdrn7IXfvAX4JLAm4pmw6YGYTAeI/D2ZiJ/kQAGuAmWY2zcxKiE0UrQq4poyz2JMffgRscfe/CbqebHD3r7n7JHefSuz/55fdPe+PCt29Bdgbv+0KwC1Avt9OZQ9wnZmVx/9bv4U8n/geYBVwT/z1PcCT52l70VJ+HkDQ3D1iZvcBzxFbKfBjd98UcFnZcAPwGWCDma2Lb/tLd18dYE2SOf8ZeCx+kLMT+HzA9WSUu79hZk8AbxNb8fYOeXpFsJn9HLgJGGtmTcADwEPAL8zsi8Tus/ZHGdm3rgQWEQmnfBgCEhGRi6AAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElAJARCSk/j/S6wenxajNdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "pyplot.rcParams['axes.linewidth'] = 2 #set the value globally\n",
    "%matplotlib inline\n",
    "\n",
    "Beta = np.arange(0.1,10,0.1)\n",
    "et = np.copy(Beta)\n",
    "n = 0\n",
    "for x in Beta:\n",
    "    p = np.exp(-w*x)\n",
    "    z = np.sum(p)\n",
    "    et[n] = np.dot(w,p)/z\n",
    "    print (x,et[n])    \n",
    "    n+=1\n",
    "\n",
    "pyplot.plot(1/Beta,et,lw=2);\n",
    "     \n",
    "    \n",
    "    \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge 12.1:\n",
    "\n",
    "Compute the energy and specific heat of the spin chain with $L=12$ as a function of temperature, for $T < 4$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A practical exact diagonalization algorithm\n",
    "\n",
    "1. Initialization: Topology of the lattice, neighbors, and signs. \n",
    "2. Contruction of a basis suitable for the problem. \n",
    "3. Constuction of the matrix element of the Hamiltonian. \n",
    "4. Diagonalization of the matrix. \n",
    "5. Calculation of observables or expectation values.\n",
    "\n",
    "As we shall see below, we are going to need the concept of “binary\n",
    "word”. A binary word is the binary representation of an integer (in\n",
    "powers of ‘two’, i.e. $n=\\sum_{i}b_{i}.2^{i}$), and consist of a\n",
    "sequence of ‘bits’. A bit $b_{i}$ in the binary basis correspond to our\n",
    "digits in the decimal system, and can assumme the values ‘zero’ or\n",
    "‘one’. At this stage we consider appropriate to introduce the logical\n",
    "operators `AND, OR, XOR`. The binary operators act between bits and\n",
    "their multiplication tables are listed below\n",
    "\n",
    "**AND** | 0 | 1\n",
    ":---:|:---:|:---:\n",
    "**0** | 0 | 0\n",
    "**1** | 0 | 1\n",
    "\n",
    "**OR** | 0 | 1\n",
    ":---:|:---:|:---:\n",
    "**0** | 0 | 1\n",
    "**1** | 1 | 1\n",
    "\n",
    "**XOR** | 0 | 1\n",
    ":---:|:---:|:---:\n",
    "**0** | 0 | 1\n",
    "**1** | 1 | 0\n",
    "\n",
    "\n",
    "The bit manipulation is a very useful tool, and most of the programming\n",
    "languages provide the needed commands.\n",
    "\n",
    "## Initilalization and definitions\n",
    "\n",
    "In the program, it is convenient to store in arrays all those quantities\n",
    "that will be used more frequently. In particular, we must determine the\n",
    "topology of the cluster, labeling the sites, and storing the components\n",
    "of the lattice vectors. We must also generate arrays with the nearest\n",
    "neighbors, and the next-nearest neighbors, according to our needs.\n",
    "\n",
    "## Construction of the basis\n",
    "\n",
    "Memory limitations impose severe restrictions on the size of the\n",
    "clusters that can be studied with this method. To understand this point,\n",
    "note that although the lowest energy state can be written in the\n",
    "$\\{|\\phi\n",
    "_{n}\\rangle \\}$ basis as\n",
    "$|\\psi _{0}\\rangle =\\sum_{m}c_{m}|\\phi _{m}\\rangle $, this expression is\n",
    "of no practical use unless $|\\phi _{m}\\rangle $ itself is expressed in a\n",
    "convenient basis to which the Hamiltonian can be easily applied. A\n",
    "natural orthonormal basis for fermion systems is the occupation number\n",
    "representation, describing all the possible distributions of $N_{e}$\n",
    "electrons over $N$ sites, while for spin systems it is covenient to work\n",
    "in a basis where the $S_{z}$ is defined at every site, schematically\n",
    "represented as $|n\\rangle =|\\uparrow \\downarrow \\uparrow ...\\rangle $.\n",
    "The size of this type of basis set\n",
    "grows exponentially with the system size. In practice this problem can be\n",
    "considerably alleviated by the use of symmetries of the Hamiltonian that\n",
    "reduces the matrix to a block-form. The most obvious symmetry is the\n",
    "number of particles in the problem which is usually conserved at least\n",
    "for fermionic problems. The total projection of the spin\n",
    "$S_{total}^{z}$, is also a good quantum number. For translationally\n",
    "invariant problems, the total momentum $\\mathbf{k%\n",
    "}$ of the system is also conserved introducing a reduction of $1/N$ in\n",
    "the number of states (this does not hold for models with open boundary\n",
    "conditions or explicit disorder). In addition, several Hamiltonians have\n",
    "additional symmetries. On a square lattice, rotations in $\\pi /2$ about\n",
    "a given site, spin inversion, and reflections with respect to the\n",
    "lattice axis are good quantum numbers (although care must be taken in\n",
    "their implementation since some of these operations are combinations of\n",
    "others and thus not independent). \n",
    "\n",
    "In the following we shall consider a spin-1/2 Heisenberg chain as a\n",
    "practical example. In this model it is useful to represent the spins\n",
    "pointing in the ‘up’ direction by a digit ‘1’, and the down-spins by a\n",
    "‘0’. Following this rule, a state in the $S^{z}$ basis can be\n",
    "represented by a sequence of ones and zeroes, i.e., a “binary word”.\n",
    "Thus, two Néel configurations in the 4-site chain can be seen as\n",
    "$$\\begin{aligned}\n",
    "\\mid \\uparrow \\downarrow \\uparrow \\downarrow \\rangle \\equiv |1010\\rangle , \\\\\n",
    "\\mid \\downarrow \\uparrow \\downarrow \\uparrow \\rangle \\equiv |0101\\rangle .\\end{aligned}$$\n",
    "Once the up-spins have been placed the whole configuration has been\n",
    "uniquely determined since the remaining sites can only be occupied by\n",
    "down-spins. The resulting binary number can be easily converted into\n",
    "integers $i\\equiv\n",
    "\\sum_{l}b_{l}.2^{l}$, where the summation is over all the sites of the\n",
    "lattice, and $b_{l}$ can be $1$ or $0$. For example:\n",
    "$$\\begin{array}{lllll}\n",
    "2^{3} & 2^{2} & 2^{1} & 2^{0} &  \\\\\n",
    "1 & 0 & 1 & 0 & \\rightarrow 2^{1}+2^{3}=10, \\\\\n",
    "0 & 1 & 0 & 1 & \\rightarrow 2^{0}+2^{2}=5.\n",
    "\\end{array}$$\n",
    "\n",
    "Using the above convention, we can construct the whole basis for the\n",
    "given problem. However, we must consider the memory limitations of our\n",
    "computer, introducing some symetries to make the problem more tractable.\n",
    "The symetries are operations that commute with the Hamiltonian, allowing\n",
    "us to divide the basis in subspaces with well defined quantum numbers.\n",
    "By means of similarity transformations we can generate a sequence of\n",
    "smaller matrices along the diagonal (i.e. the Hamiltonian matrix is\n",
    "“block diagonal”), and each of them can be diagonalized independently.\n",
    "For fermionic systems, the simplest symmetries are associated to the\n",
    "conservation of the number of particles and the projection\n",
    "$S_{total}^{z}$ of the total spin in the $z$ direction. In the spin-1/2\n",
    "Heisenberg model, a matrix with $2^{N}\\times 2^{N}$ elements can be\n",
    "reduced to $2S+1$ smaller matrices, corresponding to the projections $%\n",
    "m=(N_{\\uparrow }-N_{\\downarrow })/2=-S,-S+1,...,S-1,S$, with $S=N/2$.\n",
    "\n",
    "The dimension of each of these subspaces is obtained from the\n",
    "combinatory number $\\left(\n",
    "\\begin{array}{c}\n",
    "N \\\\\n",
    "N_{\\uparrow }\n",
    "\\end{array}\n",
    "\\right) = \\frac{N!}{N_{\\uparrow }!N_{\\downarrow }!}$.\n",
    "\n",
    "#### Example: 4-site Heisenberg chain\n",
    "\n",
    "The total basis has\n",
    "$2^{4}=16$ states that can be grouped in\n",
    "$4+1=5$ subspaces with well defined values of the quantum\n",
    "number $m=S^z$. The dimensions of these subspaces are:\n",
    "\n",
    " m | dimension \n",
    ":--:|:--:\n",
    "-2 | 1 \n",
    "-1 | 4 \n",
    "0 | 6 \n",
    "2 | 4 \n",
    "1 | 1 \n",
    " \n",
    "Since we know that the ground state of the Hilbert\n",
    "chain is a singlet, we are only interested in the subspace with\n",
    "</span>$S_{total}^{z}=m=0$<span>. For our example, this subspace is\n",
    "given by (in increasing order): </span> \n",
    "$$\\begin{eqnarray}\n",
    "\\mid \\downarrow \\downarrow \\uparrow \\uparrow \\rangle & \\equiv &|0011\\rangle\n",
    "\\equiv |3\\rangle , \\\\\n",
    "\\mid \\downarrow \\uparrow \\downarrow \\uparrow \\rangle &\\equiv &|0101\\rangle\n",
    "\\equiv |5\\rangle , \\\\\n",
    "\\mid \\downarrow \\uparrow \\uparrow \\downarrow \\rangle &\\equiv &|0110\\rangle\n",
    "\\equiv |6\\rangle , \\\\\n",
    "\\mid \\uparrow \\downarrow \\downarrow \\uparrow \\rangle &\\equiv &|1001\\rangle\n",
    "\\equiv |9\\rangle , \\\\\n",
    "\\mid \\uparrow \\downarrow \\uparrow \\downarrow \\rangle &\\equiv &|1010\\rangle\n",
    "\\equiv |10\\rangle , \\\\\n",
    "\\mid \\uparrow \\uparrow \\downarrow \\downarrow \\rangle &\\equiv &|1010\\rangle\n",
    "\\equiv |12\\rangle .\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Generating the possible configurations of $N$ up-spins (or\n",
    "spinless fermions) in $L$ sites is equivalent to the problem of\n",
    "generating the corresponding combinations of zeroes and ones in\n",
    "lexicographic order.\n",
    "\n",
    "## Construction of the Hamiltonian matrix elements\n",
    "\n",
    "We now have to address the problem of generating all the Hamiltonian matrix elements\n",
    ". These are obtained by the application of the\n",
    "Hamiltonian on each state of the basis $|\\phi \\rangle $, generating all\n",
    "the values\n",
    "$H_{\\phi ,\\phi ^{\\prime }}=\\langle \\phi ^{\\prime }|{{\\hat{H}}}%\n",
    "|\\phi \\rangle $. We illustrate this procedure in the following\n",
    "pseudocode:\n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;**for** each term ${\\hat{H}_i}$ of the Hamiltonian ${\\hat{H}}$\n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;**for** all the states $|\\phi\\rangle $ in the basis\n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;$\\hat{H}_i|\\phi\\rangle = \\langle \\phi^{\\prime }|{\\hat{H}_i}|\\phi\\rangle |\\phi^{\\prime }\\rangle  $ \n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;\n",
    "$H_{\\phi,\\phi^{\\prime }}=H_{\\phi,\\phi^{\\prime }} + \\langle \\phi^{\\prime }|{\\hat{H}_i}|\\phi\\rangle $\n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;&ensp;**end for**\n",
    "\n",
    "&ensp;&ensp;&ensp;&ensp;**end for**\n",
    "\n",
    "As the states are representated by binary words, the spin operators (as\n",
    "well as the creation and destructions operators) can be be easily\n",
    "implemented using the logical operators `AND, OR, XOR`.\n",
    "\n",
    "As a practical example let us consider the spin-1/2 Heisenberg\n",
    "Hamiltonian on the 4-sites chain:\n",
    "$${{\\hat{H}}}=J\\sum_{i=0}^{3}\\mathbf{S}_{i}.\\mathbf{S}_{i+1}=J%\n",
    "\\sum_{i=0}^{3}S_{i}^{z}.S_{i+1}^{z}+\\frac{J}{2}%\n",
    "\\sum_{i=0}^{3}(S_{i}^{+}S_{i+1}^{-}+S_{i}^{-}S_{i+1}^{+}), $$\n",
    "\n",
    "with $\\mathbf{S}_{4}=\\mathbf{S}_{0}$, due to the periodic boundary\n",
    "conditions. To evaluate the matrix elements in the basis (\\[basis0\\]) we\n",
    "must apply each term of ${{\\hat{H}}}$ on such states. The first term is\n",
    "called Ising term, and is diagonal in this basis (also called Ising\n",
    "basis). The last terms, or fluctuation terms, give strictly off-diagonal\n",
    "contributions to the matrix in this representation (when symmetries are\n",
    "considered, they can also have diagonal contributions). These\n",
    "fluctuations cause an exchange in the spin orientation between neighbors\n",
    "with opposite spins, e.g. $\\uparrow \\downarrow \\rightarrow $\n",
    "$\\downarrow \\uparrow $. The way to implement spin-flip operations of\n",
    "this kind on the computer is defining new ‘masks’. A mask for this\n",
    "operation is a binary number with ‘zeroes’ everywhere, except in the\n",
    "positions of the spins to be flipped, e.g 00..0110...00. Then, the\n",
    "logical operator `XOR` is used between the initial state and the mask to\n",
    "invert the bits (spins) at the positions indicated by the mask. For\n",
    "example, (0101)$\\mathtt{.}$`XOR.`(1100)=(1001), i.e., 5`.XOR.`12=9. It\n",
    "is useful to store all the masks for a given geometry in memory,\n",
    "generating them immediately after the tables for the sites and neighbors.\n",
    "\n",
    "<span>In the Heisenberg model the spin flips can be implemented using\n",
    "masks, with ‘zeroes’ everywhere, except in the positions of the spins to\n",
    "be flipped. To illustrate this let us show the effect of the\n",
    "off-diagonal terms on one of the Néel configurations:</span>\n",
    "$$\\begin{eqnarray}\n",
    "S_{0}^{+}S_{1}^{-}\\,+S_{0}^{-}S_{1}^{+}\\,|0101\\rangle &\\equiv &(0011)\\mathtt{%\n",
    ".XOR.}(0101)=3\\mathtt{.XOR.}5=(0110)=6\\equiv \\,|0110\\rangle , \\\\\n",
    "S_{1}^{+}S_{2}^{-}+S_{1}^{-}S_{2}^{+}\\,\\,|0101\\rangle &\\equiv &(0110)\\mathtt{%\n",
    ".XOR.}(0101)=6\\mathtt{.XOR.}5=(0011)=3\\equiv \\,|0011\\rangle , \\\\\n",
    "S_{2}^{+}S_{3}^{-}+S_{2}^{-}S_{3}^{+}\\,\\,|0101\\rangle &\\equiv &(1100)\\mathtt{%\n",
    ".XOR.}(0101)=12\\mathtt{.XOR.}5=(1001)=9\\equiv \\,|1001\\rangle , \\\\\n",
    "S_{3}^{+}S_{0}^{-}\\,+S_{3}^{-}S_{0}^{+}|0101\\rangle &\\equiv &(1001)\\mathtt{%\n",
    ".XOR.}(0101)=9\\mathtt{.XOR.}5=(1100)=12\\equiv \\,|1100\\rangle .%\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "<span>After applying the Hamiltonian (\\[h4\\]) on our basis states, the\n",
    "reader can verify as an exercise that we obtain </span>\n",
    "$$\\begin{eqnarray}\n",
    "{{\\hat{H}}}\\,|0101\\rangle &=&-J\\,|0101\\rangle +\\frac{J}{2}\\left[\n",
    "\\,|1100\\rangle +\\,|1001\\rangle +\\,|0011\\rangle +\\,|0110\\rangle \\right] , \\\\\n",
    "{{\\hat{H}}}\\,|1010\\rangle &=&-J\\,|1010\\rangle +\\frac{J}{2}\\left[\n",
    "\\,|1100\\rangle +\\,|1001\\rangle +\\,|0011\\rangle +\\,|0110\\rangle \\right] , \\\\\n",
    "{{\\hat{H}}}\\,|0011\\rangle &=&\\frac{J}{2}\\left[ \\,\\,|0101\\rangle\n",
    "+\\,|1010\\rangle \\right] , \\\\\n",
    "{{\\hat{H}}}\\,|0110\\rangle &=&\\frac{J}{2}\\left[ \\,\\,|0101\\rangle\n",
    "+\\,|1010\\rangle \\right] , \\\\\n",
    "{{\\hat{H}}}\\,|1001\\rangle &=&\\frac{J}{2}\\left[ \\,|0101\\rangle +\\,|1010\\rangle\n",
    "\\right] , \\\\\n",
    "{{\\hat{H}}}\\,|1100\\rangle &=&\\frac{J}{2}\\left[ \\,|0101\\rangle +\\,|1010\\rangle\n",
    "\\right] .%\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "<span>The resulting matrix for the operator </span>$\\hat{H}$<span>in the\n",
    "</span>$S^{z}$<span> representation is: </span> $$H=J\\left(\n",
    "\\begin{array}{llllll}\n",
    "0 & 1/2 & 0 & 0 & 1/2 & 0 \\\\\n",
    "1/2 & -1 & 1/2 & 1/2 & 0 & 1/2 \\\\\n",
    "0 & 1/2 & 0 & 0 & 1/2 & 0 \\\\\n",
    "0 & 1/2 & 0 & 0 & 1/2 & 0 \\\\\n",
    "1/2 & 0 & 1/2 & 1/2 & -1 & 1/2 \\\\\n",
    "0 & 1/2 & 0 & 0 & 1/2 & 0\n",
    "\\end{array}\n",
    "\\right) .  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Basis:\n",
      "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255]\n",
      "[[2.  0.  0.  ... 0.  0.  0. ]\n",
      " [0.  1.  0.5 ... 0.  0.  0. ]\n",
      " [0.  0.5 1.  ... 0.  0.  0. ]\n",
      " ...\n",
      " [0.  0.  0.  ... 1.  0.5 0. ]\n",
      " [0.  0.  0.  ... 0.5 1.  0. ]\n",
      " [0.  0.  0.  ... 0.  0.  2. ]]\n",
      "===================================================================================================================\n",
      "[-3.65109341e+00 -3.12841906e+00 -3.12841906e+00 -3.12841906e+00\n",
      " -2.69962815e+00 -2.45873851e+00 -2.45873851e+00 -2.45873851e+00\n",
      " -2.45873851e+00 -2.45873851e+00 -2.45873851e+00 -2.14514837e+00\n",
      " -2.14514837e+00 -2.14514837e+00 -2.14514837e+00 -2.14514837e+00\n",
      " -2.14514837e+00 -1.85463768e+00 -1.85463768e+00 -1.85463768e+00\n",
      " -1.85463768e+00 -1.85463768e+00 -1.85463768e+00 -1.80193774e+00\n",
      " -1.80193774e+00 -1.80193774e+00 -1.80193774e+00 -1.80193774e+00\n",
      " -1.70710678e+00 -1.70710678e+00 -1.61803399e+00 -1.61803399e+00\n",
      " -1.26703510e+00 -1.26703510e+00 -1.26703510e+00 -1.26703510e+00\n",
      " -1.26703510e+00 -1.26703510e+00 -1.26703510e+00 -1.26703510e+00\n",
      " -1.26703510e+00 -1.26703510e+00 -1.20163968e+00 -1.20163968e+00\n",
      " -1.20163968e+00 -1.14412281e+00 -1.14412281e+00 -1.14412281e+00\n",
      " -1.14412281e+00 -1.14412281e+00 -1.14412281e+00 -1.14412281e+00\n",
      " -1.14412281e+00 -1.14412281e+00 -1.14412281e+00 -1.00000000e+00\n",
      " -1.00000000e+00 -1.00000000e+00 -8.58923550e-01 -8.58923550e-01\n",
      " -8.58923550e-01 -8.58923550e-01 -8.58923550e-01 -8.58923550e-01\n",
      " -7.60876722e-01 -7.26109445e-01 -6.28378427e-01 -6.28378427e-01\n",
      " -6.28378427e-01 -6.28378427e-01 -6.28378427e-01 -6.28378427e-01\n",
      " -5.96968283e-01 -5.96968283e-01 -5.96968283e-01 -5.96968283e-01\n",
      " -5.96968283e-01 -5.96968283e-01 -5.12683203e-01 -5.12683203e-01\n",
      " -5.12683203e-01 -5.12683203e-01 -5.12683203e-01 -5.12683203e-01\n",
      " -4.45041868e-01 -4.45041868e-01 -4.45041868e-01 -4.45041868e-01\n",
      " -4.45041868e-01 -4.37016024e-01 -4.37016024e-01 -4.37016024e-01\n",
      " -4.37016024e-01 -4.37016024e-01 -4.37016024e-01 -4.37016024e-01\n",
      " -4.37016024e-01 -4.37016024e-01 -4.37016024e-01 -2.92893219e-01\n",
      " -2.92893219e-01 -2.58652023e-01 -2.58652023e-01 -2.58652023e-01\n",
      " -2.58652023e-01 -2.58652023e-01 -2.58652023e-01 -2.58652023e-01\n",
      " -2.58652023e-01 -2.58652023e-01 -2.58652023e-01 -7.75689638e-16\n",
      " -7.63208968e-16 -4.57134676e-16 -4.13645713e-16 -3.39947086e-16\n",
      " -3.06926362e-16 -2.34643233e-16 -1.74475331e-16 -1.72272455e-16\n",
      " -1.38134367e-16 -1.37645901e-16 -1.32021661e-16 -7.69140168e-17\n",
      "  4.75886108e-17  5.83588051e-17  7.00756471e-17  9.12342182e-17\n",
      "  1.26508026e-16  2.24592706e-16  2.86790224e-16  4.82216812e-16\n",
      "  5.84928975e-16  9.02098163e-16  2.65452712e-01  2.65452712e-01\n",
      "  2.65452712e-01  2.65452712e-01  2.65452712e-01  2.65452712e-01\n",
      "  2.92893219e-01  2.92893219e-01  2.92893219e-01  2.92893219e-01\n",
      "  2.92893219e-01  2.92893219e-01  2.92893219e-01  2.92893219e-01\n",
      "  2.92893219e-01  2.92893219e-01  2.92893219e-01  2.92893219e-01\n",
      "  2.92893219e-01  2.92893219e-01  3.77202854e-01  4.37016024e-01\n",
      "  4.37016024e-01  4.37016024e-01  4.37016024e-01  4.37016024e-01\n",
      "  4.37016024e-01  4.37016024e-01  4.37016024e-01  4.37016024e-01\n",
      "  4.37016024e-01  4.51605963e-01  4.51605963e-01  4.51605963e-01\n",
      "  4.51605963e-01  4.51605963e-01  4.51605963e-01  6.18033989e-01\n",
      "  6.18033989e-01  8.92693358e-01  8.92693358e-01  8.92693358e-01\n",
      "  8.92693358e-01  8.92693358e-01  8.92693358e-01  1.00000000e+00\n",
      "  1.00000000e+00  1.00000000e+00  1.00000000e+00  1.00000000e+00\n",
      "  1.00000000e+00  1.00000000e+00  1.00000000e+00  1.00000000e+00\n",
      "  1.00000000e+00  1.00000000e+00  1.00000000e+00  1.00000000e+00\n",
      "  1.00000000e+00  1.00000000e+00  1.00000000e+00  1.00000000e+00\n",
      "  1.00000000e+00  1.00000000e+00  1.14412281e+00  1.14412281e+00\n",
      "  1.14412281e+00  1.14412281e+00  1.14412281e+00  1.14412281e+00\n",
      "  1.14412281e+00  1.14412281e+00  1.14412281e+00  1.14412281e+00\n",
      "  1.24697960e+00  1.24697960e+00  1.24697960e+00  1.24697960e+00\n",
      "  1.24697960e+00  1.33005874e+00  1.33005874e+00  1.33005874e+00\n",
      "  1.44572599e+00  1.44572599e+00  1.44572599e+00  1.44572599e+00\n",
      "  1.44572599e+00  1.44572599e+00  1.46050487e+00  1.52568712e+00\n",
      "  1.52568712e+00  1.52568712e+00  1.52568712e+00  1.52568712e+00\n",
      "  1.52568712e+00  1.52568712e+00  1.52568712e+00  1.52568712e+00\n",
      "  1.52568712e+00  1.70710678e+00  1.70710678e+00  1.70710678e+00\n",
      "  1.70710678e+00  1.70710678e+00  1.70710678e+00  1.70710678e+00\n",
      "  1.70710678e+00  1.70710678e+00  1.70710678e+00  1.70710678e+00\n",
      "  1.70710678e+00  1.70710678e+00  1.70710678e+00  2.00000000e+00\n",
      "  2.00000000e+00  2.00000000e+00  2.00000000e+00  2.00000000e+00\n",
      "  2.00000000e+00  2.00000000e+00  2.00000000e+00  2.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "class BoundaryCondition:\n",
    "    RBC, PBC = range(2)\n",
    "\n",
    "class Direction:\n",
    "    RIGHT, TOP, LEFT, BOTTOM = range(4)\n",
    "    \n",
    "L = 8\n",
    "Nup = 2\n",
    "maxdim = 2**L\n",
    "bc = BoundaryCondition.PBC\n",
    "\n",
    "#  Tables for energy\n",
    "\n",
    "hdiag = np.zeros(4) # Diagonal energies\n",
    "hflip = np.zeros(4) # off-diagonal terms\n",
    "\n",
    "hdiag[0] = +0.25\n",
    "hdiag[1] = -0.25\n",
    "hdiag[2] = -0.25\n",
    "hdiag[3] = +0.25\n",
    "\n",
    "hflip[0] = 0.\n",
    "hflip[1] = 0.5\n",
    "hflip[2] = 0.5\n",
    "hflip[3] = 0.\n",
    "\n",
    "#hflip *= 2\n",
    "#hdiag *= 0.\n",
    "\n",
    "# Lattice geometry  (1D chain)\n",
    "\n",
    "nn = np.zeros(shape=(L,4), dtype=np.int16)\n",
    "for i in range(L):\n",
    "    nn[i, Direction.RIGHT] = i-1\n",
    "    nn[i, Direction.LEFT] = i+1\n",
    "\n",
    "if(bc == BoundaryCondition.RBC):   # Open Boundary Conditions\n",
    "    nn[0, Direction.RIGHT] = -1    # This means error\n",
    "    nn[L-1, Direction.LEFT] = -1\n",
    "else:                              # Periodic Boundary Conditions\n",
    "    nn[0, Direction.RIGHT] = L-1   # We close the ring\n",
    "    nn[L-1, Direction.LEFT] = 0\n",
    "  \n",
    "# We build basis\n",
    "\n",
    "basis = []\n",
    "\n",
    "dim = 0\n",
    "for state in range(maxdim):\n",
    "    basis.append(state)\n",
    "    dim += 1\n",
    "\n",
    "print (\"Basis:\")\n",
    "print (basis)\n",
    "\n",
    "# We build Hamiltonian matrix\n",
    "\n",
    "H = np.zeros(shape=(dim,dim))\n",
    "\n",
    "def IBITS(n,i):\n",
    "    return ((n >> i) & 1)\n",
    "\n",
    "for i in range(dim):\n",
    "    state = basis[i]\n",
    "\n",
    "    # Diagonal term\n",
    "    for site_i in range(L):\n",
    "        site_j = nn[site_i, Direction.RIGHT]\n",
    "        if(site_j != -1):   # This would happen for open boundary conditions\n",
    "            two_sites  = IBITS(state,site_i) | (IBITS(state,site_j) << 1)\n",
    "            value = hdiag[two_sites]\n",
    "            H[i,i] += value\n",
    "\n",
    "        \n",
    "for i in range(dim):\n",
    "    state = basis[i]\n",
    "\n",
    "    # Off-diagonal term\n",
    "    for site_i in range(L):\n",
    "        site_j = nn[site_i, Direction.RIGHT]\n",
    "\n",
    "        if(site_j != -1):\n",
    "            mask = (1 << site_i) | (1 << site_j)\n",
    "            two_sites  = IBITS(state,site_i) | (IBITS(state,site_j) << 1)\n",
    "            value = hflip[two_sites]\n",
    "            if(value != 0.):\n",
    "                new_state = (state ^ mask)\n",
    "                j = new_state\n",
    "                H[i,j] += value\n",
    "\n",
    "print(H)\n",
    "d, v = np.linalg.eigh(H) #Diagonalize the matrix \n",
    "print(\"===================================================================================================================\")\n",
    "print(d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0\n",
      "1 1\n",
      "2 1\n",
      "3 2\n",
      "4 1\n",
      "5 2\n",
      "6 2\n",
      "7 3\n",
      "8 1\n",
      "9 2\n",
      "10 2\n",
      "11 3\n",
      "12 2\n",
      "13 3\n",
      "14 3\n",
      "15 4\n",
      "Dim= 6\n",
      "Basis:\n",
      "[3, 5, 6, 9, 10, 12]\n",
      "[[ 0.25  0.5   0.    0.    0.    0.  ]\n",
      " [ 0.5  -0.75  0.5   0.5   0.    0.  ]\n",
      " [ 0.    0.5  -0.25  0.    0.5   0.  ]\n",
      " [ 0.    0.5   0.   -0.25  0.5   0.  ]\n",
      " [ 0.    0.    0.5   0.5  -0.75  0.5 ]\n",
      " [ 0.    0.    0.    0.    0.5   0.25]]\n",
      "[-1.6160254  -0.95710678 -0.25        0.1160254   0.45710678  0.75      ] -1.6160254037844386\n",
      "[ 0.14942925 -0.55767754  0.40824829  0.40824829 -0.55767754  0.14942925]\n"
     ]
    }
   ],
   "source": [
    "# Using quantum number conservation : Sz\n",
    "\n",
    "basis = []\n",
    "\n",
    "dim = 0\n",
    "for state in range(maxdim):\n",
    "    n_ones = 0\n",
    "    for bit in range(L):\n",
    "        n_ones += IBITS(state,bit)\n",
    "    print (state,n_ones)\n",
    "    if(n_ones == Nup):\n",
    "        basis.append(state)\n",
    "        dim += 1\n",
    "\n",
    "print (\"Dim=\",dim)\n",
    "print (\"Basis:\")\n",
    "print (basis)\n",
    "\n",
    "#hflip[:] *= 2\n",
    "# We build Hamiltonian matrix\n",
    "\n",
    "H = np.zeros(shape=(dim,dim))\n",
    "\n",
    "def IBITS(n,i):\n",
    "    return ((n >> i) & 1)\n",
    "\n",
    "for i in range(dim):\n",
    "    state = basis[i]\n",
    "\n",
    "    # Diagonal term\n",
    "    for site_i in range(L):\n",
    "        site_j = nn[site_i, Direction.RIGHT]\n",
    "        if(site_j != -1):   # This would happen for open boundary conditions\n",
    "            two_sites  = IBITS(state,site_i) | (IBITS(state,site_j) << 1)\n",
    "            value = hdiag[two_sites]\n",
    "            H[i,i] += value\n",
    "\n",
    "\n",
    "def bisect(state, basis):  \n",
    "    # Binary search, only works in a sorted list of integers\n",
    "    # state : integer we seek\n",
    "    # basis : list of sorted integers in increasing order\n",
    "    # ret_val : return value, position on the list, -1 if not found\n",
    "    \n",
    "    ret_val = -1\n",
    "    dim = len(basis)\n",
    "    \n",
    "#    for i in range(dim):\n",
    "#        if(state == basis[i]):\n",
    "#            return i\n",
    "\n",
    "    origin = 0\n",
    "    end = dim-1\n",
    "    middle = (origin+end)/2\n",
    "\n",
    "    while(1):\n",
    "        index_old = middle\n",
    "\n",
    "        middle = (origin+end)//2\n",
    "\n",
    "        if(state < basis[middle]):\n",
    "            end = middle\n",
    "        else:\n",
    "            origin = middle\n",
    "\n",
    "        if(basis[middle] == state):\n",
    "            break\n",
    "\n",
    "        if(middle == index_old):\n",
    "            if(middle == end):\n",
    "                end = end - 1\n",
    "            else:\n",
    "                origin = origin + 1\n",
    "\n",
    "    ret_val = middle\n",
    "    return ret_val           \n",
    "\n",
    "            \n",
    "for i in range(dim):\n",
    "    state = basis[i]\n",
    "\n",
    "    # Off-diagonal term\n",
    "    for site_i in range(L):\n",
    "        site_j = nn[site_i, Direction.RIGHT]\n",
    "\n",
    "        if(site_j != -1):\n",
    "            mask = (1 << site_i) | (1 << site_j)\n",
    "            two_sites  = IBITS(state,site_i) | (IBITS(state,site_j) << 1)\n",
    "            value = hflip[two_sites]\n",
    "            if(value != 0.):\n",
    "                new_state = (state ^ mask)\n",
    "                j = bisect(new_state, basis)\n",
    "                H[i,j] += value\n",
    "\n",
    "print(H)\n",
    "d, v = np.linalg.eigh(H) #Diagonalize the matrix \n",
    "print(d,np.min(d))\n",
    "print(v[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtaining the ground-state: Lanczos diagonalization\n",
    "---------------------------------------------------\n",
    "\n",
    "Once we have a superblock matrix, we can apply a library routine to\n",
    "obtain the ground state of the superblock $|\\Psi\\rangle$. The two\n",
    "algorithms widely used for this purpose are the Lanczos and Davidson\n",
    "diagonalization. Both are explained to great extent in Ref.@noack, so we\n",
    "refer the reader to this material for further information. In these\n",
    "notes we will briefly explain the Lanczos procedure.\n",
    "\n",
    "The basic idea of the Lanczos method is that a special basis can be constructed\n",
    "where the Hamiltonian has a tridiagonal representation. This is carried\n",
    "out iteratively as shown below. First, it is necessary to select an\n",
    "arbitrary seed vector $|\\phi _{0}\\rangle $ in the Hilbert space of the\n",
    "model being studied. If we are seeking the ground-state of the model,\n",
    "then it is necessary that the overlap between the actual ground-state\n",
    "$|\\psi\n",
    "_{0}\\rangle $, and the initial state $|\\phi _{0}\\rangle $ be nonzero. If\n",
    "no “a priori” information about the ground state is known, this\n",
    "requirement is usually easily satisfied by selecting an initial state\n",
    "with *randomly* chosen coefficients in the working basis that is being\n",
    "used. If some other information of the ground state is known, like its\n",
    "total momentum and spin, then it is convenient to initiate the\n",
    "iterations with a state already belonging to the subspace having those\n",
    "quantum numbers (and still with random coefficients within this\n",
    "subspace).\n",
    "\n",
    "After $|\\phi _{0}\\rangle $ is selected, define a new vector by applying\n",
    "the Hamiltonian ${{\\hat{H}}}$, over the initial state. Subtracting the\n",
    "projection over $|\\phi _{0}\\rangle $, we obtain\n",
    "$$|\\phi _{1}\\rangle ={{\\hat{H}}}|\\phi _{0}\\rangle -\\frac{{\\langle }\\phi _{0}|{{%\n",
    "\\hat{H}}}|{\\phi }_{0}{\\rangle }}{\\langle \\phi _{0}|\\phi _{0}\\rangle }|\\phi\n",
    "_{0}\\rangle ,$$ that satisfies $\\langle \\phi _{0}|\\phi _{1}\\rangle =0$.\n",
    "Now, we can construct a new state that is orthogonal to the previous two\n",
    "as,\n",
    "$$|\\phi _{2}\\rangle ={{\\hat{H}}}|\\phi _{1}\\rangle -{\\frac{{\\langle \\phi _{1}|{%\n",
    "\\hat{H}}|\\phi _{1}\\rangle }}{{\\langle \\phi _{1}|\\phi _{1}\\rangle }}}|\\phi\n",
    "_{1}\\rangle -{\\frac{{\\langle \\phi _{1}|\\phi _{1}\\rangle }}{{\\langle \\phi\n",
    "_{0}|\\phi _{0}\\rangle }}}|\\phi _{0}\\rangle .$$ It can be easily checked\n",
    "that $\\langle \\phi _{0}|\\phi _{2}\\rangle\n",
    "=\\langle \\phi _{1}|\\phi _{2}\\rangle =0$. The procedure can be\n",
    "generalized by defining an orthogonal basis recursively as,\n",
    "$$|\\phi _{n+1}\\rangle ={{\\hat{H}}}|\\phi _{n}\\rangle -a_{n}|\\phi _{n}\\rangle\n",
    "-b_{n}^{2}|\\phi _{n-1}\\rangle ,$$ where $n=0,1,2,...$, and the\n",
    "coefficients are given by\n",
    "$$a_{n}={\\frac{{\\langle \\phi _{n}|{\\hat{H}}|\\phi _{n}\\rangle }}{{\\langle \\phi\n",
    "_{n}|\\phi _{n}\\rangle }}},\\qquad b_{n}^{2}={\\frac{{\\langle \\phi _{n}|\\phi\n",
    "_{n}\\rangle }}{{\\langle \\phi _{n-1}|\\phi _{n-1}\\rangle }}},$$\n",
    "supplemented by $b_{0}=0$, $|\\phi _{-1}\\rangle =0$. In this basis, it\n",
    "can be shown that the Hamiltonian matrix becomes, $$H=\\left(\n",
    "\\begin{array}{lllll}\n",
    "a_{0} & b_{1} & 0 & 0 & ... \\\\\n",
    "b_{1} & a_{1} & b_{2} & 0 & ... \\\\\n",
    "0 & b_{2} & a_{2} & b_{3} & ... \\\\\n",
    "0 & 0 & b_{3} & a_{3} & ... \\\\\n",
    "\\vdots\n",
    "{} &\n",
    "\\vdots\n",
    "&\n",
    "\\vdots\n",
    "&\n",
    "\\vdots\n",
    "&\n",
    "\\end{array}\n",
    "\\right)$$ i.e. it is tridiagonal as expected. Once in this form the\n",
    "matrix can be diagonalized easily using standard library subroutines.\n",
    "However, note that to diagonalize completely a Hamiltonian on a finite\n",
    "cluster, a number of iterations equal to the size of the Hilbert space\n",
    "(or the subspace under consideration) are needed. In practice this would\n",
    "demand a considerable amount of CPU time. However, one of the advantages\n",
    "of this technique is that accurate enough information about the ground\n",
    "state of the problem can be obtained after a small number of iterations\n",
    "(typically of the order of $\\sim 100$ or less).\n",
    "\n",
    "Another way to formulate the problem is by obtaining the tridiagonal\n",
    "form of the Hamiltonian starting from a Krylov basis, which is spanned\n",
    "by the vectors\n",
    "$$\\left\\{|\\phi_0\\rangle,\\hat{H}|\\phi_0\\rangle,\\hat{H}^2|\\phi_0\\rangle,...,\\hat{H}^n|\\phi_0\\rangle\\right\\}$$\n",
    "and asking that each vector be orthogonal to the previous two. Notice\n",
    "that each new iteration of the process requires one application of the\n",
    "Hamiltonian. Most of the time this simple procedure works for practical\n",
    "purposes, but care must be payed to the possibility of losing\n",
    "orthogonality between the basis vectors. This may happen due to the\n",
    "finite machine precision. In that case, a re-orthogonalization procedure\n",
    "may be required.\n",
    "\n",
    "Notice that the new super-Hamiltonian matrix has dimensions\n",
    "$D_L D_R d^2 \\times D_L D_R d^2$. This could be a large matrix. In\n",
    "state-of-the-art simulations with a large number of states, one does not\n",
    "build this matrix in memory explicitly, but applies the operators to the\n",
    "state directly in the diagonalization routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lanczos(m, seed, maxiter, tol, use_seed, calc_gs, force_maxiter = False):\n",
    "    x1 = seed\n",
    "    x2 = seed\n",
    "    gs = seed\n",
    "    a = np.zeros(100)\n",
    "    b = np.zeros(100)\n",
    "    z = np.zeros((100,100))\n",
    "    lvectors = []\n",
    "    control_max = maxiter;\n",
    "\n",
    "    if(maxiter == -1):\n",
    "        force_maxiter = False\n",
    "\n",
    "    if(control_max == 0):\n",
    "        gs = 1\n",
    "        maxiter = 1\n",
    "        return(e0,gs)\n",
    "    \n",
    "    x1[:] = 0\n",
    "    x2[:] = 0\n",
    "    gs[:] = 0\n",
    "    maxiter = 0\n",
    "    a[:] = 0.0\n",
    "    b[:] = 0.0\n",
    "    if(use_seed):\n",
    "        x1 = seed\n",
    "    else :\n",
    "        x1 = np.random.random(x1.shape[0])*2-1.\n",
    "\n",
    "    b[0] = np.sqrt(np.dot(x1,x1))\n",
    "    x1 = x1 / b[0]\n",
    "    x2[:] = 0\n",
    "    b[0] = 1.\n",
    "\n",
    "    e0 = 9999\n",
    "    nmax = min(99, gs.shape[0])\n",
    "    for iter in range(1,nmax+1):\n",
    "        eini = e0\n",
    "        if(b[iter - 1] != 0.):\n",
    "            aux = x1\n",
    "            x1 = -b[iter-1] * x2\n",
    "            x2 = aux / b[iter-1]\n",
    "\n",
    "        aux = np.dot(m,x2)\n",
    "        x1 = x1 + aux\n",
    "        a[iter] = np.dot(x1, x2)\n",
    "        x1 = x1 - x2*a[iter]\n",
    "        b[iter] = np.sqrt(np.dot(x1, x1))\n",
    "        lvectors.append(x2)                                          \n",
    "\n",
    "        z.resize((iter+1,iter+1))\n",
    "        z[:,:] = 0\n",
    "        for i in range(0,iter-1):\n",
    "            z[i,i+1] = b[i+1]\n",
    "            z[i+1,i] = b[i+1]\n",
    "            z[i,i] = a[i+1]\n",
    "        z[iter-1,iter-1]=a[iter]\n",
    "        d, v = np.linalg.eig(z)\n",
    "        col = 0\n",
    "        n = 0\n",
    "        e0 = 9999\n",
    "        for e in d:\n",
    "            if(e < e0):\n",
    "                e0 = e\n",
    "                col = n\n",
    "            n+=1\n",
    "        e0 = d[col]\n",
    "        print (\"Iter = \",iter,\" Ener = \",e0)\n",
    "        if((force_maxiter and iter >= control_max) or (iter >= gs.shape[0] or iter == 99 or abs(b[iter]) < tol) or \\\n",
    "            ((not force_maxiter) and abs(eini-e0) <= tol)):\n",
    "            # converged\n",
    "            gs = 0.\n",
    "            for i in range(0,iter):\n",
    "                gs += v[i,col]*lvectors[i]\n",
    "            print (\"E0 = \", e0, np.sqrt(np.dot(gs,gs)))\n",
    "            maxiter = iter\n",
    "            return(e0,gs) # We return with ground states energy\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = np.zeros(H.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter =  1  Ener =  -0.6848415206892902\n",
      "Iter =  2  Ener =  -1.9335584929208487\n",
      "Iter =  3  Ener =  -1.9939600915319768\n",
      "Iter =  4  Ener =  -1.999999999999999\n",
      "E0 =  -1.999999999999999 0.9999999999999997\n"
     ]
    }
   ],
   "source": [
    "e0, gs = lanczos(H,seed,6,1.e-5,False,False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.28867513 -0.57735027  0.28867513  0.28867513 -0.57735027  0.28867513]\n"
     ]
    }
   ],
   "source": [
    "print(gs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999999999994\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(gs,gs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}