Therm limit, a la MZ & OG

gistfile1.txt

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "import mcba.db as db\n",
      "from mcba.models.impurity import gamma, initial_q, k_F, SingleImpurity, fsPairs\n",
      "from mcba.models.impurity.BA import buckets_from_pairs as buckets_from_pairs, is_valid as is_valid\n",
      "from mcba.models.impurity.BA import Lieb_Wu_a\n",
      "import mcba.models.impurity.matrix_elements as mx\n",
      "\n",
      "import postproc as pp\n",
      "\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fname = \"/home/br/sweethome/1Dimpur/data/q25/all/N45g6mq27.sqlite\"\n",
      "handle, = db.get_handles(fname)\n",
      "par = db.get_param(handle)\n",
      "\n",
      "print par, initial_q(par)/k_F(par)\n",
      "\n",
      "cs = [row[\"c\"] for row in db.row_iterator(handle) if np.abs(row[\"c\"]) != np.inf]\n",
      "cs = np.asarray(cs)\n",
      "print len(cs), min(cs), max(cs)\n",
      "\n",
      "n, bin, patches = plt.hist(cs[np.abs(cs) < 5], 50, normed=1, facecolor='green', alpha=0.75)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Par(N=45, L=135, V=1.0, m_q=27) 1.2\n",
        "89562"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " -69.6944864348 109.755550406\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFcVJREFUeJzt3X9sU9fdx/GPeZKqatlDAQVC4mhefqwxa5NlSshohcgk\nQmq2RiKqpqioq0bGLKSs6qZuiH9GQOqPbEITI9qW8UuqWkX5b0EbcxFa3XYRwd0GYlJKFdaEOm7T\nsJWUAltDzHn+6FOvJo7jEDvXOX6/JEu+N+eY75Xaj6/PPedelzHGCABgnSVOFwAAyAwCHgAsRcAD\ngKUIeACwFAEPAJYi4AHAUrMGfCAQUGVlpSoqKtTZ2Tnt7y+//LKqq6tVVVWlhx9+WOfPn4/9zePx\nqKqqSjU1NVq3bl16KwcAJOVKNg8+Go3q/vvv16lTp1RcXKy6ujr19PTI6/XG2pw+fVpr167VsmXL\nFAgE1NHRoYGBAUnSl770Jf31r3/VihUrMn8kAIA4Sc/gQ6GQysvL5fF4lJ+fr9bWVvX19cW1Wb9+\nvZYtWyZJqq+v1+joaNzfWUcFAM5IGvCRSEQlJSWxbbfbrUgkMmP7I0eOaMuWLbFtl8ulTZs2qba2\nVocOHUpDuQCAVOUl+6PL5Ur5g1599VUdPXpU/f39sX39/f1as2aNLl++rMbGRlVWVmrDhg13Xi0A\nIGVJA764uFjhcDi2HQ6H5Xa7p7U7f/68duzYoUAgoOXLl8f2r1mzRpJUUFCgrVu3KhQKTQv48vJy\n/eMf/5jXQQBArikrK9PFixeTNzJJ3Lx505SWlprh4WHzySefmOrqajM4OBjX5tKlS6asrMycPn06\nbv/169fN1atXjTHGXLt2zTz00EPmlVdemfZvzFLCordnzx6nS8gojm/xsvnYjLH/+FLJzqRn8Hl5\neerq6lJTU5Oi0aja2trk9XrV3d0tSfL7/dq3b5+uXLminTt3SpLy8/MVCoU0NjamlpYWSdLU1JS2\nbdumzZs3z/tbCwCQmqQBL0k+n08+ny9un9/vj70/fPiwDh8+PK1faWmpzp07l4YSAQB3gpWsGdbQ\n0OB0CRnF8S1eNh+bZP/xpSLpQqcFKcDlYq48AMxRKtnJGTwAWIqABwBLEfAAYCkCHgAsRcADgKUI\neACwFAEPAJYi4AHAUgQ8AFiKgAcASxHwAGApAh4ALEXAA4ClCHgAsBQBDwCWmvWJTkCmbWzaqPEP\nx6ftX7VilV575TUHKgLsQMDDceMfjquwvXDa/rGuMQeqAezBEA0AWIqABwBLEfAAYCkCHgAsRcAD\ngKUIeACwFNMkkXMSzbtnzj1sRMAj5ySad8+ce9iIIRoAsBQBDwCWIuABwFIEPABYioAHAEsR8ABg\nKaZJAkkwZx6LGQEPJMGceSxmDNEAgKVmDfhAIKDKykpVVFSos7Nz2t9ffvllVVdXq6qqSg8//LDO\nnz+fcl8gXTY2bZS3zhv32ti00emyAEclHaKJRqNqb2/XqVOnVFxcrLq6OjU3N8vr9cbalJaW6vXX\nX9eyZcsUCAT0/e9/XwMDAyn1BdKFoRRguqQBHwqFVF5eLo/HI0lqbW1VX19fXEivX78+9r6+vl6j\no6Mp9wVswgVZZJukAR+JRFRSUhLbdrvdOnPmzIztjxw5oi1bttxRX2Cx41cEsk3SgHe5XCl/0Kuv\nvqqjR4+qv79/zn07Ojpi7xsaGtTQ0JByXwDIBcFgUMFgcE59kgZ8cXGxwuFwbDscDsvtdk9rd/78\nee3YsUOBQEDLly+fU18pPuABANPdfvK7d+/eWfsknUVTW1uroaEhjYyMaHJyUr29vWpubo5r8+67\n76qlpUUvvfSSysvL59QXAJA5Sc/g8/Ly1NXVpaamJkWjUbW1tcnr9aq7u1uS5Pf7tW/fPl25ckU7\nd+6UJOXn5ysUCs3YFwCwMGZdyerz+eTz+eL2+f3+2PvDhw/r8OHDKfcFwIwbLAxuVQA4gBk3WAjc\nqgAALEXAA4ClCHgAsBQBDwCWIuABwFLMosGikmh6oSRdCl9SoQoT9AByFwGPRSXR9EJJeueZdxyo\nBshuDNEAgKUIeACwFAEPAJYi4AHAUlxkBRahmWYTccMyfB4BDyxCM80m4oZl+DyGaADAUgQ8AFiK\ngAcASxHwAGApLrJiQSWa/cF9ZIDMIOCxoBLN/uA+MkBmMEQDAJYi4AHAUgQ8AFiKgAcASxHwAGAp\nAh4ALEXAA4ClCHgAsBQBDwCWIuABwFIEPABYioAHAEsR8ABgKQIeACxFwAOApWYN+EAgoMrKSlVU\nVKizs3Pa3y9cuKD169fr7rvv1v79++P+5vF4VFVVpZqaGq1bty59VSOrbGzaKG+dN+61sWmj02UB\nOS/pAz+i0aja29t16tQpFRcXq66uTs3NzfJ6vbE2K1eu1MGDB/W73/1uWn+Xy6VgMKgVK1akv3Is\nuERPY5I+fSJTfWd93L6xrrGFKgvADJIGfCgUUnl5uTwejySptbVVfX19cQFfUFCggoIC/eEPf0j4\nGcaY9FULRyV6GpPEE5mAbJV0iCYSiaikpCS27Xa7FYlEUv5wl8ulTZs2qba2VocOHbrzKgEAc5b0\nDN7lcs3rw/v7+7VmzRpdvnxZjY2Nqqys1IYNG+b1mcgdl0YuyVvnjd/HA7qBlCUN+OLiYoXD4dh2\nOByW2+1O+cPXrFkj6dNhnK1btyoUCiUM+I6Ojtj7hoYGNTQ0pPxvwF5RRXlAN/D/gsGggsHgnPok\nDfja2loNDQ1pZGRERUVF6u3tVU9PT8K2t4+137hxQ9FoVF/4whd0/fp1nTx5Unv27EnY9/MBDwCY\n7vaT3717987aJ2nA5+XlqaurS01NTYpGo2pra5PX61V3d7ckye/3a2xsTHV1dbp69aqWLFmiAwcO\naHBwUOPj42ppaZEkTU1Nadu2bdq8efM8Dg8AMBdJA16SfD6ffD5f3D6/3x97X1hYGDeM85mlS5fq\n3LlzaSgRAHAnZg144E4kukAqcZEUWEgEPDIi0QVSaWEvkvIlg1xHwCOhRKtWF1swZsOXDOAkAh4J\nJVq1SjACiwt3kwQASxHwAGApAh4ALEXAA4ClCHgAsBSzaACLJJr7v2rFKr32ymsOVQQnEfCARRLN\n/efpWrmLgAfEqlfYiYAHxKpX2ImAB+aIs30sFgQ8MEec7WOxYJokAFiKgAcASxHwAGApxuCBLDHT\nxVsWKuFOEfBAlpjp4m22LlRK9FAYvoyyCwEP4I4keihMtn4Z5SrG4AHAUpzB55BEP6klflYDtiLg\nc0iin9QSP6uReYzXO4OABzKI2xp8ivF6ZxDwQAal47YGib4kcu0LAneGgAeyXKIvCe57g1QwiwYA\nLEXAA4ClGKIBLMctEHIXAQ9YbqYLvWeeOcMDui1HwAM5igd0248xeACwFAEPAJYi4AHAUgQ8AFhq\n1oAPBAKqrKxURUWFOjs7p/39woULWr9+ve6++27t379/Tn0BAJmTNOCj0aja29sVCAQ0ODionp4e\nvfXWW3FtVq5cqYMHD+qZZ56Zc18AQOYkDfhQKKTy8nJ5PB7l5+ertbVVfX19cW0KCgpUW1ur/Pz8\nOfcFAGRO0oCPRCIqKSmJbbvdbkUikZQ+eD59AQDzlzTgXS7XHX/wfPoCAOYv6UrW4uJihcPh2HY4\nHJbb7U7pg+fSt6OjI/a+oaFBDQ0NKf0bAJArgsGggsHgnPokDfja2loNDQ1pZGRERUVF6u3tVU9P\nT8K2xpg77vv5gAcATHf7ye/evXtn7ZM04PPy8tTV1aWmpiZFo1G1tbXJ6/Wqu7tbkuT3+zU2Nqa6\nujpdvXpVS5Ys0YEDBzQ4OKilS5cm7AsAWBiz3mzM5/PJ5/PF7fP7/bH3hYWFcUMxs/UFACwMVrIC\ngKUIeACwFAEPAJYi4AHAUgQ8AFiKR/YByHobmzZq/MPxaft5hmxyBDyArDf+4XjCB4fzDNnkGKIB\nAEtxBr/I8dMVwEwI+EWOn64AZsIQDQBYioAHAEsR8ABgKQIeACxFwAOApZhFY6lLI5fkrYt/wMql\n8CUVavqMGwB2IuAtFVV02vTJd555x6FqADiBIRoAsBRn8ABmlWjFNEN+2Y+ABzCrRCumGfLLfgzR\nAIClOIPPQtxADEA6EPBZiBuIAUgHAh4AUrAYf1kT8GBRFJCCxfjLmoDPgETf9Nn8Lc+iKMBOBHwG\nJPqmz+ZveQB2YpokAFiKgAcASxHwAGApAh4ALEXAA4ClmEXjMO7SByBTCHiHcZc+AJnCEA0AWIqA\nBwBLzRrwgUBAlZWVqqioUGdnZ8I2Tz31lCoqKlRdXa2zZ8/G9ns8HlVVVammpkbr1q1LX9UAgFkl\nHYOPRqNqb2/XqVOnVFxcrLq6OjU3N8vr/e+NqU6cOKGLFy9qaGhIZ86c0c6dOzUwMCBJcrlcCgaD\nWrFiRWaPAgAwTdKAD4VCKi8vl8fjkSS1traqr68vLuCPHz+uJ598UpJUX1+viYkJffDBB1q9erUk\nyRiTodIBIDNsmd2WdIgmEomopKQktu12uxWJRFJu43K5tGnTJtXW1urQoUPprBsAMuaz2W2ff0Wj\nUafLmrOkZ/AulyulD5npLP3Pf/6zioqKdPnyZTU2NqqyslIbNmyY1q6joyP2vqGhQQ0NDSn9uwCQ\nK4LBoILB4Jz6JA344uJihcPh2HY4HJbb7U7aZnR0VMXFxZKkoqIiSVJBQYG2bt2qUCg0a8ADAKa7\n/eR37969s/ZJOkRTW1uroaEhjYyMaHJyUr29vWpubo5r09zcrBdffFGSNDAwoPvuu0+rV6/WjRs3\n9PHHH0uSrl+/rpMnT+rBBx+c6zEBWECfPd3r9tel8CWnS8MdSHoGn5eXp66uLjU1NSkajaqtrU1e\nr1fd3d2SJL/fry1btujEiRMqLy/Xvffeq2PHjkmSxsbG1NLSIkmamprStm3btHnz5gwfDoD5SPR0\nL4nV1YvVrLcq8Pl88vl8cfv8fn/cdldX17R+paWlOnfu3DzLAwDcKVayAoClCHgAsBQBDwCWIuAB\nwFIEPABYigd+zEOi+1VIi/OeFQDsQ8DPQ6KnMUmZmzP82SKUuH18mQCYAQG/iCRahMICFAAzYQwe\nACxFwAOApRiiAbBoJboutWrFKr32ymsOVZRdCHgAi1ai61JjXWMOVZN9CPgU2fIILwC5g4BPUaIp\nkcxgAZDNCHgAWYVfy+lDwN+G1amAs/i1nD4E/G0WenUqAGQK8+ABwFKcwQPIWbYPyRLwAHLCTBdv\n6zvrp7W1ZUiWgAeQE3Lx4i1j8ABgKQIeACxFwAOApQh4ALAUAQ8AliLgAcBSOTtN8tatW7p169a0\n/UbGgWoAIP1yNuB3/3S3jp88LpfLFdtnjNH4+LjWaI2DlQFAeuRswF+euKx7fPfof0v/N7bv2rvX\nFOmOOFgVAKQPY/AAYKmcPYMHYKdED+KW7LmB2FwQ8ACskuhB3JL9951JhIBfIJxVAFhoBPwC4awC\nuYATmewy60XWQCCgyspKVVRUqLOzM2Gbp556ShUVFaqurtbZs2fn1BeAPT47kbn9FY1GnS4tJyUN\n+Gg0qvb2dgUCAQ0ODqqnp0dvvfVWXJsTJ07o4sWLGhoa0m9/+1vt3Lkz5b654MpbV5wuIaNsP75b\nk9MXw9nC5mOT7P9vMxVJh2hCoZDKy8vl8XgkSa2trerr65PX+9+fYMePH9eTTz4pSaqvr9fExITG\nxsY0PDw8a99cMHFhQsu9y50uI2MmLkw4XUJG3bppbwg6fWyZHs6x/f+9VCQN+EgkopKSkti22+3W\nmTNnZm0TiUT03nvvzdoXQO7iulTmJQ34zy/jT8aYxXf/lvz/ydd/3vyPbg3+9yxm8sZkyscMAFnP\nJHH69GnT1NQU237uuefMCy+8ENfG7/ebnp6e2Pb9999vxsbGUuprjDFlZWVGEi9evHjxmsOrrKws\nWXwbY4xJegZfW1uroaEhjYyMqKioSL29verp6Ylr09zcrK6uLrW2tmpgYED33XefVq9erZUrV87a\nV5IuXryYrAQAwB1KGvB5eXnq6upSU1OTotGo2tra5PV61d3dLUny+/3asmWLTpw4ofLyct177706\nduxY0r4AgIXhMotxAB0AMKusuZvkwYMH5fV69cADD2jXrl1Ol5MR+/fv15IlS/Thhx86XUra/PjH\nP5bX61V1dbVaWlr00UcfOV1SWti8SC8cDusb3/iGvvKVr+iBBx7QL3/5S6dLyohoNKqamho9+uij\nTpeSVhMTE3rsscfk9Xq1du1aDQwMzNx41lH6BfCnP/3JbNq0yUxOThpjjBkfH3e4ovR79913TVNT\nk/F4POZf//qX0+WkzcmTJ000GjXGGLNr1y6za9cuhyuav6mpKVNWVmaGh4fN5OSkqa6uNoODg06X\nlTbvv/++OXv2rDHGmI8//th8+ctftur4PrN//37z+OOPm0cffdTpUtLqO9/5jjly5IgxxpibN2+a\niYmJGdtmxRn8r3/9a+3evVv5+fmSpIKCAocrSr8f/ehH+tnPfuZ0GWnX2NioJUs+/c+ovr5eo6Oj\nDlc0f59f4Jefnx9bpGeLwsJCffWrX5UkLV26VF6vV++9957DVaXX6OioTpw4oe9973uLchr3TD76\n6CO98cYb2r59u6RPr3UuW7ZsxvZZEfBDQ0N6/fXX9fWvf10NDQ36y1/+4nRJadXX1ye3262qqiqn\nS8moo0ePasuWLU6XMW8zLd6z0cjIiM6ePav6+nqnS0mrH/7wh/r5z38eO/mwxfDwsAoKCvTd735X\nX/va17Rjxw7duHFjxvYLdjfJxsZGjY2NTdv/7LPPampqSleuXNHAwIDefPNNffvb39Y77yyu1WzJ\nju/555/XyZMnY/sW2xnFTMf23HPPxcY3n332Wd111116/PHHF7q8tMuVxW7Xrl3TY489pgMHDmjp\n0qVOl5M2v//977Vq1SrV1NQoGAw6XU5aTU1N6W9/+5u6urpUV1enp59+Wi+88IL27duXuMPCjBol\n98gjj5hgMBjbLisrM//85z8drCh9/v73v5tVq1YZj8djPB6PycvLM1/84hfNBx984HRpaXPs2DHz\n0EMPmX//+99Ol5IWqS7SW8wmJyfN5s2bzS9+8QunS0m73bt3G7fbbTwejyksLDT33HOPeeKJJ5wu\nKy3ef/994/F4YttvvPGG+eY3vzlj+6wI+N/85jfmpz/9qTHGmLffftuUlJQ4XFHm2HaR9Y9//KNZ\nu3atuXz5stOlpM3NmzdNaWmpGR4eNp988ol1F1lv3bplnnjiCfP00087XUrGBYNB861vfcvpMtJq\nw4YN5u233zbGGLNnzx7zk5/8ZMa2WfHAj+3bt2v79u168MEHddddd+nFF190uqSMse3n/w9+8ANN\nTk6qsbFRkrR+/Xr96le/criq+bF9kV5/f79eeuklVVVVqaamRpL0/PPP65FHHnG4ssyw7f+5gwcP\natu2bZqcnFRZWVlscWkiLHQCAEvZdYkZABBDwAOApQh4ALAUAQ8AliLgAcBSBDwAWIqABwBLEfAA\nYKn/A1OmpzMJXIiiAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x28efe10>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def dx_dc(x, a):\n",
      "    \"\"\"dx/dc for x a solution of BA_eq.\"\"\"\n",
      "    s = np.sin(x)**2\n",
      "    return s/(a*s + 1.)\n",
      "\n",
      "\n",
      "def OG_sum(c, x, model, it):\n",
      "    \"\"\"Check OG/MZ therm limit hypothesis. \n",
      "    \n",
      "    Calculate the relevant bit of Eq. (10) of OG's notes for a fixed c & x.\n",
      "    \n",
      "    Parameters\n",
      "    ==========\n",
      "    c, x : self-explanatory\n",
      "    model : SingleImpurity type thing\n",
      "    it: iterator over the DB (or walker-generated cnf-s maybe)\n",
      "\n",
      "    Returns\n",
      "    =======\n",
      "    sum_all, num_terms\n",
      "\n",
      "    \"\"\" \n",
      "\n",
      "    sum_all = 0.\n",
      "\n",
      "    for j, row in enumerate(it):\n",
      "        # given c & fs_pairs, calculate the overlaps\n",
      "        fs_pairs = row[\"partition\"]\n",
      "       \n",
      "        # c= row[\"c\"]   # XXX: fix c instead\n",
      "        roots = model.BA_solver.find_roots(fs_pairs, c)\n",
      "        roots[\"FSfq\"] = mx.overlap(roots, model.basis_state, model.par)\n",
      "        overlap = roots[\"FSfq\"]\n",
      "    \n",
      "        # phase factor\n",
      "        phase_resc = np.pi * model.par.N / 2.    # OG uses momenta [-1, 1]\n",
      "        phase = np.sum(roots[\"roots\"]) * x / phase_resc\n",
      "        exp_factor = np.exp(1.j * phase)\n",
      "          \n",
      "        # inv jacobian\n",
      "        a = Lieb_Wu_a(model.par)\n",
      "        inv_jac = np.sum(dx_dc(roots[\"roots\"], a))\n",
      "        \n",
      "        # sum up\n",
      "        sum_all += inv_jac * exp_factor * overlap**2\n",
      "    \n",
      "        #print fs_pairs, phase, roots[\"FSfq\"]**2 * exp_factor * overlap**2\n",
      "    return sum_all, j"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 130
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compare to OG / MZ therm limit \n",
      "# fix x, c\n",
      "\n",
      "c = 2.\n",
      "\n",
      "model = SingleImpurity(par)\n",
      "par = model.par\n",
      "\n",
      "xx = np.linspace(0, 5, 11)\n",
      "\n",
      "for xp in xx:\n",
      "    print xp, OG_sum(c, xp, model, db.row_iterator(handle))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((11.297484422377682+0j), 89562)\n",
        "0.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((8.8275781056965617+6.6183611692877387j), 89562)\n",
        "1.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((3.0040186406510609+9.8360743842614902j), 89562)\n",
        "1.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-2.7883069269549532+8.7472560714180343j), 89562)\n",
        "2.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-6.061832149208942+5.0955916372704717j), 89562)\n",
        "2.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-6.6080690368521386+1.2824412354575272j), 89562)\n",
        "3.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-5.6268814944087318-1.4936966630086166j), 89562)\n",
        "3.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-4.0900318600431227-3.3171322230184659j), 89562)\n",
        "4.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-2.1872293767310134-4.4522309402491125j), 89562)\n",
        "4.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((0.023412697726102806-4.7149547288096603j), 89562)\n",
        "5.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((2.0458205565919014-3.8905541646449766j), 89562)\n"
       ]
      }
     ],
     "prompt_number": 143
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma(par), k_F(par) / initial_q(par)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 135,
       "text": [
        "(3.0, 0.8333333333333333)"
       ]
      }
     ],
     "prompt_number": 135
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compare to OG / MZ therm limit \n",
      "# fix x, c\n",
      "\n",
      "c = 0.\n",
      "\n",
      "model = SingleImpurity(par)\n",
      "par = model.par\n",
      "\n",
      "xx = np.linspace(0, 5, 11)\n",
      "\n",
      "for xp in xx:\n",
      "    print xp, OG_sum(c, xp, model, db.row_iterator(handle))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((2.6789980235870399+0j), 89562)\n",
        "0.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((2.249043145208534+1.2677629542686879j), 89562)\n",
        "1.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((1.2019297600217291+1.9647537893241183j), 89562)\n",
        "1.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((0.098594514146543202+1.8786281197089263j), 89562)\n",
        "2.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.55648047701278458+1.2547555733839399j), 89562)\n",
        "2.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.61810892649314464+0.57980817512545779j), 89562)\n",
        "3.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.32254494842755255+0.23180444650980084j), 89562)\n",
        "3.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.047020472164262257+0.26226445753043426j), 89562)\n",
        "4.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.018809210705005994+0.4478567181850685j), 89562)\n",
        "4.5 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.195755133601797+0.52760586004253129j), 89562)\n",
        "5.0 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "((-0.37571189979252195+0.41258685784268379j), 89562)\n"
       ]
      }
     ],
     "prompt_number": 142
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.linspace(0, 5, 21)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 141,
       "text": [
        "array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ,\n",
        "        2.25,  2.5 ,  2.75,  3.  ,  3.25,  3.5 ,  3.75,  4.  ,  4.25,\n",
        "        4.5 ,  4.75,  5.  ])"
       ]
      }
     ],
     "prompt_number": 141
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}