{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# oeu\n\n\n#@author: \u00d3scar N\u00e1jera\n#Created on Wed Sep 24 14:55:15 2014\neu\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "from __future__ import division, absolute_import, print_function\n\nimport slavespins.storage as st\ndef follow_cf(save, Uspan, target_cf, nup, n_tot=5.0, slsp=None):\n    \"\"\"Calculates the quasiparticle weight in single\n       site spin hamiltonian under with N degenerate half-filled orbitals \"\"\"\n\n    if slsp == None:\n        slsp = Spinon(slaves=6, orbitals=3, avg_particles=n_tot,\n                       hopping=[0.5]*6, populations = np.asarray([n_tot]*6)/6)\n\n    zet, lam, mu, mean_f = [], [], [], []\n\n    for co in Uspan:\n        print('U=', co, 'del=', target_cf)\n        res=root(targetpop, nup[-1],(co,target_cf,slsp, n_tot))\n        print(res.x)\n        if res.x>nup[-1]: break\n\n        nup.append(res.x)\n        slsp.param['populations']=population_distri(nup[-1])\n        mean_f.append(slsp.mean_field())\n        zet.append(slsp.quasiparticle_weight())\n        lam.append(slsp.param['lambda'])\n        mu.append(orbital_energies(slsp.param, zet[-1]))\n\n#    plt.plot(np.asarray(zet)[:,0], label='d={}, zl'.format(str(target_cf)))\n#    plt.plot(np.asarray(zet)[:,5], label='d={}, zh'.format(str(target_cf)))\n\n    case = save.createGroup('cf={}'.format(target_cf))\n    varis = st.setgroup(case)\n    st.storegroup(varis, Uspan[:len(zet)], zet, lam, mu, nup[1:],target_cf,mean_f)\n#    save.close()\n\ndef targetpop(upper_density, coul, target_cf, slsp, n_tot):\n    \"\"\"restriction on finding the right populations that leave the crystal\n    field same\"\"\"\n    if upper_density < 0.503: return 0.\n    trypops=population_distri(upper_density, n_tot)\n    slsp.set_filling(trypops)\n    slsp.selfconsistency(coul,0)\n    efm_free = dos_bethe_find_crystalfield(trypops, slsp.param['hopping'])\n    orb_ener = slsp.param['lambda']+ slsp.quasiparticle_weight()*efm_free\n    obtained_cf = orb_ener[5] - orb_ener[0]\n    return target_cf - obtained_cf\n\ndef population_distri(nup, n_tot=5.0):\n    nup=float(nup)\n    ndw=(n_tot-2*nup)/4\n    return np.asarray([ndw]*4+[nup]*2)\n\n#### Trying on lucas statement\n#But it doesn't work, fixing the fermion orbital energy from the free case\n#solution is bad when raising the interaction. The lagrange multiplier has\n#to high values compared to this energy and the band populations change to dras\n#tically even at low interaction\n#def targetpop(density,coul,hund, slsp):\n#    \"\"\"restriction on finding the right populations that leave the crystal\n#    field same\"\"\"\n#    ndw,nup = density\n#\n#    trypops=np.asarray([ndw]*4+[nup]*2)\n#    slsp.set_filling(trypops)\n#    slsp.selfconsistency(coul,hund)\n##    orb_ener = orbital_energies(slsp.param, slsp.quasiparticle_weight())\n##    obtained_cf = orb_ener[5] - orb_ener[0]\n#    efermi = slsp.param['lambda'] - slsp.param['orbital_e_shift']\n#    res_pops = fermion_avg(efermi, slsp.param['hopping'], 'ocupation')\n##    return [(trypops-res_pops)[-1],ndw-(5-2*nup)/4.]\n#    return [ndw-res_pops[0], nup-res_pops[5]]\n\nif __name__ == \"__main__\":\n    save = st.Dataset('dop3b5_02e.nc', 'w', format='NETCDF4')\n    deltacf = np.arange(0., 2.51, 0.1)\n    coul = -1.84*deltacf+3.5\n    for cf, ul in zip(deltacf, coul):\n#        Uspan = np.concatenate((np.arange(0, ul, 0.1), \\\n#                                np.arange(ul, 4.651, 0.02)))\n        Uspan = np.arange(0, 6, 0.1)\n        nglo = 5.02\n        nhi = [nglo/6.]\n        follow_cf(save, Uspan, cf, nhi, nglo)\n\n    save.close()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.6", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}