#!/usr/bin/env python #ReadDOSCAR: Sung Sakong, Institute of Theoretical Chemistry, Ulm University def ReadDOSCAR(filename,eref): #input #filename: path_filename of DOSCAR #eref: reference energy, it will be the Fermi level for False import numpy as np try: with open(filename,'r') as f: lines=f.readlines() f.close() except IOError: print('DOSCAR is not accessable.\n') return False, False tmp=lines[0].split() nions=np.int(tmp[0]) ndata=np.int(tmp[1]) jobpar=np.int(tmp[2]) kblockr=np.int(tmp[3]) tmp=lines[5].split() emax=np.float(tmp[0]) emin=np.float(tmp[1]) ndos=np.int(tmp[2]) if eref==False: efermi=np.float(tmp[3]) else: efermi=eref dim=len(lines[6].split()) tdos=np.zeros((dim,ndos)) for i in range(6,ndos+6,1): tmp=lines[i].split() for j in range(dim): tdos[j,i-6]=np.float(tmp[j]) tdos[0,i-6]=tdos[0,i-6]-efermi if jobpar==0: pdos=False tdos[0,:]=tdos[0,:]-efermi if dim==5: Spin='T' elif dim==3: Spin='F' print('Atoms=',nions,'NEDos=',ndos,'PDos= F Spin=',Spin,'\n') return tdos, pdos elif jobpar>0: dim=len(lines[ndos+7].split()) pdos=np.zeros((nions,dim,ndos)) for k in range(nions): start=6+(k+1)*(ndos+1) for i in range(start,start+ndos,1): tmp=lines[i].split() for j in range(dim): pdos[k,j,i-start]=np.float(tmp[j]) pdos[k,0,i-start]=pdos[k,0,i-start]-efermi if dim==19: Spin='T' lm='T' elif dim==10: Spin='F' lm='T' elif dim==4: Spin='F' lm='F' elif dim==7: Spin='T' lm='F' elif dim==13: Spin='F' lm='F' print('Atoms=',nions,'NEDos=',ndos,'PDos= T Spin=',Spin,'lm=',lm,'\n') return tdos, pdos else: print('DOSCAR is contaminated.\n') return tdos, False #Data structure of the returned arrays: tdos=total dos, pdos=projected dos #Spin=F: tdos[0]=E-Eref,tdos[1]=tdos,tdos[2]=integrated tdos #Spin=T: tdos[0]=E-Eref,tdos[1]=tdos_up,tdos[2]=tdos_dn,tdos[3]=integrated tdos_up,tdos[3]=integrated tdos_dn # #PDos=T: #pdos[atom,orbital] #atom=atom index 0,...,nions-1 #Spin=F,lm=F: #pdos[atom,0]=E-Eref,pdos[atom,1]=s,pdos[atom,2]=p,pdos[atom,3]=d #Spin=T,lm=F: #pdos[atom,0]=E-Eref,pdos[atom,1]=s_up,pdos[atom,2]=s_dn, # pdos[atom,3]=p_up,pdos[atom,4]=p_dn,pdos[atom,5]=d_up,pdos[atom,6]=d_dn #Spin=F,lm=T: #pdos[atom,0]=E-Eref,pdos[atom,1]=s,pdos[atom,2]=py,pdos[atom,3]=pz,pdos[atom,4]=px, # pdos[atom,5]=dxy,pdos[atom,6]=dyz,pdos[atom,7]=dz2-r2,pdos[atom,8]=dxz,pdos[atom,9]=dx2-y2 #Spin=T,lm=T: #pdos[atom,0]=E-Eref,pdos[atom,1]=s_up,pdos[atom,2]=s_dn,... #Usage example: # import matplotlib.pyplot as plt plt.rcParams.update({'font.family': 'arial'}) plt.rcParams.update({'font.size': 16}) plt.figure(figsize=(8,6)) #plt.xlim(-25,7.5) #plt.ylim(-0.2,0.2) plt.xlabel(r'Energy $E-E_\mathrm{F}$ (eV)') plt.ylabel(r'Local density of states (Arbitrary Unit)') plt.yticks([]) tdos,pdos=ReadDOSCAR('DOSCAR', False) plt.plot(tdos[0],tdos[1]) plt.plot(tdos[0],-tdos[2]) #example of summed p-dos of atom 100 for spin up and down states of lm decomposed DOS #atom=99 #plt.plot(tdos[0],pdos[atom,3]+pdos[atom,5]+pdos[atom,7],c='blue') #plt.plot(tdos[0],-pdos[atom,4]-pdos[atom,6]-pdos[atom,8],c='red') plt.show() #plt.savefig('Fig_DOS.pdf', bbox_inches='tight', transparent=True)