#DOS Plotter #DOS_PLOT_PARAMETER file format: # #Cd Te elements #4 4 num atoms for each element #s p d orbitals to plot for element 1, in order spdf (can skip orbitals, see below). #s p d orbitals to plot for element 2 etc #-7 7 x axis range #-1 1 y axis range # #if no orbitals required for element use an x on element's orbital line, eg: #x #s d #will plot s and d orbitals of element 2 only. import matplotlib.pyplot as plt params = {'legend.fontsize': 'x-large', #sets text sizes 'axes.labelsize': 'x-large', 'axes.titlesize':'x-large', 'xtick.labelsize':'x-large', 'ytick.labelsize':'x-large'} plt.rcParams.update(params) input_file=open("DOSCAR","r") p_file=open("DOS_PLOT_PARAMETERS.txt") plottable=open("DOS_DATA.txt","w") parameters=p_file.readlines() #read parameters file elements=parameters[0].split() num_el=len(elements) #find num elements num_els=parameters[1].split() for i in range(len(parameters)): dummy=parameters[i].split() if dummy[0]=='Total': plot_total=True break else: plot_total=False colours=['C0','C1','C2','C3','C4','C5','C6','C7','C8','C9'] # specify own colour list to keep up+down same orbitals=[] cum_sum=[] #build cummulative list to find new element print points pp=0 for element in range(num_el): pp=pp+int(num_els[element]) cum_sum.append(pp) total_atoms=0 for i in range(0, len(num_els)): #find total number atoms total_atoms=total_atoms+int(num_els[i]) line=input_file.readlines() # read in whole file line by line dummy=line[5].split() # split parameter line num_states=int(dummy[2]) # find number of states in each atom block fermi=float(dummy[3]) # find fermi level plottable.write(str(num_el)+' '+str(num_states)+'\n') energy=[] #define empty lists for data su=[0.0]*num_states sd=[0.0]*num_states pu=[0.0]*num_states pd=[0.0]*num_states du=[0.0]*num_states dd=[0.0]*num_states tu=[0.0]*num_states td=[0.0]*num_states totalup=[0.0]*num_states totaldown=[0.0]*num_states leg=[] x=0 for l in range(6, 6+num_states): #append to energy list only once! dummy=line[l].split() energy.append(float(dummy[0])-fermi) x=x+1 start_line=7+num_states #start and end points lines of first atom block defined end_line=start_line+num_states y=0 #variable for print points orbital_plot=2 c=0 b=0 for atom in range(0,total_atoms): #loop through all atoms x=0 for l in range(start_line,end_line): #loop through each block dummy=line[l].split() su[x]=su[x]+(float(dummy[1])/float(num_els[b])) #spd orbitals added to lists sd[x]=sd[x]-(float(dummy[2])/float(num_els[b])) pu[x]=pu[x]+(float(dummy[3])/float(num_els[b]))+(float(dummy[5])/float(num_els[b]))+(float(dummy[7])/float(num_els[b])) pd[x]=pd[x]-(float(dummy[4])/float(num_els[b]))-(float(dummy[6])/float(num_els[b]))-(float(dummy[8])/float(num_els[b])) du[x]=du[x]+(float(dummy[9])/float(num_els[b]))+(float(dummy[11])/float(num_els[b]))+(float(dummy[13])/float(num_els[b]))+(float(dummy[15])/float(num_els[b]))+(float((dummy[17]))/float(num_els[b])) dd[x]=dd[x]-(float(dummy[10])/float(num_els[b]))-(float(dummy[12])/float(num_els[b]))-(float(dummy[14])/float(num_els[b]))-(float(dummy[16])/float(num_els[b]))-(float(dummy[18])/float(num_els[b])) tu[x]=su[x]+pu[x]+du[x] td[x]=sd[x]+pd[x]+dd[x] x=x+1 start_line=start_line+num_states+1 #move to next atom block end_line=end_line+num_states+1 if atom==cum_sum[y]-1: #write data to plottable when element changes > print point found plottable.write(elements[y]+"\n") b=b+1 for state in range(0,num_states): plottable.write(str(energy[state])+' '+str(su[state])+' '+str(sd[state])+' '+str(pu[state])+' '+str(pd[state])+' '+str(du[state])+' '+str(dd[state])+' '+str(tu[state])+' '+str(td[state])+"\n") orbitals=parameters[orbital_plot].split() # read parameters file for spd orbitals, only plot when specified for o in range(0,len(orbitals)): if orbitals[o]=='s': plt.plot(energy,su,color=colours[c], label=elements[y]+' s') #Plot s up using color c, label in legend plt.plot(energy,sd,color=colours[c], label=None) #plot s down using same colour as up, no legend label c=c+1 #next colour if c==10: #loop around colour list, modulo better? c=c-10 #how to do modulus? elif orbitals[o]=='p': plt.plot(energy,pu,color=colours[c],label=elements[y]+' p') #repeat above for p and d plt.plot(energy,pd,color=colours[c],label=None) c=c+1 if c==10: c=c-10 elif orbitals[o]=='d': plt.plot(energy,du,color=colours[c],label=elements[y]+' d') plt.plot(energy,dd,color=colours[c],label=None) c=c+1 if c==10: c=c-10 elif orbitals[o]=='t': plt.plot(energy,tu,color=colours[c],label=elements[y]) plt.plot(energy,td,color=colours[c],label=None) c=c+1 if c==10: c=c-10 orbital_plot=orbital_plot+1 #next atom spd line y=y+1 #next element name for e in range(0,len(tu)): totalup[e]=totalup[e]+tu[e] totaldown[e]=totaldown[e]+td[e] su=[0.0]*num_states #reset lists back to 0 for next element sd=[0.0]*num_states pu=[0.0]*num_states pd=[0.0]*num_states du=[0.0]*num_states dd=[0.0]*num_states tu=[0.0]*num_states td=[0.0]*num_states #plottable.write('END OF THE DATA') #write to end data, needs 4 words to match above line search plottable.write('Total\n') for i in range(num_states): plottable.write(str(energy[i])+' '+str(totalup[i])+' '+str(totaldown[i])+'\n') if plot_total==True: plt.plot(energy,totalup,color=colours[c],label='Total') plt.plot(energy,totaldown,color=colours[c],label=None) c=c+1 if c==10: c=c-10 input_file.close() plottable.close() axis=parameters[2+num_el].split() #read axis range from parameter file and set for plot plt.xlim(float(axis[0]),float(axis[1])) axis=parameters[2+num_el+1].split() #read axis range from parameter file and set for plot plt.ylim(float(axis[0]),float(axis[1])) plt.xlabel('E-E(Fermi) (eV)') #axis label axis=parameters[2+num_el+1].split() 0 plt.ylabel('Normalised DOS (States/eV)') p_file.close() plt.legend(loc='upper right') #show legend plt.savefig('DOS.png') #show plot