Revision 344d0dcb53a7ad0b30f0a5a2953b3e7f22630088 authored by uthpalah on 15 September 2019, 18:10:41 UTC, committed by uthpalah on 15 September 2019, 18:10:41 UTC
1 parent 5f15144
Raw File
scriptCompareBands.py
from .utilsprocar import UtilsProcar
from .procarparser import ProcarParser
from .procarselect import ProcarSelect
from .procarplotcompare import ProcarPlotCompare
import numpy as np
import matplotlib.pyplot as plt
import re




def bandscompare(file,file2,mode='plain',abinit_output=None,abinit_output2=None,spin=0,spin2=0,atoms=None,atoms2=None,orbitals=None,orbitals2=None,fermi=None,fermi2=None,elimit=None,mask=None,markersize=0.02,markersize2=0.02,cmap='jet',vmax=None,vmin=None,vmax2=None,vmin2=None,grid=True,marker=',',marker2=',',permissive=False,human=False,savefig=None,kticks=None,knames=None,title=None,outcar=None,outcar2=None,color='r',color2='g',legend='PROCAR1',legend2='PROCAR2',kpointsfile=None):
  #First handling the options, to get feedback to the user and check
  #that the input makes sense.
  #It is quite long
  if atoms is None:
    atoms = [-1]
    if human is True:
      print("WARNING: `--human` option given without atoms list!")
      print("--human will be set to False (ignored)\n ")
      human = False
  
  #repeat for 2nd data set
  if atoms2 is None:
    atoms2 = [-1]         
      
  if orbitals is None:
    orbitals = [-1]
  #repeat for 2nd data set  
  if orbitals2 is None:
    orbitals2 = [-1]  
    

  print("Script initiated")
  print("input file 1   : ", file)
  print("input file 2   : ", file2) #2nd file
  print("Mode           : ", mode)
  
  print("spin comp.  #1 : ", spin)
  print("spin comp.  #2 : ", spin2)
  print("atoms list. #1 : ", atoms)
  print("atoms list. #2 : ", atoms2)
  print("orbs. list. #1 : ", orbitals)
  print("orbs. list  #2 : ", orbitals2)

  if fermi is None and outcar is None and abinit_output is None:
    print("WARNING: Fermi Energy not set! ")
    print("You should use '-f' or '--outcar'\n Are you using Abinit Procar?\n")
    print("The zero of energy is arbitrary\n")
    fermi = 0
    
  if fermi2 is None and outcar2 is None and abinit_output2 is None:
    print("WARNING: Fermi Energy #2 not set! ")
    print("You should use '-f' or '--outcar'\n Are you using Abinit Procar?\n")
    print("The zero of energy is arbitrary\n")
    fermi2 = 0

  if kpointsfile is None:
    print("No KPOINTS file present. Please set knames and kticks manually.")    


###################reading abinit output (added by uthpala) ##########################

  if abinit_output:
  	print("Abinit output used")

  #reading abinit output file
  	rf = open(abinit_output,'r')
  	data = rf.read()
  	rf.close()

  	fermi = float(re.findall('Fermi\w*.\(\w*.HOMO\)\s*\w*\s*\(\w*\)\s*\=\s*([0-9.+-]*)',data)[0])
   
   
  if abinit_output2:
  	print("Abinit output #2 used")

  #reading abinit output file
  	rf2 = open(abinit_output2,'r')
  	data2 = rf2.read()
  	rf2.close()

  	fermi2 = float(re.findall('Fermi\w*.\(\w*.HOMO\)\s*\w*\s*\(\w*\)\s*\=\s*([0-9.+-]*)',data2)[0])


####################################################################  

 
  print("Fermi Ener. #1  : ", fermi)
  print("Fermi Ener. #2  : ", fermi2)
  print("Energy range    : ", elimit)

  if mask is not None:
    print("masking thres.: ", mask) 
    
  print("Colormap        : ", cmap)
  print("MarkerSize #1   : ", markersize)
  print("MarkerSize #2   : ", markersize2)
    
  print("Permissive      : ", permissive)
  if permissive:
    print("INFO: Permissive flag is on! Be careful")
  print("vmax            : ", vmax)
  print("vmin            : ", vmin)
  print("vmax #2         : ", vmax2)
  print("vmin #2         : ", vmin2)
  print("grid enabled    : ", grid) 
  if human is not None:
    print("human          : ", human)
  print("Savefig         : ", savefig)
  print("kticks          : ", kticks)
  if kpointsfile is None:
    print("knames          : ", knames)
    print("title           : ", title)

  print("outcar #1       : ", outcar)
  print("outcar #2       : ", outcar2)
  
  print("legend #1       : ",legend)
  print("legend #2       : ",legend2)

#If KPOINTS file is given:
  if kpointsfile is not None:
    #Getting the high symmetry point names from KPOINTS file
    f = open(kpointsfile)
    KPread = f.read()
    f.close()

    KPmatrix = re.findall('reciprocal[\s\S]*',KPread)
    tick_labels = np.array(re.findall('!\s(.*)',KPmatrix[0]))
    knames=[]
    knames=[tick_labels[0]]

    for i in range(len(tick_labels)-1):
      if tick_labels[i] !=tick_labels[i+1]:
        knames.append(tick_labels[i+1])

    knames = [str("$"+latx+"$") for latx in knames] 

    #getting the number of grid points from the KPOINTS file
    f2 = open(kpointsfile)
    KPreadlines = f2.readlines()
    f2.close()
    numgridpoints = int(KPreadlines[1].split()[0])

    kticks=[0]
    gridpoint=0
    for kt in range(len(knames)-1):
      gridpoint=gridpoint+numgridpoints
      kticks.append(gridpoint-1)
    print("knames        : ", knames)
    print("kticks        : ", kticks)  


  #If ticks and names are given by user manually:
  if kticks is not None and knames is not None:
    ticks = list(zip(kticks,knames))
  elif kticks is not None:
    ticks = list(zip(kticks,kticks))
  else:
    ticks = None
  
  #The spin argument should be a number (index of an array), or
  #'st'. In the last case it will be handled separately (later)

  #Commented out by Uthpala to preserve the spin=1 tagging convention
  #spin = {'0':0, '1':1, '2':2, '3':3, 'st':'st'}[spin]
  #spin2 = {'0':0, '1':1, '2':2, '3':3, 'st':'st'}[spin2]


  #The second part of this function is parse/select/use the data in
  #OUTCAR (if given) and PROCAR

  #first parse the outcar if given, to get Efermi and Reciprocal lattice
  recLat = None 
  recLat2 = None
  
  if outcar and outcar2:
    outcarparser = UtilsProcar()
    if fermi is None:
      fermi = outcarparser.FermiOutcar(outcar)
      print("INFO: Fermi energy found in outcar file = " + str(fermi))
    if fermi2 is None:  
      fermi2 = outcarparser.FermiOutcar(outcar2)
      print("INFO: Fermi energy #2 found in outcar file = " + str(fermi2))
      
      
    recLat = outcarparser.RecLatOutcar(outcar)
    recLat2 = outcarparser.RecLatOutcar(outcar2)

  # parsing the PROCAR file
  procarFile = ProcarParser()
  procarFile.readFile(file, permissive, recLat)
  
  # parsing the PROCAR file #2
  procarFile2 = ProcarParser()
  procarFile2.readFile(file2, permissive, recLat2)

  # processing the data, getting an instance of the class that reduces the data
  data = ProcarSelect(procarFile, deepCopy=True)
  data2 = ProcarSelect(procarFile2,deepCopy=True)
  
  #handling the spin, `spin='st'` is not straightforward, needs
  #to calculate the k vector and its normal. Other `spin` values
  #are trivial.
  if spin is 'st':
    #two `ProcarSelect` instances, to store temporal values: spin_x, spin_y
    dataX = ProcarSelect(procarFile, deepCopy=True)
    dataX.selectIspin([1])
    dataX.selectAtoms(atoms, fortran=human)
    dataX.selectOrbital(orbitals)  
    dataY = ProcarSelect(procarFile, deepCopy=True)
    dataY.selectIspin([2])
    dataY.selectAtoms(atoms, fortran=human)
    dataY.selectOrbital(orbitals)
    #getting the signed angle of each K-vector
    angle = np.arctan2(dataX.kpoints[:,1], (dataX.kpoints[:,0]+0.000000001))
    sin = np.sin(angle)
    cos = np.cos(angle)
    sin.shape = (sin.shape[0],1)
    cos.shape = (cos.shape[0],1)
    ##print sin, cos
    #storing the spin projection into the original array
    data.spd = -sin*dataX.spd + cos*dataY.spd
  else:
    data.selectIspin([spin])
    data.selectAtoms(atoms, fortran=human)
    data.selectOrbital(orbitals)
    
  #repeat for dataset #2

  if spin2 is 'st':
    #two `ProcarSelect` instances, to store temporal values: spin_x, spin_y
    dataX2 = ProcarSelect(procarFile2, deepCopy=True)
    dataX2.selectIspin([1])
    dataX2.selectAtoms(atoms2, fortran=human)
    dataX2.selectOrbital(orbitals2)  
    dataY2 = ProcarSelect(procarFile2, deepCopy=True)
    dataY2.selectIspin([2])
    dataY2.selectAtoms(atoms2, fortran=human)
    dataY2.selectOrbital(orbitals2)
    #getting the signed angle of each K-vector
    angle2 = np.arctan2(dataX2.kpoints[:,1], (dataX2.kpoints[:,0]+0.000000001))
    sin2 = np.sin(angle2)
    cos2 = np.cos(angle2)
    sin2.shape = (sin2.shape[0],1)
    cos2.shape = (cos2.shape[0],1)
    ##print sin, cos
    #storing the spin projection into the original array
    data2.spd = -sin2*dataX2.spd + cos2*dataY2.spd
  else:
    data2.selectIspin([spin2])
    data2.selectAtoms(atoms2, fortran=human)
    data2.selectOrbital(orbitals2)  
    
 
  # Plotting the data
  data.bands = (data.bands.transpose() - np.array(fermi)).transpose()    
  # Plotting the data for data #2
  data2.bands = (data2.bands.transpose() - np.array(fermi2)).transpose()
  plot = ProcarPlotCompare(data.bands, data2.bands, data.spd, data2.spd, data.kpoints, data2.kpoints)
  
  
  ###### start of mode dependent options #########

  if mode == "scatter":
    plot.scatterPlot(mask=mask,size=markersize,size2= markersize2, cmap=cmap, vmin=vmin, vmax=vmax,vmin2=vmin2, vmax2=vmax2, marker=marker, marker2=marker2,legend1=legend,legend2=legend2, ticks=ticks)
    if fermi is not None:
    	plt.ylabel(r"$E-E_f$ [eV]",fontsize=22)
    else:
    	plt.ylabel(r"Energy [eV]",fontsize=22)
    if elimit is not None:
      plt.ylim(elimit)
#
  if mode == "plain":
    plot.plotBands(size=markersize,size2= markersize2, marker=marker, marker2=marker2,color=color,color2=color2,legend1=legend,legend2=legend2, ticks=ticks)
    if fermi is not None:
    	plt.ylabel(r"$E-E_f$ [eV]",fontsize=22)
    else:
    	plt.ylabel(r"Energy [eV]",fontsize=22)
    if elimit:
      plt.ylim(elimit)
      
  if mode == "parametric":
    plot.parametricPlot(cmap=cmap, vmin=vmin, vmax=vmax,vmin2=vmin2, vmax2=vmax2, marker='solid', marker2='dashed', legend1=legend,legend2=legend2,ticks=ticks)
    if fermi is not None:
    	plt.ylabel(r"$E-E_f$ [eV]",fontsize=22)
    else:
    	plt.ylabel(r"Energy [eV]",fontsize=22)
    if elimit is not None:
      plt.ylim(elimit)

  elif mode == "atomic":
    plot.atomicPlot(cmap=cmap, vmin=vmin, vmax=vmax,vmin2=vmin2, vmax2=vmax2)
    if fermi is not None:
    	plt.ylabel(r"$E-E_f$ [eV]",fontsize=22)
    else:
    	plt.ylabel(r"Energy [eV]",fontsize=22)
    if elimit is not None:
      plt.ylim(elimit)

  ##### end of mode dependent options ###########
  plt.tight_layout()
  if grid:
    plt.grid()
  
  if title:
    plt.title(title,fontsize=22)

  if savefig:
    plt.savefig(savefig,bbox_inches='tight')
  else:
    plt.show()

  return
    
back to top