In [1]:

    

#from __future__ import division 
%pylab
import string
# %matplotlib qt 
#%matplotlib notebook
#%matplotlib nbagg
import matplotlib.ticker as ticker
import os

txt = 'run'
#txt = 'OldData/run'
txt2 = '/rho_'
inp ='/input_nrg.dat'
densPath = '/home/cifucito/nrgcode/TwoChNRG/Images/Density/run'



rcParams.update({'font.size': 30})
rc("text", usetex = True)
rc("font", family = "serif")
#rcParams['figure.subplot.hspace'] = 0.3
rcParams['figure.subplot.wspace'] = 0.1
from scipy.interpolate import griddata
from scipy.optimize import curve_fit
from scipy import interpolate
mpl.rcParams['lines.linewidth'] = 4
Gamma = 2.82691*0.01
pi*Gamma


    

    




Using matplotlib backend: TkAgg
Populating the interactive namespace from numpy and matplotlib






    
Out[1]:





0.08880999688359521

In [2]:

    

#read files for folder ini and dot 
#return
def readfile(ini, dot):
        output = txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow_zEQAVG.00.dat'
        
        #output = txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow.dat'
        #output = 'NotSoOld_Data/'+txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow.dat'
        #output = 'OldData/'+txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow.dat'
        infile = open(os.path.abspath(output), 'r')
        text = infile.readlines()
        
        #Data managemente
        vec = [];        vec2 = []  
        for x in text:
            #print(list(x))
            a =x.split(' ') ;            vec.append(float(a[-3]))
            vec2.append(float(a[0]))
        #if j%20 == 1 : vec = vec[::-1];
        DOS= array(vec) ; w = array(vec2) 
        return DOS , w
    
zvec = ['0.60','0.80','1.00','1.20','1.40']    
def readfileZ(ini, dot, Ztrick):
        output = txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow_zEQ'+zvec[Ztrick]+'.dat'
        #output = txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow.dat'
        infile = open(os.path.abspath(output), 'r')
        text = infile.readlines()
        #Data managemente
        vec = [];        vec2 = []  
        for x in text:
            #print(list(x))
            a =x.split(' ') ;            vec.append(float(a[-3]))
            vec2.append(float(a[0]))
        #if j%20 == 1 : vec = vec[::-1];
        DOS= array(vec) ; w = array(vec2) 
        return DOS , w    
    
#returns the parameter value of file ini  at the given line    
def readfile_Param(ini, line):
        #output = txt+repr(ini)+txt2 +repr(dot)+'_'+repr(dot)+'_OmegaRhow_zEQAVG.00.dat'
        output = txt+repr(ini)+ inp 
        #output = '/home/cifucito/nrgcode/TwoChNRG/src/Main/Run/Run2DOtM/rho_2_2_OmegaRhow.dat'
        infile = open(os.path.abspath(output), 'r')
        text = infile.readlines()
        return  float(text[line])    
    
    
#Define the log transform. 
#Inputs : w -- vector to transform ,  maxi- max(-log(w)) usually 10
def LogTransform(w, maxi):
    logw = zeros(size(w)) ; logw[w<0] = log10(-w[w<0]) ; logw[w>0] = log10(w[w>0]) 
    logw = -logw ;  #print "b" , logw
    logw[size(logw)/2:] = -logw[size(logw)/2:]+ 2*maxi
    logw =  logw -maxi
    return logw


    

In [3]:

    

# def find_nearest(array, value):
#     array = asarray(array)
#     idx = (abs(array - value)).argmin()
#     return idx

# def integrate(N, w, rho , A):
#     r = random.rand( 3, N) ;
#     sum = 0 
#     for n in range(N):
#         i1 =  find_nearest( w, A*r[0,n]) 
#         add = 1 if r[1,n]<rho[i1] else 0
#         sum = add + sum
#     return  A*sum/float(N)   

# rho, w = readfile(450, 0)
# rho1, w1 = readfile(550, 0)
# rho = 0.5*(rho +rho1)
# rho = rho*Gamma*pi ; w = w/Gamma
# plot(w, rho)


# N = 100000


# plot(w,rho)
# for i in range(20):
#     print integrate(N, w, rho,100)


    

In [4]:

    

from mpl_toolkits.axes_grid1 import make_axes_locatable


def LogColorMap(Plus , initials , dot , line , f2,ax , Method , maxplot ,PrintAxis , mult):
    W = [] ;T = [] ; Rho_values = [] ; t2=[]
    
    for ini in initials:
        param =  readfile_Param(ini, line) ; t2.append(param)
        
    maxt2 = max(t2) ; mint2 = min(t2); t2 = []
    print maxt2 -mint2
    for ini in initials:
        param =  mult*readfile_Param(ini, line)/(maxt2 - mint2)
        t2.append(param)
        rho , w = readfile(ini, dot) ; 
        #if ini == 185:
        #    rho = rho + rho[::-1]
        if Plus :
        #KEY TO READ TAKE THE AVERAGE
            rho1, w1 = readfile(ini+10, dot)
            #print size(rho) , size(rho1)
            rho = 0.5*(rho +rho1) 
            print ini
        rho = rho*pi*Gamma ; w = w/Gamma
        maxi = max(-log10(w[w>0]))
        logw = LogTransform(w, maxi) ; w = logw
        
        for term in range(size(w)):
            W.append(w[term]) ; T.append(param) ;  Rho_values.append(rho[term])
            
    points = zeros((size(W),2)) ; points[:,0] = array(W) ; points[:,1] = array(T) ; t2 = array(t2); w = array(w)
    Rho_values = array(Rho_values) #/max(Rho_values) 
    grid_x, grid_y = mgrid[min(w):max(w):1000j, min(t2):max(t2):1000j]
    grid_zA = griddata(points, Rho_values, (grid_x, grid_y), method=Method)
    #grid_zA = griddata(points, Rho_values, (grid_x, grid_y), method='linear')
    im = ax.imshow(grid_zA.T, extent=(min(w),max(w),min(t2),max(t2)),vmin=0, vmax=maxplot,cmap="gist_stern", origin='lower')
    
    divider = make_axes_locatable(ax)
    #cax = divider.append_axes("right", size="5%", pad=0.05);
    #if colorAxis: 
    #ax.set_xlim(-wlims,wlims) ;
    t2 = sort(t2) ; t2 = t2[::2]  
    if PrintAxis :
        paramTicks_scaled = linspace(min(t2),max(t2),3) ; paramTicks = paramTicks_scaled*(maxt2 - mint2)/mult
        yticks(paramTicks_scaled, (str(paramTicks[0]) ,str(paramTicks[1]) ,str(paramTicks[2]) ))
        #print 'ParamTicks =' , paramTicks_scaled , paramTicks
        
        ticks = array([-1,-10**(-3) ,-10**(-6) ,10**(-maxi), 10**(-6) ,  10**(-3),1])
        logticks = LogTransform(ticks , maxi)
        xticks(logticks , ('$-1$','$-10^{-3}$', '$-10^{-6}$','$0$', '$10^{-6}$','$10^{-3}$', '$1$'))     
        cbar = f2.colorbar(im)
    #ticks_y = ticker.FuncFormatter(lambda t2, pos: '{0:g}'.format(t2*maxt2/mult))
    #print t2 , ticks_x
    #ax.yaxis.set_major_formatter(ticks_y)
    
    #y1 = mult*0.0025/maxt2
    #ax.axhline(y=y1,linewidth=4, color='r', ls = '--')
    
    #y2 = mult*0.02/maxt2
    #ax.axhline(y=y2,linewidth=4, color='g', ls = '--')
    
    return maxi , im 

def colorMap( Plus , initials , dot , line , wlims, f2,ax , Method , maxplot ,PrintAxis , mult):
    W = [] ;T = [] ; Rho_values = [] ; t2=[]
    #mult = 80 ; mult = 0.04 
    for ini in initials:
        param =  readfile_Param(ini, line) ; t2.append(param)
        
    maxt2 = max(t2) ; mint2 = min(t2); t2 = []    
    for ini in initials:
        param =  mult*readfile_Param(ini, line)/(maxt2 - mint2)
        t2.append(param)
        rho , w = readfile(ini, dot) 
        
        #KEY TO READ TAKE THE AVERAGE
        if Plus :
        #KEY TO READ TAKE THE AVERAGE
            rho1, w1 = readfile(ini+10, dot)
            #print size(rho) , size(rho1)
            rho = 0.5*(rho +rho1) 
        rho = rho*pi*Gamma ; w = w/Gamma 
        rho = rho[abs(w)<wlims*5];  w = w[abs(w)<wlims*5]/wlims ;
        for term in range(size(w)):
            W.append(w[term]) ; T.append(param) ;  Rho_values.append(rho[term])
    points = zeros((size(W),2)) ; points[:,0] = array(W) ; points[:,1] = array(T) ; t2 = array(t2); w = array(w)
    Rho_values = array(Rho_values) #/max(Rho_values) 
    grid_x, grid_y = mgrid[min(w):max(w):1000j, min(t2):max(t2):1000j]
    grid_zA = griddata(points, Rho_values, (grid_x, grid_y), method=Method)
    #grid_zA = griddata(points, Rho_values, (grid_x, grid_y), method='linear')
    im = ax.imshow(grid_zA.T, extent=(min(w),max(w),min(t2),max(t2)),vmin=0, vmax=maxplot,cmap="gnuplot2", origin='lower')
    
    divider = make_axes_locatable(ax)
    #cax = divider.append_axes("right", size="5%", pad=0.05);
    #if colorAxis: 
    #

    ax.set_xlim(-wlims,wlims) ; t2 = sort(t2) ; t2 = t2[::2]  
    if PrintAxis :
        paramTicks_scaled = linspace(min(t2),max(t2),3) ; paramTicks = paramTicks_scaled*(maxt2 - mint2)/mult
        yticks(paramTicks_scaled, (str(paramTicks[0]) ,str(paramTicks[1]) ,str(paramTicks[2]) ))
        #print 'ParamTicks =' , paramTicks_scaled , paramTicks
        wTicks_scaled = linspace(-1,1,3) ;  wTicks = wTicks_scaled*wlims
        xticks(wTicks_scaled, (str(wTicks[0]) ,str(wTicks[1]) ,str(wTicks[2])))  
        cbar = f2.colorbar(im, ax=ax)
        cbar.set_ticks([0, 0.5*maxplot, maxplot])
        cbar.set_ticklabels(['$0$', '$'+str(0.5*maxplot)+'$', '$>'+str(maxplot)+'$'])
        #print 'wTicks =' , wTicks_scaled , wTicks    
        
    #ticks_y = ticker.FuncFormatter(lambda t2, pos: '{0:g}'.format(t2*maxt2/mult))
    #print t2 , ticks_x
    #ax.yaxis.set_major_formatter(ticks_y)
    
    #y1 = mult*0.0025/maxt2
    #ax.axhline(y=y1,linewidth=4, color='r', ls = '--')
    
    #y2 = mult*0.02/maxt2
    #ax.axhline(y=y2,linewidth=4, color='g', ls = '--')
    #return max(Rho_values) , im ]
    return im 


#Plots 2 or 4 Dos plots
#Input : LogPlot (boolean) : true - logPlot , false - Normal plot
#        numplots (1,2) : 1 - plots only dots 0 and 1 , 2 - plots 4 dots
#        initials - array of file numbers
#        line - line of changing parameter 
#        wlims - defines the w limits if logPlot = false 
#        paramName - name of changing param 
#        Method: Method of approximation - 'Linear ' , 'Nearest'
#        MaxRho: Maximum Permited DOS 
#        mult : Allows to change the dimensions of the plot if required
def color4Dots( Plus ,  LogPlot, numPlots,  initials ,  line , wlims ,paramName, Method , MaxRho , mult ):  
    f2, axarr = subplots(numPlots, 2, sharex='col', sharey='row', figsize=(12, 2.5*numPlots))
    maxplot = 0 ; dots = range(2,4) ;  dots = range(2*numPlots) 
    t2 = [] 
    for ini in initials:
        param =  readfile_Param(ini, line)/Gamma ;  t2.append(param)
    print t2    
    maxt2 = (max(t2) -min(t2)) ; t2 = sort(mult*array(t2)/maxt2 ) ;  t2 = t2[::2] 
    print maxt2 , t2
    #for dot in dots:
    #    for ini in initials:
    #        rho , w = readfile(ini, dot) 
    #        maxrho = max(rho[abs(w)<wlims]) 
    #        #print maxrho , ini , dot
    #        maxplot = maxrho if maxrho>maxplot else maxplot 
    for dot in dots:
        spin = '\uparrow'
        if dot%2 == 1 : spin = '\downarrow'; 
        ax = axarr[dot-dots[0]] if numPlots == 1 else axarr[dot/2 , dot%2 ]
        
        if LogPlot:
            maxi , im = LogColorMap( Plus , initials , dot , line , f2, ax, Method, MaxRho , False, mult)
            ticks = array([-10,-10**(-4),10**(-maxi),  10**(-4),10])
            wTicks_scaled = LogTransform(ticks , maxi)
            wTicksLabels = ['$-10$','$-10^{-4}$','$0$', '$10^{-4}$', '$10$']
        else:
            im = colorMap(Plus , initials , dot , line , wlims,f2, ax, Method ,MaxRho  , False , mult)
            wTicks_scaled = linspace(-1,1,3) ;        wTicks = wTicks_scaled*wlims
            wTicksLabels = ['$'+str(wTicks[0])+'$', '$'+str(wTicks[1])+'$', '$'+str(wTicks[2])+'$']
        
        paramTicks_scaled = linspace(min(t2),max(t2),3) ;
#         paramTicks_scaled = linspace(min(t2)+ 0.25/Gamma ,max(t2)+ 0.25/Gamma ,3) 
       
        paramTicks =  0.25/Gamma+paramTicks_scaled*maxt2/mult   
        paramTicks =  paramTicks_scaled*maxt2/mult 
        print  paramTicks_scaled, paramTicks
        
        setp(axarr, xticks = wTicks_scaled, xticklabels= wTicksLabels ,
        yticks=paramTicks_scaled ,yticklabels=['$'+"{0:.1f}".format(paramTicks[0])+'$', '$'+"{0:.1f}".format(paramTicks[1])+'$', '$'+"{0:.1f}".format(paramTicks[2])+'$'] )
        
        #setting ticks
        
        title = '$\\rho_{1'+ spin +'}  =\\rho_{2'+ spin +'} $ ' if numPlots == 1 else '$\\rho_{'+str(dot/2 + 1)+ spin +'}$'
        ax.set_title(title)
        if dot >= (numPlots*2-2) : ax.set_xlabel('$\omega/\Gamma_1'+ '$')
        if dot%2 == 0 : ax.set_ylabel('$'+paramName+'$') 
            
    cbar = f2.colorbar(im, ax=axarr.ravel().tolist() )
    cbar.set_ticks([0, 0.5*MaxRho, MaxRho])
    cbar.set_ticklabels(['$0$', '$'+str(0.5*MaxRho)+'$', '$>'+str(MaxRho)+'$'])
    #cbar.ax.set_yticklabels(['< -1', '0', '> 1'])
    #f2.subplots_adjust(right=0.8) ticks=[0, 0.5*MaxRho, MaxRho]
    #cbar_ax = f2.add_axes([0.85, 0.15, 0.05, 0.7])
    #f2.colorbar(im, cax=cbar_ax)


    

In [ ]:

    




    

In [69]:

    

dot=4
line = 13
Method = 'linear' #; Method = 'nearest'
initials = array([450,451,452,453,454,455,456,457,458,459])
#initials = range(400,409)

# initials = array([180,181,182,183,184,190 ,191,192,193,194,195,186,187,188,189])
f2, ax = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 8))

maxplot = 0
for ini in initials:
     rho , w = readfile(ini, dot) ;     maxrho = max(rho[abs(w)<wlims*2]) ;    maxplot = maxrho if maxrho>maxplot else maxplot 
LogColorMap( initials , dot , line , f2, ax, Method, maxrho , True , 10)
print maxrho  , maxplot


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-69-639f454b596c> in <module>()
     12 for ini in initials:
     13      rho , w = readfile(ini, dot) ;     maxrho = max(rho[abs(w)<wlims*2]) ;    maxplot = maxrho if maxrho>maxplot else maxplot
---> 14 LogColorMap( initials , dot , line , f2, ax, Method, maxrho , True , 10)
     15 print maxrho  , maxplot

TypeError: LogColorMap() takes exactly 10 arguments (9 given)

In [37]:

    

wlims = 0.01
dot=0
line = 12
Method = 'linear' #; Method = 'nearest'
initials = array([450,451,452,453,454,455,456,457,458,459])
initials = range(820,829)
#initials = range(100,109)
f2, ax = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 8))

maxplot = 0
for ini in initials:
     rho , w = readfile(ini, dot) ;     maxrho = max(rho[abs(w)<wlims*2]) ;    maxplot = maxrho if maxrho>maxplot else maxplot 
colorMap( True, initials , dot , line , wlims,f2, ax, Method, 0.2 , True , 1.5)
# colorMap( Plus , initials , dot , line , wlims, f2,ax , Method , maxplot ,PrintAxis , mult)
print maxrho  , maxplot

#color4Dots( initials ,  line , wlims)


    

    




7.98560804909 7.98560804909

In [44]:

    

Method = 'linear' #; Method = 'nearest';  
line = 4
xtitle = [None]*23
xtitle[4] = '\Delta\epsilon_{1}/\Gamma_1' ; xtitle[12] = 't_1/\Gamma = t_2/\Gamma' ; #xtitle[12] = 't_1/\Gamma' ;
xtitle[16] = '\Delta\epsilon_{2}/\Gamma_1'; xtitle[15] = '\epsilon_{d2}/\Gamma = \\frac{-U_2}{\Gamma}';xtitle[13] = 't_2/\Gamma'
xtitle[17] = 't_{dots}/\Gamma ' ; xtitle[18] = '\epsilon_m/\Gamma ' ;xtitle[19] = '\phi_1' ;  xtitle[20] = '\phi_2' ;xtitle[21] = '\eta_1' ;
xtitle[22] = '\eta_1' 
rcParams.update({'font.size': 30})
initials = array([450, 812,814,815,451,452,453,454,455,456,457,458,459,817,818])
#initials = array([460,461,462,463,464,465,466,467,468,469])
#initials = array([470,471,472,473,474,475,476,477,478,479])

# initials = range(430,439)
# initials = range(820,829)
# initials = range(860,869)
# initials = range(960,970)
# initials = range(830,839)
#initials = range(830,837)
#initials = range(460,469)
#initials = array([185,190 ,191,192,193,194,186,187,188,189])
#initials = range(195,199)
#initials = range(800,809)
#initials = array([180,181,182,183,184,185,186,187,188,189])
#initials = range(4)
#initials = range(100,109)
#initials = range(420,429) ; initials = concatenate(initials,range(480,484))
#initials = range(400,409)
#initials = array([420,421,422,423,424,425,426,427,428,429,480])
#initials = range(130,134)
#initials = range(410,415)
#initials = array([100,130,131,132,133,134])
#initials = range(140,144)
#initials = range(110,119)
#initials = range(160,164)
#initials = range(420,429)
#initials = range(430,434)
initials = range(380,389)
# initials = range(140,144)
line = 16
wlims =0.5
#Normal 2D
color4Dots(False , False, 2,initials ,  line , wlims, xtitle[line], Method ,1,1.5)

#Log 2D
color4Dots(False , True, 2,initials ,  line , wlims, xtitle[line], Method ,1,12)


    

    




[-8.843578324035786, -7.843578324035786, -6.843578324035785, -5.843578324035785, -4.843578324035785, -3.843578324035785, -2.843578324035785, -1.8435783240357846, -0.8435783240357847]
8.0 [-1.65817094 -1.28317094 -0.90817094 -0.53317094 -0.15817094]
[-1.65817094 -0.90817094 -0.15817094] [-8.84357832 -4.84357832 -0.84357832]
[-1.65817094 -0.90817094 -0.15817094] [-8.84357832 -4.84357832 -0.84357832]
[-1.65817094 -0.90817094 -0.15817094] [-8.84357832 -4.84357832 -0.84357832]
[-1.65817094 -0.90817094 -0.15817094] [-8.84357832 -4.84357832 -0.84357832]
[-8.843578324035786, -7.843578324035786, -6.843578324035785, -5.843578324035785, -4.843578324035785, -3.843578324035785, -2.843578324035785, -1.8435783240357846, -0.8435783240357847]
8.0 [-13.26536749 -10.26536749  -7.26536749  -4.26536749  -1.26536749]
0.2261528
[-13.26536749  -7.26536749  -1.26536749] [-8.84357832 -4.84357832 -0.84357832]
0.2261528
[-13.26536749  -7.26536749  -1.26536749] [-8.84357832 -4.84357832 -0.84357832]
0.2261528
[-13.26536749  -7.26536749  -1.26536749] [-8.84357832 -4.84357832 -0.84357832]
0.2261528
[-13.26536749  -7.26536749  -1.26536749] [-8.84357832 -4.84357832 -0.84357832]

In [33]:

    

j = 458 ; i = 1
rho , vec2 = readfile(j,i)
vec0 = [] ;         maxi = []         
            #KEY TO READ TAKE THE AVERAGE
rho1, w1 = readfile(j+100, i)
vec = 0.5*(rho +rho1) 
#vec = 0.5*(rho +rho[::-1])
vec0mas = vec[ vec2 > 0] ; vec0mas = vec0mas[0]
vec0menos = vec[ vec2 < 0] ; vec0menos = vec0menos[-1]
vec0.append(0.5*(vec0mas + vec0menos)) ; maxi.append(max(vec))
print maxi


    

    




[1.3075933222493286]

In [35]:

    

2.6086911593015083


    

    
Out[35]:





'13.95'

In [45]:

    

from matplotlib import colors as mcolors
rcParams.update({'font.size': 30})
f2, ax = subplots(1, 1, sharex='col', sharey='row', figsize=(4, 2))
t2 = []    
initials = array([458,461,462,463,464,465,466,467])
initials = range(330,338)
# initials = range(110,120)
initials = range(380,389)
# initials = array([458,461,462,463,464,465,466,467,468,469])
line = 16
# initials = range(360,368)
Markers = ["^","v"]

for ini in initials: 
    t2.append(readfile_Param(ini, line))
t2 = array(t2)
print t2
for i in range(4):
        vec0 = [] ;         maxi = [] ;no = False
        #for j in range(40,46):
        for j in initials:
            rho , vec2 = readfile(j, i)
            
            #KEY TO READ TAKE THE AVERAGE
            if no:
#                 print i , no
                rho1, w1 = readfile(j+10, i)
                vec = 0.5*(rho +rho1) 
                no = False
                
                
            vec = 0.5*(rho +rho[::-1])
            vec0mas = vec[ vec2 > 0] ; vec0mas = vec0mas[0]
            vec0menos = vec[ vec2 < 0] ; vec0menos = vec0menos[-1]
            vec0.append(0.5*(vec0mas + vec0menos)) ; maxi.append(max(vec))
        vec0 = array(vec0)
        spin = '\uparrow' if i%2 == 0 else '\downarrow'; 
#         ax.semilogy(t2/Gamma + 0.25/Gamma  ,vec0*pi*Gamma, label= '$'+spin+'$',  marker=Markers[i%2] , markersize=10, linewidth=3.0)
             
        title = '$\\rho_{'+str(i/2 +1 )+ spin +'}$' 
        ax.plot(t2/Gamma + 0.25/Gamma , vec0*pi*Gamma , label= title + '$(0)$', marker=Marks[i/2], markersize=11, linewidth=3.0)
        if i%2 == 0 : 
            new0 = vec0
        else:
            spin = '\uparrow' if i%2 == 0 else '\downarrow'; 
            
            title = '$\\rho_{'+str(i/2 +1 )+ spin +'}$'    
            ydata = array(new0)/array(vec0)
            ydata = array(new0) ;  t2_new = sort(t2) ;  ydata_new = ydata[argsort(t2)] ;
            Marks = ["s","h"]
            cols = ["teal" , "indigo"]
#             ax.plot(t2_new/Gamma + 0.249/Gamma  ,ydata_new, label= 'Dot' + str(i/2 + 1),  marker=Marks[i/2],
#                     markersize=11, linewidth=3.0 , color = cols[i/2] )
        
            ydata_new = ydata[argsort(t2)] ;
ax.set_xlabel('$\Delta \epsilon_2/\Gamma_1 $') 
# ax.set_xlim(0,7) ; ax.set_ylim(0,3)


# ax.set_ylabel('$\\rho_{\uparrow}(0)/\\rho_{\downarrow}(0)$')
ax.set_ylabel('$\\rho(0)\pi \Gamma_1$')
# ax.axvline(x=5.02, color = 'r', linestyle = ':')
# ax.axvline(x=.02, color = 'r', linestyle = '--')
ax.legend(loc='upper left', fontsize=18)
#ax.set_xlim(0,5.4)


    

    




[-0.25      -0.2217309 -0.1934618 -0.1651927 -0.1369236 -0.1086545
 -0.0803854 -0.0521163 -0.0238472]






    
Out[45]:





<matplotlib.legend.Legend at 0x7f9d8c5b4a10>

In [72]:

    

ax.legend(loc='upper center', fontsize=18)


    

    
Out[72]:





<matplotlib.legend.Legend at 0x7f7c481a7350>

In [ ]:

    




    

In [28]:

    

Markers = ["^","v"]
def plotFermi(initials,line,limsx, limsy, Param_name):
    f2, axarr2 = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 8))
    #plot the relation between both DOS (up/Down)
    f3, axarr3 = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 8))
    t2 = []    
    for ini in initials: 
        t2.append(readfile_Param(ini, line))
    t2 = array(t2)
    
    print t2
    for i in range(2):
        vec0 = [] ;         maxi = []
        #for j in range(40,46):
        for j in initials:
            rho , vec2 = readfile(j, i)
            
            #KEY TO READ TAKE THE AVERAGE
            rho1, w1 = readfile(j+10, i)
            vec = 0.5*(rho +rho1) ; vec2 = vec2/Gamma
            vec = 0.5*(rho +rho[::-1])
            vec0mas = vec[ vec2 > 0] ; vec0mas = vec0mas[0]
            vec0menos = vec[ vec2 < 0] ; vec0menos = vec0menos[-1]
            vec0.append(0.5*(vec0mas + vec0menos)) ; maxi.append(max(vec))
        vec0 = array(vec0)
        spin = '\uparrow' if i%2 == 0 else '\downarrow';  
#         title = '$'+str(i/2 +1 )+ spin +'}$'    
        title = '$'+ spin +'$'    
        axarr2.plot(vec2/Gamma + 0.25/Gamma,vec0*pi*Gamma , label= title + '$(0)$', marker=Markers[i] , markersize=10, linewidth=3.0)
        if i%2 == 0 : 
            new0 = vec0
        else: 
            ydata = array(new0)/array(vec0) ;  t2_new = sort(t2) ;  ydata_new = ydata[argsort(t2)] ; 
            t2_new =   t2_new/Gamma - 0.25/Gamma
            axarr3.plot(t2_new/Gamma + 0.25/Gamma ,ydata_new, label= 'Dot' + str(i/2 + 1),  marker="o" , markersize=10, linewidth=3.0)
            ydata_new = ydata[argsort(t2)] ; 
            #axarr3.plot(t2_new,ydata_new,'-', label= 'Dot' + str(i/2 + 1) )
            #tnew =  np.linspace(t2[0], t2[-1], 50) ; f = interpolate.interp1d(t2, ydata, kind='linear')
            #axarr3.plot(tnew,f(tnew),'-' )

    axarr2.set_xlabel('$'+Param_name+'$')
    axarr2.set_ylabel('$\\rho_1(0)$')
    axarr2.legend(fontsize=20 , loc='upper right' )
    axarr2.set_ylim([0,1])
    #axarr2.set_ylim([0.0,4.5])
    #axarr3.set_xlabel('$\epsilon_{d2}$',fontsize=25)
    axarr3.set_xlabel('$'+Param_name+'$')
    axarr3.set_ylabel('$\\rho_{\uparrow}(0)/\\rho_{\downarrow}(0)$')     
    axarr3.legend(fontsize=27, loc='upper right')  
    axarr3.set_xlim([t2[0]/Gamma + 0.25/Gamma,t2[-1]/Gamma + 0.25/Gamma])
    t2 = sort(t2) 
    paramTicks = linspace(t2[1]/Gamma + 0.25/Gamma,t2[-2]/Gamma + 0.25/Gamma,3) 
    axarr3.set_xticks(paramTicks) ;    axarr3.set_xticklabels(['$'+str(paramTicks[0])+'$', '$'+str(paramTicks[1])+'$', '$'+str(paramTicks[2])+'$'])
    axarr3.set_yticks(array([1,2,3,4,5])) ;  axarr3.set_yticklabels( ['$1.0$', '$2.0$', '$3.0$', '$4.0$' ,'$5.0$'])
    #axarr3.set_yticks(linspace(1,1.5,2))
    #axarr3.set_yticklabels(['$1$', '$1.5$', '$2$'])
    #axarr3.set_xticks(linspace(-2,2,5), ('$-2$','$-1$', '$0$','$1$', '$2$'))
    #axarr3.set_ylim([0,3])


    

In [29]:

    

#xtitle = 't_1/D = t_2/D' ; xtitle = '\epsilon_{d2}/D'; #xtitle = '\epsilon_{d2}/D = \\frac{-U_2}{2D}';# xtitle = 't_2/D'
rcParams.update({'font.size': 45})
initials = array([450,451,452,453,454,455,456,457,458,459])
initials = array([460,461,462,463,464,465,466,467,468])
initials = range(400,409)
initials = range(420,429)
initials = range(110,119)
initials = range(120,129)
initials = array([420,421,422,423,424,425,426,427,428,429,480,481,482,483])
#initials = array([420,421,422,423,424,425,426,427,428,429])
#initials = array([470,471,472,473,474,475,476,477,478,479])
initials = range(180,189)
initials = range(640,649)
initials = range(960,970)
initials = range(460,469)
initials = range(330,338)
initials = range(370,379)
# initials = range(110,119)
line = 13
plotFermi(initials,line,0.1,2.5,xtitle[line])


    

    




[ 0.0565382  0.0565382  0.0565382  0.0565382  0.0565382  0.0565382
  0.0565382  0.0565382  0.0565382]






    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-29-39b19712b431> in <module>()
     18 # initials = range(110,119)
     19 line = 13
---> 20 plotFermi(initials,line,0.1,2.5,xtitle[line])

<ipython-input-28-59cf8df73f0e> in plotFermi(initials, line, limsx, limsy, Param_name)
     27 #         title = '$'+str(i/2 +1 )+ spin +'}$'
     28         title = '$'+ spin +'$'
---> 29         axarr2.plot(vec2/Gamma + 0.25/Gamma,vec0*pi*Gamma , label= title + '$(0)$', marker=Markers[i] , markersize=10, linewidth=3.0)
     30         if i%2 == 0 :
     31             new0 = vec0

/usr/lib/python2.7/dist-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1716                     warnings.warn(msg % (label_namer, func.__name__),
   1717                                   RuntimeWarning, stacklevel=2)
-> 1718             return func(ax, *args, **kwargs)
   1719         pre_doc = inner.__doc__
   1720         if pre_doc is None:

/usr/lib/python2.7/dist-packages/matplotlib/axes/_axes.pyc in plot(self, *args, **kwargs)
   1370         kwargs = cbook.normalize_kwargs(kwargs, _alias_map)
   1371 
-> 1372         for line in self._get_lines(*args, **kwargs):
   1373             self.add_line(line)
   1374             lines.append(line)

/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in _grab_next_args(self, *args, **kwargs)
    402                 this += args[0],
    403                 args = args[1:]
--> 404             for seg in self._plot_args(this, kwargs):
    405                 yield seg
    406 

/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in _plot_args(self, tup, kwargs)
    382             x, y = index_of(tup[-1])
    383 
--> 384         x, y = self._xy_from_xy(x, y)
    385 
    386         if self.command == 'plot':

/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in _xy_from_xy(self, x, y)
    241         if x.shape[0] != y.shape[0]:
    242             raise ValueError("x and y must have same first dimension, but "
--> 243                              "have shapes {} and {}".format(x.shape, y.shape))
    244         if x.ndim > 2 or y.ndim > 2:
    245             raise ValueError("x and y can be no greater than 2-D, but have "

ValueError: x and y must have same first dimension, but have shapes (492,) and (9,)

In [ ]:

    




    

In [22]:

    

#Define the log transform. 
#Inputs : w -- vector to transform ,  maxi- max(-log(w)) usually 10
def LogTransform(w, maxi):
    logw = zeros(size(w)) ; logw[w<0] = log10(-w[w<0]) ; logw[w>0] = log10(w[w>0]) 
    logw = -logw ;  #print "b" , logw
    logw[size(logw)/2:] = -logw[size(logw)/2:]+ 2*maxi
    logw =  logw -maxi
    return logw 
        #    return logVector
def LogPlotDOS(PlusMinus,Critical,initials , dot , param , param_name , MajoranaSig):
    for i in initials:
        rho, w = readfile(i, dot)
        #rho = rho*Gamma*pi*1.06; w=w/Gamma
        #rho1 = rho[::-1] ; rho = (rho +rho1)*0.5
        if PlusMinus:
            rho1, w1 = readfile(i+10, dot)
            #rho1 = rho1*Gamma*pi*1.06; w1=w1/Gamma
            rho = 0.5*(rho +rho1) ; w = 0.5*(w+w1)
        #maxi = max(-log(w[w>0])) ; print maxi
        maxi = max(-log10(w[w>0]))
        logw = LogTransform(w , maxi)
        ticks = array([-1,-10**(-3) ,-10**(-6) ,10**(-maxi), 10**(-6) ,  10**(-3),1])
        #ticks = array([-1,-10**(-3) ,-10**(-6) ,-10**(-9),10**(-maxi),10**(-9), 10**(-6) ,  10**(-3),1])
        logticks = LogTransform(ticks , maxi)
        xticks(logticks , ('$-1$','$-10^{-3}$', '$-10^{-6}$','$0$', '$10^{-6}$','$10^{-3}$', '$1$'))
        #xticks(logticks , ('$-1$','$-10^{-3}$', '$-10^{-6}$','$-10^{-9}$','$0$','$10^{-9}$', '$10^{-6}$','$10^{-3}$', '$1$'))
        #rhoDw, w3 = readfile(i, dot+1) ; rhoDw1, w3 = readfile(i+100, dot+1) ;
        #rhoDW1 = rhoDw[::-1] ; rhoDw = (rhoDw +rhoDW1)*0.5 ; rhoDw = rhoDw[::-1] 
       # rhoDw1, w4 = readfile(i+100, dot+1)
        #rhoDw = 0.5*(rhoDw +rhoDw1) 
        #logw = zeros(size(w)) ; logw[w<0] = log10(-w[w<0]) ; logw[w>0] = log10(w[w>0]) 
        #logw = -logw ; 
         #print "b" , logw
        spin = '\uparrow' if dot%2==0 else '\downarrow'
        pr = readfile_Param(i, param)
        lb = 'QD$'+str(dot/2 + 1)+spin +','+ param_name +'= '+str(pr)#+ 'D$ '
        #lb = spin +','+ param_name +'= '+str(pr)#+ 'D$ '
        lb = 'QD$'+str(dot/2 + 1)+spin +','
        style = '-' if dot%2==0 else '--'
        if Critical:
            deltaW = logw[-1]-logw[-2]
            print abs(deltaW)
            logrho = log10(rho)
#             deltaRho = (logrho[1:]-logrho[:-1])/abs(deltaW)
#             plot(logw[:-1], deltaRho, label = lb , ls = style)
            semilogy(logw, rho, label = lb , ls = style)
            siz = size(logw)/2 ;  siz1= siz+1
#             print "Critical: " ,deltaRho[siz:]
            #print (logrho[siz1:]-logrho[siz:-1])/deltaW
        else:
            if MajoranaSig:
                plot(logw, rho/rhoDw, label = lb , ls = style)
            else :    
            #if pr == 0:
                print 'a'
                plot(logw, rho, label = lb , ls = style)
            #else:    
            #plot(logw, rho, label = lb , ls = style)
            #lb = 'QD'+str(dot/2 +1) +','+ param_name +'= '+str(pr)#+ 'D$ '    
            #plot(logw, rho, label = lb , ls = style)
        #if pr == 0:
        #    lb = '$'+param_name +'= '+ str(pr) + 'D$ ' 
        #    plot(logw, rho,label = lb ,color = 'k', marker = 'o',  ls = 'None' , markersize=5 )
        #else:    
        
        #    mk = '^' if dot%2==0 else 'v'
        #    #mk = "*" if dot%2==0 else 'H'
        #    if i%2==0:
        #        cl = 'r'# if dot%2==0 else 'y'
        #    else:    
        #        cl = 'b'# if dot%2==0 else 'g'
        #plot(logw, rho,label = lb ,color = cl, marker = mk,  ls = 'None' , markersize=-20/log10(pr))
        #plot(logw, rho, label = lb , ls = style)
        
        #plot(logw, rhoDw, label = lb , ls = style)
        #plot(logw, rhoDw, label = lb + 'dw' , ls = style)
            #logw[size(logw)/2:]= logw[size(logw)/2:]+1
    return logticks
    
#print logticks , str(readfile_Param(i, param))
#; ylabel('$\\rho_{2}$')
#title('Dot 1  $(t_1=t_2=0.005)$')

initials = [ 451,458] ; 
initials = [ 430,432,438] ;
#initials = [451,458 ] ;
#initials = [460,466,468 ] ;
#initials = [430,434 ] ;
#initials = [400,403 ,407 ] ;
#initials = [410 , 411 , 414] ; 
#initials = [110,111]
initials = [154]
#initials = [2]
#initials = [450]
initials = [429,484]
initials = [497]
initials = [110,112,114]
#initials = [160,162,164]
initials = [369]

#initials = [873 ,879]
param = 4
xtitle = [None]*23
xtitle[4] = '\epsilon_{d1}/D' ; xtitle[12] = 't_1/D = t_2/D' ; xtitle[12] = 't_1/D ' ; 
xtitle[16] = '\epsilon_{d2}/D'; xtitle[15] = '\epsilon_{d2}/D = \\frac{-U_2}{2D}';xtitle[13] = 't_{1,2}/D'
xtitle[17] = 't_{dots}/D ' ; xtitle[18] = '\epsilon_m ' ;xtitle[19] = '\phi_1' ;  xtitle[20] = '\phi_2' ;xtitle[21] = '\eta_2' ;
xtitle[22] = '\eta_2' 
xlabel('$\omega/\Gamma_1$')
ylabel('$\\rho_{1}=\\rho_{2}$') ; ylabel('$\\rho$')
logticks = LogPlotDOS(False, True, initials,2,param , xtitle[param] , False)
logticks = LogPlotDOS(False, True ,  initials,3,param , xtitle[param] , False)
#logticks = LogPlotDOS(True, False ,  initials,2,param , xtitle[param] , False)
#logticks = LogPlotDOS(True, False ,  initials,3,param , xtitle[param] , False)
#logticks = LogPlotDOS(True,False , initials,1,param , xtitle[param])
#logticks = LogPlotDOS(True,False ,  initials,3 ,param , xtitle[param])
#logticks=LogPlotDOS(True ,False ,initials[:-1] , 1,param ,xtitle[param])
#LogPlotDOS(551 , param)

legend(fontsize=18)

#title('Dot 2 ($t_1=t_2=0.0025D$)')


    

    




0.0397940008672
0.0397940008672






    




/usr/lib/python2.7/dist-packages/matplotlib/ticker.py:2206: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
  "Data has no positive values, and therefore cannot be "






    
Out[22]:





<matplotlib.legend.Legend at 0x7f9d8c806a90>






    




/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py:1398: UserWarning: aspect is not supported for Axes with xscale=linear, yscale=log
  'yscale=%s' % (xscale, yscale))
/usr/lib/python2.7/dist-packages/matplotlib/ticker.py:2176: RuntimeWarning: overflow encountered in power
  ticklocs = b ** decades
/usr/lib/python2.7/dist-packages/matplotlib/ticker.py:1099: RuntimeWarning: invalid value encountered in double_scalars
  coeff = np.round(x / b ** exponent)

In [23]:

    

xticks(logticks , ('$-1$','$-10^{-3}$', '$-10^{-6}$','$0$', '$10^{-6}$','$10^{-3}$', '$1$'))
legend(fontsize=17, loc='upper right')
#xticks(logticks , ('$-1$','$-10^{-3}$', '$-10^{-6}$','$-10^{-9}$','$0$','$10^{-9}$', '$10^{-6}$','$10^{-3}$', '$1$'))


    

    




Exception in Tkinter callback
Traceback (most recent call last):
  File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 1544, in __call__
    return self.func(*args)
  File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 595, in callit
    func(*args)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_tkagg.py", line 323, in idle_draw
    self.draw()
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_tkagg.py", line 304, in draw
    FigureCanvasAgg.draw(self)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py", line 430, in draw
    self.figure.draw(self.renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/figure.py", line 1299, in draw
    renderer, self, artists, self.suppressComposite)
  File "/usr/lib/python2.7/dist-packages/matplotlib/image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py", line 2437, in draw
    mimage._draw_list_compositing_images(renderer, self, artists)
  File "/usr/lib/python2.7/dist-packages/matplotlib/image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/legend.py", line 768, in draw
    bbox = self._legend_box.get_window_extent(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 263, in get_window_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 385, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 456, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 385, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 456, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 829, in get_extent
    bbox, info, d = self._text._get_layout(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/text.py", line 317, in _get_layout
    ismath=ismath)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py", line 226, in get_text_width_height_descent
    s, fontsize, renderer=self)
  File "/usr/lib/python2.7/dist-packages/matplotlib/texmanager.py", line 602, in get_text_width_height_descent
    dvifile = self.make_dvi(tex, fontsize)
  File "/usr/lib/python2.7/dist-packages/matplotlib/texmanager.py", line 400, in make_dvi
    exc.output.decode("utf-8"))))
RuntimeError: LaTeX was not able to process the following string:
'QD$2\\\\uparrow,'

Here is the full report generated by LaTeX:
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=latex)
 restricted \write18 enabled.
entering extended mode
(./ee45ce73b3895efa1c12465ea0a86a74.tex
LaTeX2e <2017-04-15>
Babel <3.18> and hyphenation patterns for 84 language(s) loaded.
(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
(/usr/share/texlive/texmf-dist/tex/latex/base/size10.clo))
(/usr/share/texlive/texmf-dist/tex/latex/type1cm/type1cm.sty)
(/usr/share/texlive/texmf-dist/tex/latex/base/textcomp.sty
(/usr/share/texlive/texmf-dist/tex/latex/base/ts1enc.def))
(/usr/share/texlive/texmf-dist/tex/latex/geometry/geometry.sty
(/usr/share/texlive/texmf-dist/tex/latex/graphics/keyval.sty)
(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifpdf.sty)
(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifvtex.sty)
(/usr/share/texlive/texmf-dist/tex/generic/ifxetex/ifxetex.sty)

Package geometry Warning: Over-specification in `h'-direction.
    `width' (5058.9pt) is ignored.


Package geometry Warning: Over-specification in `v'-direction.
    `height' (5058.9pt) is ignored.

)
No file ee45ce73b3895efa1c12465ea0a86a74.aux.
(/usr/share/texlive/texmf-dist/tex/latex/base/ts1cmr.fd)
*geometry* driver: auto-detecting
*geometry* detected driver: dvips
! Extra }, or forgotten $.
l.12 ...00000}{22.500000}{\rmfamily QD$2\uparrow,}
                                                  
! Missing $ inserted.
<inserted text> 
                $
l.13 \end{document}
                   
[1] (./ee45ce73b3895efa1c12465ea0a86a74.aux) )
(\end occurred inside a group at level 1)

### simple group (level 1) entered at line 12 ({)
### bottom level
(see the transcript file for additional information)
Output written on ee45ce73b3895efa1c12465ea0a86a74.dvi (1 page, 336 bytes).
Transcript written on ee45ce73b3895efa1c12465ea0a86a74.log.
 








    
Out[23]:





<matplotlib.legend.Legend at 0x7f9d8c1496d0>






    




Exception in Tkinter callback
Traceback (most recent call last):
  File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 1544, in __call__
    return self.func(*args)
  File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 595, in callit
    func(*args)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_tkagg.py", line 323, in idle_draw
    self.draw()
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_tkagg.py", line 304, in draw
    FigureCanvasAgg.draw(self)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py", line 430, in draw
    self.figure.draw(self.renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/figure.py", line 1299, in draw
    renderer, self, artists, self.suppressComposite)
  File "/usr/lib/python2.7/dist-packages/matplotlib/image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py", line 2437, in draw
    mimage._draw_list_compositing_images(renderer, self, artists)
  File "/usr/lib/python2.7/dist-packages/matplotlib/image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/matplotlib/legend.py", line 768, in draw
    bbox = self._legend_box.get_window_extent(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 263, in get_window_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 385, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 456, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 385, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 256, in get_extent
    w, h, xd, yd, offsets = self.get_extent_offsets(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 456, in get_extent_offsets
    for c in self.get_visible_children()]
  File "/usr/lib/python2.7/dist-packages/matplotlib/offsetbox.py", line 829, in get_extent
    bbox, info, d = self._text._get_layout(renderer)
  File "/usr/lib/python2.7/dist-packages/matplotlib/text.py", line 317, in _get_layout
    ismath=ismath)
  File "/usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py", line 226, in get_text_width_height_descent
    s, fontsize, renderer=self)
  File "/usr/lib/python2.7/dist-packages/matplotlib/texmanager.py", line 602, in get_text_width_height_descent
    dvifile = self.make_dvi(tex, fontsize)
  File "/usr/lib/python2.7/dist-packages/matplotlib/texmanager.py", line 400, in make_dvi
    exc.output.decode("utf-8"))))
RuntimeError: LaTeX was not able to process the following string:
'QD$2\\\\uparrow,'

Here is the full report generated by LaTeX:
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=latex)
 restricted \write18 enabled.
entering extended mode
(./931d996306e720c622cc7e3c11201cf6.tex
LaTeX2e <2017-04-15>
Babel <3.18> and hyphenation patterns for 84 language(s) loaded.
(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
(/usr/share/texlive/texmf-dist/tex/latex/base/size10.clo))
(/usr/share/texlive/texmf-dist/tex/latex/type1cm/type1cm.sty)
(/usr/share/texlive/texmf-dist/tex/latex/base/textcomp.sty
(/usr/share/texlive/texmf-dist/tex/latex/base/ts1enc.def))
(/usr/share/texlive/texmf-dist/tex/latex/geometry/geometry.sty
(/usr/share/texlive/texmf-dist/tex/latex/graphics/keyval.sty)
(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifpdf.sty)
(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifvtex.sty)
(/usr/share/texlive/texmf-dist/tex/generic/ifxetex/ifxetex.sty)

Package geometry Warning: Over-specification in `h'-direction.
    `width' (5058.9pt) is ignored.


Package geometry Warning: Over-specification in `v'-direction.
    `height' (5058.9pt) is ignored.

)
No file 931d996306e720c622cc7e3c11201cf6.aux.
(/usr/share/texlive/texmf-dist/tex/latex/base/ts1cmr.fd)
*geometry* driver: auto-detecting
*geometry* detected driver: dvips
! Extra }, or forgotten $.
l.12 ...00000}{21.250000}{\rmfamily QD$2\uparrow,}
                                                  
! Missing $ inserted.
<inserted text> 
                $
l.13 \end{document}
                   
[1] (./931d996306e720c622cc7e3c11201cf6.aux) )
(\end occurred inside a group at level 1)

### simple group (level 1) entered at line 12 ({)
### bottom level
(see the transcript file for additional information)
Output written on 931d996306e720c622cc7e3c11201cf6.dvi (1 page, 336 bytes).
Transcript written on 931d996306e720c622cc7e3c11201cf6.log.
 

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    

def plotDOS(initials,line,limsx, limsy):
    Param_name = 't_1 = t_2'
    Param_name = '\epsilon_2'
    #Param_name = '\epsilon_2 = \\frac{-U_2}{2}'
    #Param_name = 't_2'
    
    f, axarr = subplots(2, 2, sharex='col', sharey='row', figsize=(15, 10))
    #plot the DOS at the FE
    t2 = []
    for ini in initials: 
        t2.append(readfile_Param(ini, line))
    t2 = array(t2) 
    
    #suptitle('Majorana Connected to the First QD')
    #print initials
    
    for i in range(4):
        vec0 = []
        maxi = []
        #for j in range(40,46):
        for j in initials:
            param = t2[where(initials == j)] ; param =  param[0]
            
            #Read file
            vec , vec2 = readfile(j, i)
            #plot DOS
            axarr[i/2 , i%2 ].plot(vec2,vec, '-',label= '$' + Param_name + '= '+str(param) + '$') 

        #Set graphic font
        if i%2 == 1: axarr[i/2 , i%2 ].legend(fontsize=16)
        #axarr[i/2 , i%2 ].set_ylim([0,max(maxi)])    
        #axarr[i/2 , i%2 ].set_xlim([-0.5,0.5])
        spin = '\uparrow' 
        if i%2 == 1 : spin = '\downarrow';  
        title = '$\\rho_{'+str(i/2 +1 )+ spin +'}$'
        axarr[i/2 , i%2 ].set_title(title)    
        axarr[i/2 , i%2 ].set_xlim([-limsx,limsx])
        axarr[i/2 , i%2 ].set_ylim([0,limsy])
        
        if i>1 : axarr[i/2 , i%2 ].set_xlabel('$\epsilon $',fontsize=26)
        if i%2 == 0 : axarr[i/2 , i%2 ].set_ylabel('$\\rho_{'+str(i/2 +1 )+ '}$',fontsize=26)
initials = array([451])      


def plotOneDot(i , dot):
    vec , w = readfile(i, dot)
    rho , vec2 = readfile(i+100, dot) ; vec = 0.5*(rho+vec)
    
    vec = vec#*pi*Gamma#*Hib
    spin = '\uparrow' if dot%2!=1 else '\downarrow'
            #plot DOS
    p0 = vec[w>0][0]
    print p0*pi*Gamma*2   
    #plot(w/(Gamma*0.8),vec*pi*Gamma/0.944036503513, 'o',label= '$'  +spin + '$') 
    plot(w/Gamma,vec, '-',label= '$'  +spin + '$')
#plotDOS(initials,17,0.1,4)            
plotOneDot(819, 0)
plotOneDot(819 , 1)
#plotOneDot(22 , 2)
#plotOneDot(22, 3)
#plot(linspace(-10*Gamma,10*Gamma,500),Green, label = 'Green')
#readMath('DQD')
#readMath('Gamma1=0.0282691,Gamma2=0,tdots=0.02,t1=0,t2=0.02')
xlim(-1.5,1.5)
#ylim(0,11)
legend(fontsize=17, loc='lower right')
#plotOneDot(556 , 0)
xlabel('$\omega/\Gamma$')
xticks((-10,0,10))
#yticks((0,5,10))
#ylabel('$\\rho$')


    

In [ ]:

    




    

In [ ]:

    

def plotOneDot(ax,i , dot):
    vec , w = readfile(i, dot)
    rho , vec2 = readfile(i+100, dot) ; vec = 0.5*(rho+vec)
    vec = vec#*pi*Gamma#*Hib
    spin = '\uparrow' if dot%2!=1 else '\downarrow'
    ln = 'b-' if dot%2!=1 else 'r--'
            #plot DOS
    p0 = vec[w>0][0]
    print p0*pi*Gamma*2   
    #plot(w/(Gamma*0.8),vec*pi*Gamma/0.944036503513, 'o',label= '$'  +spin + '$') 
    ax.plot(2*w/Gamma,vec, ln,label= '$'  +spin + '$')
i = 818
f, axarr = subplots(1, 2, sharex='col', sharey='row', figsize=(10, 5))
plotOneDot(axarr[0],i , 0)
plotOneDot(axarr[0],i , 1)
plotOneDot(axarr[1],i , 2)
plotOneDot(axarr[1],i , 3)
axarr[1].legend(fontsize=20, loc='upper right')
axarr[0].set_xlim(-10,10)
axarr[1].set_xlim(-10,10)
#axarr[0].set_xticks((-1.5,0,1.5))
#axarr[1].set_xticks((-1.5,0,1.5))

#axarr[0].set_ylim(0,11)
#axarr[0].set_yticks((0,1,2,3))

axarr[0].set_title('Dot 1')
axarr[1].set_title('Dot 2')

axarr[0].set_xlabel('$\omega/\Gamma$')
axarr[1].set_xlabel('$\omega/\Gamma$')
axarr[0].set_ylabel('DOS')
print readfile_Param(i, 12) , readfile_Param(i, 13)


    

In [ ]:

    

axarr[1].legend(fontsize=20, loc='upper right')


    

In [ ]:

    

def readMath(title):
    infile = open(os.path.abspath('Gamma1=0.0282691,Gamma2=0,tdots=0.02,t1=0,t2=0.02'), 'r')
    #infile = open(os.path.abspath('Gamma1=0.0282691,Gamma2=0,tdots=0.02,t1=0.02,t2=0.'), 'r')
    infile = open(os.path.abspath( title ), 'r')
    text = infile.readlines()
    Green = []
    for x in text:
         Green.append(pi*Gamma*float(x.split('\n')[0]))
    array(Green)
    #plot(linspace(-10*Gamma,10*Gamma,500),Green, label = 'Green')
    plot(linspace(-10,10,500),Green, label = 'Green')
#plot(Green)


    

In [ ]:

    

2.5*2**(0.5)
Lambda = 2.5
lambTerm = (2**(1/2))*(1+Lambda)**(-0.5)
lambTerm = (2/(1+Lambda**(-1)))**(0.5)
lambTerm
Hib = 0.5*(1+Lambda)/(-1+Lambda)*log(Lambda)


    

In [ ]:

    

2.5*Gamma*pi*(1.41)


    

In [ ]:

    

2**(0.5)


    

In [ ]:

    

from scipy.signal import argrelextrema
i = 451 ; dot = 1
def findLocalMax(i, dot):
    rho, w = readfile(i, dot) ;rho1, w1 = readfile(i+100, dot)
    rho = 0.5*(rho +rho1) ; w = 0.5*(w+w1)

    arg = array(argrelextrema(rho[1::2], greater))
    w1 = -array(w[arg])
    a = w1>0.0001 
    b = w1<0.5
    w1= w[arg[a&b]] ; rho1 = rho1[arg[a&b]]
    return w1,rho1


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    

xtitle = 't_1/D = t_2/D' ; xtitle = '\epsilon_{d2}/D';# xtitle = '\epsilon_{d2}/D = \\frac{-U_2}{2D}'; xtitle = 't_2/D'

#initials = array([8,10])

i = 451
rho, w = readfile(i, dot)
rho1, w1 = readfile(i+100, dot)
rho = 0.5*(rho +rho1) ; w = 0.5*(w+w1)
rho2, w2 = readfile(ini, 1) ; rho3, w = readfile(ini+1, 1) ; rho2 = rho2 +rho3
rho2 = 0.5*(rho2 +rho1)
Gamma = 2.82691*0.01
print pi , Gamma
plot( w , rho  ,color='b', marker="D" , linewidth=2.0)
#plot( w , rho2  ,color='g', marker="D" , linewidth=2.0)
#plot(w2/Gamma , rho2/maxi , label = '$QSz \\rho_\downarrow$',color='g', marker="D" ,linewidth=2.0)

#rho, w = readfile(10, 0)
#maxi = max(rho)
#rho2, w2 = readfile(10, 1)
#Gamma = 2.82691*0.01
#print pi , Gamma
#plot( w/Gamma , rho/maxi ,  label = 'NupPdn $\\rho_\uparrow$' ,color='b', marker="s" , linewidth=2.0)
#plot(w2/Gamma , rho2/maxi , label = 'NupPdn $\\rho_\downarrow$',color='g', marker="s" ,linewidth=2.0)

legend(fontsize=26)
xlim(-2,2)
#xlabel('$\epsilon/\Gamma$') ; ylabel('Density of States $(\pi \Gamma \\rho)$')
yticks(linspace(0,1,3), ('$0$', '$0.5$', '$1$'))
xticks(linspace(-2,2,5), ('$-2$','$-1$', '$0$','$1$', '$2$'))
#plotDOS(initials,16,1,4)
print maxi


    

In [ ]:

    




    

In [ ]:

    

def plotPH_Sym(initials , param , mini,maxi):
    f2, ax = subplots(size(initials), 1, sharex='col', sharey='row', figsize=(10, 4*(size(initials))))
    #suptitle('Majorana Connected to the First QD')
    for i in initials:
        num = where(initials == i)[0] ; num = num[0]
        for j in array([0]):
            dosUp , wUp = readfile(i, j)
            dosUp1 , wUp1 = readfile(i+100, j) 
            dosUp = dosUp1+ dosUp
            dosDown , wDown = readfile(i, j+1)
            dosDown1 , wSown1 = readfile(i+100, j+1) 
            dosDown = dosDown1+ dosDown
            #plot(wUp , dosUp, wUp , dosDown)
            vec1 = dosUp[ wUp > 0] + dosDown[ wDown > 0]
            vec2 = dosUp[ wUp < 0] + dosDown[ wDown < 0]
            #vec1 = dosUp[ wUp > 0] - dosUp[ wUp < 0][::-1] 
            #vec2 = dosDown[ wDown > 0][::-1] - dosDown[ wDown < 0]

            #print vec1 - vec2[::-1]
            #plot(wUp[ wUp > 0],vec1, wUp[ wUp < 0],vec2 )
            #ax.semilogx( wUp[ wUp > 0] , vec1/vec2[::-1] , 'o',label = 'D '+ str(j/2+1)+', i '+str(i))
            #if num%3 == 0 :
            #ax[num].semilogx( wUp[ wUp > 0] , vec2[::-1]/vec1, 'o',label = 'Dot '+str(j/2+1))
            #else :
            ax[num].loglog( wUp[ wUp > 0] , vec1/vec2[::-1], 'o',label = 'Dot '+str(j/2+1))
            #plot( range(size(vec1))[::-1] , vec1/vec2[::-1] , 'o',label = 'Dot '+ str(j/2 +1))
            #print vec1[-1], vec2[0], vec1[-1]/vec2[0] 
            #plot( log(wUp[ wUp > 0]), vec1 , 'o',label = ''+ str(j/2))
            #plot(  (-1)*log(wUp[ wUp > 0]) , vec2[::-1] , 'o',label = ''+ str(j/2))
            #xlim([-10,10])
            #ylim([1.0-0.00002,1.0+0.00002])
            #ax[num].set_title('$\Gamma =' + str(readfile_Param(i, 3))+', N_s = ' +repr(readfile_Param(i, 1))+', \epsilon_{d2} = ' +repr(readfile_Param(i, 16))+'$'  , fontsize=22)
            ax[num].set_title(' h_z = ' +repr(readfile_Param(i, param))+'$'  , fontsize=22)
        
        #ax[num].text(1, 1, r'$E=mc^2$', fontsize=15) 
        #ax[num].set_ylim(mini,maxi)
        ax[num].set_ylabel('$\\frac{\\rho_\uparrow(\omega^+)+ \\rho_\downarrow(\omega^+)}{\\rho_\uparrow(\omega^-) +\\rho_\downarrow(\omega^-)}$')
            #print vec1 - vec2[::-1]
    ax[num].legend(fontsize=20 , loc='upper left' )    
    ax[num].set_xlabel( '$\omega$' )  
                

        
def plotDOS(initials,line,limsx, limsy):
    Param_name = 't_1 = t_2'
    Param_name = '\epsilon_2'
    #Param_name = '\epsilon_2 = \\frac{-U_2}{2}'
    #Param_name = 't_2'
    
    f, axarr = subplots(2, 2, sharex='col', sharey='row', figsize=(15, 10))
    #plot the DOS at the FE
    t2 = []
    for ini in initials: 
        t2.append(readfile_Param(ini, line))
    t2 = array(t2) 
    
    #suptitle('Majorana Connected to the First QD')
    #print initials
    
    for i in range(4):
        vec0 = []
        maxi = []
        #for j in range(40,46):
        for j in initials:
            param = t2[where(initials == j)] ; param =  param[0]
            
            #Read file
            vec , vec2 = readfile(j, i)
            #plot DOS
            axarr[i/2 , i%2 ].plot(vec2,vec, '-',label= '$' + Param_name + '= '+str(param) + '$') 

        #Set graphic font
        if i%2 == 1: axarr[i/2 , i%2 ].legend(fontsize=16)
        #axarr[i/2 , i%2 ].set_ylim([0,max(maxi)])    
        #axarr[i/2 , i%2 ].set_xlim([-0.5,0.5])
        spin = '\uparrow' 
        if i%2 == 1 : spin = '\downarrow';  
        title = '$\\rho_{'+str(i/2 +1 )+ spin +'}$'
        axarr[i/2 , i%2 ].set_title(title)    
        axarr[i/2 , i%2 ].set_xlim([-limsx,limsx])
        axarr[i/2 , i%2 ].set_ylim([0,limsy])
        
        if i>1 : axarr[i/2 , i%2 ].set_xlabel('$\epsilon $',fontsize=26)
        if i%2 == 0 : axarr[i/2 , i%2 ].set_ylabel('$\\rho_{'+str(i/2 +1 )+ '}$',fontsize=26)




    
def logplotDOS(initials , line):
    #plot the 4 DOS
    f, axarr = subplots(2, 2, sharex='col', sharey='row', figsize=(15, 10))
    #plot the DOS at the FE
    f2, axarr2 = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 5))
    #plot the relation between both DOS (up/Down)
    f3, axarr3 = subplots(1, 1, sharex='col', sharey='row', figsize=(8, 5))
    
    t2 = []
    for j in initials: 
        t2.append(readfile_Param(j, line))
    t2 = array(t2) 
    for i in range(4):
        vec0 = []
        #for j in range(40,46):
        for j in initials:
            param = t2[where(initials == j)] ; param =  param[0]
            #Read file
            vec , vec2 = readfile(j, i)
            #get 0 value of the DOS
            if j%10 ==0 : vec = vec[::-1];
            vec0mas = vec[ vec2 > 0] ; vec0mas = vec0mas[0]
            vec0menos = vec[ vec2 < 0] ; vec0menos = vec0menos[-1]

            vec0.append(0.5*(vec0mas + vec0menos))

            #plot DOS
            #axarr[i/2 , i%2 ].plot(vec2,vec, label= '$\epsilon_{d2} ='+str(-0.25+(j-initials[0])*0.05 )+'$') 
            axarr[i/2 , i%2 ].semilogx(vec2[ vec2 > 0],vec[ vec2 > 0],'o',  label= '$'+ str()+'+$') 
            axarr[i/2 , i%2 ].semilogx(-vec2[ vec2 < 0],vec[ vec2 < 0],'o' ,label= '$'+ str(param)+',-$') 
        #Set graphic font
        if i%2 == 1: axarr[i/2 , i%2 ].legend(fontsize=12)
        #axarr[i/2 , i%2 ].set_xlim([-0.5,0.5])
        spin = '\uparrow' 
        if i%2 == 1 : spin = '\downarrow';  
        title = '$\\rho_{'+str(i/2 +1 )+ spin +'}$'
        print title
        axarr[i/2 , i%2 ].set_title(title)    
        #axarr[i/2 , i%2 ].set_ylim([-0.0,9.0])
        if i>1 : axarr[i/2 , i%2 ].set_xlabel('$\epsilon $',fontsize=25)
        if i%2 == 0 : axarr[i/2 , i%2 ].set_ylabel('$\\rho_{'+str(i/2 +1 )+ '}$',fontsize=25)

        #t2 = -0.25 + double(initials-initials[0]) * 0.05  
        #t2[-1] = -0.06 ; t2[-2] = -0.08 
        axarr2.plot(t2,vec0 ,'o', label= title + '$(0)$')


        if i%2 == 0 : 
            new0 = vec0
        else: 
            print array(vec0)/array(new0) 
            axarr3.plot(t2,array(new0)/array(vec0) ,'o', label= 'Dot' + str(i/2 + 1))

    axarr2.set_xlabel('$\epsilon_{d2}$',fontsize=25)
    axarr2.set_ylabel('$\\rho(0)$',fontsize=25)
    axarr2.legend(fontsize=24 , loc='upper right' )
    #axarr2.set_xlim([-0.25,-0.05])
    #axarr2.set_ylim([0.0,4.5])

    axarr3.set_xlabel('$\epsilon_{d2}$',fontsize=25)
    axarr3.set_ylabel('$\log_2(\\rho_\uparrow(0)/\\rho_\downarrow(0))$',fontsize=25)        
    axarr3.legend(fontsize=20, loc='upper right')        
    #axarr3.set_xlim([-0.25,-0.05])
    #axarr3.set_ylim([0.0,1.8])
    #save = directory +'/'+txt2 +repr(i)+'_'+repr(i)+'_OmegaRhow.png'
    #savefig(save)    

zvec2 = ['0.6','0.8','1','1.2','1.4']    
def readCost(ini, Ztrick):
        output = txt+repr(ini)+'/output_DMNRG_zEQ'+zvec2[Ztrick]+'_Dir'+repr(ini)+'.txt'
        #output = txt+repr(ini)+'/output_NRG_Code_DM_NRG_Dir'+repr(ini)+'.txt'
        infile = open(os.path.abspath(output), 'r')
        text = infile.readlines()
        w= []
        costi = []
        for x in text:
                #print(list(x))
            a=x.split('(')
            if a[0]=='RhoCosti0_0':
                omega = a[1].split('=')[1] ; omega = omega.split(' ')[1] 
                cost = a[1].split('=')[2] 
                w.append(float(omega)) ; costi.append(float(cost))
                #print float(omega) ,cost
                #Data managemente
        w = array(w) ; costi = array(costi)  
        costi = costi[argsort(w)] ;         w = w[argsort(w)]
        return w , costi