# Laboratorio de Fisica Molecular Aplicada (FAMA), UAM-A, Mexico # ======================================================================= # Dr. Oscar Olvera Neria # Script para sumar los orbitales de ciertos atomos seleccionados # # Cambio Usuario Descripcion # ========== =============== =========== # 27/05/2015 J.C Glez-Torres Script original # 10/05/2017 J.C Glez-Torres Adiciona la suma totalDcript original # 13/05/2017 J.C Glez-Torres Simplifica la entrada de archivos ######################################################################################### # use: # python sum_dos.py 1 n >> archivo.dat # 1 stands for the first atom to the n atom # # First you have to split the DOSCAR with the script # utility provided by Utexas split_dos. # WARNING: WE HAVE TO CHECK THAT the first line of the N DOSx # files do not contain a comment. # # This will generate the .dat file which contains the add states of the atoms. # Then you can plot this files. # ######################################################################################### from pylab import * #this includes numpy as np! from scipy.optimize import leastsq from math import log, exp import matplotlib.pyplot as plt import sys # initialize some variable to be lists energia = [] sup = [] sdown = [] pxup = [] pxdown = [] pyup = [] pydown = [] pzup = [] pzdown = [] d1up = [] d1down = [] d2up = [] d2down = [] d3up = [] d3down = [] d4up = [] d4down = [] d5up = [] d5down = [] msup = [],[] msdown = [],[] E = [],[] EM = [],[] fdos = [] totalmatrix = [],[] ttotalmatrix = [],[] tdosiup = [] tdosidown = [] tdospup = [] tdospdown = [] tdosdup = [] tdosddown = [] #Ciclo para captura de nombres de archivos #----------------------------------------------- a=sys.argv[1]; b=sys.argv[2] name='DOS' for i in range (int(a), int(b)+1): fdos.append(name + str(i)) rangox=(len(fdos)) #print(rangox) #print(fdos) #Captura los rangos de energia de la DOS #------------------------------------------------- f2=open(fdos[0]) lines= f2.readlines() f2.close() for line in lines: p = line.split() energia.append(str(p[0])) E = np.vstack(energia) rangoy=(len(E)) #Numero de filas totales que tiene el archivo #------------------------------------------------- # Termina captura de archivos #Captura de datos por columna #------------------------------------------------- for i in range(len(fdos)): # print (fdos[i]) f2 = open(fdos[i]) lines = f2.readlines() f2.close() for line in lines: p = line.split() #Orbitales S sup.append(float(p[1])) sdown.append(float(p[2])) #Orbitales P pxup.append(float(p[3])) pxdown.append(float(p[4])) pyup.append(float(p[5])) pydown.append(float(p[6])) pzup.append(float(p[7])) pzdown.append(float(p[8])) #Orbitales D d1up.append(float(p[9])) d1down.append(float(p[10])) d2up.append(float(p[11])) d2down.append(float(p[12])) d3up.append(float(p[13])) d3down.append(float(p[14])) d4up.append(float(p[15])) d4down.append(float(p[16])) d5up.append(float(p[17])) d5down.append(float(p[18])) #Se almacenan los vectores de orbitales en #Matrices de orbitales #Matriz de orbitales S msup = np.vstack(sup) msdown = np.vstack(sdown) #Matriz de orbitales P mpxup = np.vstack(pxup) mpxdown = np.vstack(pxdown) mpyup = np.vstack(pyup) mpydown = np.vstack(pydown) mpzup = np.vstack(pzup) mpzdown = np.vstack(pzdown) #Matriz de orbitales D md1up = np.vstack(d1up) md1down = np.vstack(d1down) md2up = np.vstack(d2up) md2down = np.vstack(d2down) md3up = np.vstack(d3up) md3down = np.vstack(d3down) md4up = np.vstack(d4up) md4down = np.vstack(d4down) md5up = np.vstack(d5up) md5down = np.vstack(d5down) #------------------------------------------------- #Convierto los vectores en matrices #------------------------------------------------- shape = ( rangox, rangoy ) #orbitales S msup.shape = shape msdown.shape = shape #orbitales P mpxup.shape = shape mpxdown.shape = shape mpyup.shape = shape mpydown.shape = shape mpzup.shape = shape mpzdown.shape = shape #orbitales D md1up.shape = shape md1down.shape = shape md2up.shape = shape md2down.shape = shape md3up.shape = shape md3down.shape = shape md4up.shape = shape md4down.shape = shape md5up.shape = shape md5down.shape = shape #Transpuesta de las matrices para despues sumar sobre los renglores #------------------------------------------------------------------ #Orbitales S tsup = msup.transpose() tsdown = msdown.transpose() #Orbitales P tpxup = mpxup.transpose() tpxdown = mpxdown.transpose() tpyup = mpyup.transpose() tpydown = mpydown.transpose() tpzup = mpzup.transpose() tpzdown = mpzdown.transpose() #Orbitales D td1up = md1up.transpose() td1down = md1down.transpose() td2up = md2up.transpose() td2down = md2down.transpose() td3up = md1up.transpose() td3down = md1down.transpose() td4up = md1up.transpose() td4down = md1down.transpose() td5up = md1up.transpose() td5down = md1down.transpose() #Suma todos los renglones de cada matriz y obtengo el total #------------------------------------------------------------ #Orbitales S totalsup = tsup.sum(axis=1) totalsdown = tsdown.sum(axis=1) #Orbitales P totalpxup = tpxup.sum(axis=1) totalpyup = tpyup.sum(axis=1) totalpzup = tpzup.sum(axis=1) totalpxdown = tpxdown.sum(axis=1) totalpydown = tpydown.sum(axis=1) totalpzdown = tpzdown.sum(axis=1) #Orbitales D totald1up = td1up.sum(axis=1) totald2up = td2up.sum(axis=1) totald3up = td3up.sum(axis=1) totald4up = td4up.sum(axis=1) totald5up = td5up.sum(axis=1) totald1down = td1down.sum(axis=1) totald2down = td2down.sum(axis=1) totald3down = td3down.sum(axis=1) totald4down = td4down.sum(axis=1) totald5down = td5down.sum(axis=1) #A continuacion inserto una serie de controles #Para validar datos #----------------------------------------------- #totalmatrix = np.append(energia) #print (totalmatrix), "\t" #totalmatrix = np.append(totalsup) #print (totalmatrix), "\t" #totalmatrix = np.append(totalsdown) #print (totalmatrix), "\t" #print "Matriz S spin up\t" #print (msup) #print () #print () #print "Matriz S spin down\t" #print (msdown) #print () #print () #print "Matriz transpuesta S spin up\t" #print (tsup) #print () #print () #print "Matriz transpuesta S spin down\t" #print (tsdown) #print () #print () #print "Matriz transpuesta D4 spin up\t" #print (td4up) #print "Matriz transpuesta D4 spin down\t" #print (td4down) #print () #print () #print () #print () #print (totalsup) #print () #print () #print (totalsdown) #print () #print () #print (totalmatrix) #np.round(totalsup) #TERMINA LA PARTE DE CONTROL #------------------------------------------------------------------- #Impresion de los vectores de todos los orbitales con spin up y down #--------------------------------------------------------------------- np.set_printoptions(precision=6) #print ("Energia Sup total Sdown total Px up Px down Py up Py down Pz up Pz down") #for i in range(len(energia)): # x0=energia[i] # x1=totalsup[i] # x2=totalsdown[i] # x3=totalpxup[i] # x4=totalpxdown[i] # x5=totalpyup[i] # x6=totalpydown[i] # x7=totalpzup[i] # x8=totalpzdown[i] # x9=totald1up[i] # x10=totald1down[i] # x11=totald2up[i] # x12=totald2down[i] # x13=totald3up[i] # x14=totald3down[i] # x15=totald4up[i] # x16=totald4down[i] # x17=totald5up[i] # x18=totald5down[i] # print '%4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f' % (float(x0), float(x1), float(x2), float(x3), float(x4), float(x5), float(x6), float(x7), float(x8), float(x9),float(x10), float(x11), float(x12), float(x13), float(x14), float(x15), float(x16), float(x17), float(x18) ) #OBTENCION DE LOS ORBITALES TOTALES ES DECIR LA SUMA DE TODAS #SUS DIRECCIONES #--------------------------------------------------------------- tdospup = totalpxup + totalpyup + totalpzup tdospdown = totalpxdown + totalpydown + totalpzdown tdosdup= 0 tdosddown=0 tdosdup = totald1up + totald2up + totald3up + totald4up + totald5up tdosddown = totald1down + totald2down + totald3down + totald4down + totald5down tdosiup = totalsup + tdospup + tdosdup tdosidown = totalsdown + tdospdown + tdosddown print ("#Energia S up S down P up P down D up D down Total up Total down ") for i in range(len(energia)): x0=energia[i] x1=totalsup[i] x2=totalsdown[i] x3=tdospup[i] x4=tdospdown[i] x5=tdosdup[i] x6=tdosddown[i] x7=tdosiup[i] x8=tdosidown[i] print '%4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f %4.8f ' % ( float(x0), float(x1), float(x2), float(x3), float(x4), float(x5), float(x6), float(x7), float(x8)) #Graficamos los datos generados #-------------------------------------------------------------------- plt.fill(energia, tdosiup, label='Total (up)', color='#C5C1C7') plt.fill(energia, tdosidown, label='Total (down)', color='#C5C1C7') plt.plot(energia , tdospup, label='P (up)') plt.plot(energia , tdospdown, label = 'P (down)') plt.plot(energia , tdosdup, label='D (up)') plt.plot(energia , tdosddown, label='D (down)') plt.legend(loc='upper left') show()