#!/usr/bin/env python # -*- coding: utf-8 -*- """ Abinit Post Process Application author: Martin Alexandre last edited: May 2013 """ import os, sys, time import string, math, re import commands from scipy.io import netcdf import numpy as np #Utility import utility.analysis as Analysis try: from PyQt4 import QtGui,QtCore except : pass; class MolecularDynamicFile: Ha_eV = 27.2113845e0 # Conversion factor 1 Hartree = 27,2113845 eV kb_eVK = 8.617342e-5 # Boltzman constant in eV/K = (1.380658e-23)/(1.60217733e-19) #eV/K amu_emass = 1.660538782e-27/9.10938215e-31 # 1 atomic mass unit, in electronic mass fact = 29421.033e0 # Conversion factor hartree/bohr^3 en GPa namefile1 = "" namefile2 = "" ni = 0 nf = 0 Nbtime = 0 goodFile= False type_of_file = "" #-----------------------------# #-------CONSTRUCTOR-----------# #-----------------------------# def __init__(self, pnamefile): # Take the last HIST/OUT.nc file in the directory (when APPA is launch) if pnamefile == "": OUT_list = [] filesdir = os.listdir(os.getcwd()) for files in filesdir : if os.path.isfile(files): stats = os.stat(files) fic_tuple = time.localtime(stats[8]), files lili = string.split(files, '.') # if files.find('_OUT.nc') == len(files)-7 and len(files)>7 : # OUT_list.append(fic_tuple) if files.find('HIST.nc') == len(files)-7 and len(files)>7 : OUT_list.append(fic_tuple) if files.find('HIST.nc') == len(files)-7 and len(files)>7 : OUT_list.append(fic_tuple) if len(OUT_list) > 0 : OUT_list.sort() ; OUT_list.reverse() fic_HIST=OUT_list[0][1].replace('_OUT.nc','_HIST') fic_HIST=OUT_list[0][1] if os.path.exists(fic_HIST) : self.namefile1 = fic_HIST self.namefile2 = fic_HIST # self.namefile2 = fic_HIST.replace('_HIST','_OUT.nc') self.type_of_file = 'netcdf' # or the file in parameter : else: if str(pnamefile).find('_OUT.nc') != -1: self.namefile2 = str(pnamefile) self.namefile1 = str(pnamefile).replace('_OUT.nc','_HIST') self.type_of_file = 'netcdf' elif str(pnamefile).find('HIST.nc') != -1 : self.namefile1 = str(pnamefile) if(os.path.exists(self.namefile1.replace('HIST.nc','OUT.nc'))): self.namefile2 = str(pnamefile).replace('HIST.nc','OUT.nc') else: self.namefile2 = str(pnamefile) self.type_of_file = 'netcdf' self.namefile1 = str(pnamefile) self.type_of_file = 'netcdf' elif str(pnamefile).find('HIST') != -1 : self.namefile1 = str(pnamefile) if(os.path.exists(self.namefile1.replace('HIST','OUT.nc'))): self.namefile2 = str(pnamefile).replace('HIST','OUT.nc') else: self.namefile2 = str(pnamefile) self.type_of_file = 'netcdf' else: self.type_of_file = 'ASCII_output' self.namefile2 = str(pnamefile) if (self.type_of_file == 'ASCII_output'): self.read_ascii() if (self.type_of_file == 'netcdf'): self.read_netcdf() #-----------------------------# #--------METHODS--------------# #-----------------------------# def read_ascii(self): if 1==1: #try: #Read the number of image: message = "Unable to read the number of image " try: temp = commands.getoutput('grep \"nimage =\" ' + self.namefile2 + ' | awk \'{print $3}\'') self.n_image = int(temp) except: self.n_image = 1 #Reading potential energy: message = "Unable to read potential energy " if self.n_image == 1 : temp = commands.getoutput('grep \"Total energy\" ' + self.namefile2 + ' | awk \'{print $5}\'') else: temp = commands.getoutput('grep \"Potential energy\" ' + self.namefile2 + ' | awk \'{print $4}\'') self.E_pot = np.array(temp.split('\n'), dtype=float) self.ni = 1 #Set the initial step to 1 temp = len(self.E_pot) self.nf = temp self.Nbtime = temp #Reading the number of atom message = "Unable to read the number of atoms" temp = commands.getoutput('grep \"natom\" ' + self.namefile2 + ' | awk \'{print $6}\'') self.natom = int(temp) #Reading acell: message = "Unable to read acell" temp = commands.getoutput('grep \" acell \" '+self.namefile2+' | awk \'{print $2,$3,$4 }\' ') if len(np.array(temp.split(),dtype=float)) == 3: acell = np.reshape(np.array(temp.split(),dtype=float),3) self.acell = np.array([[acell,]*self.n_image,]*self.Nbtime) #Reading the primitive vectors: message = "Unable to read the primitive vectors" temp = commands.getoutput(' grep -A 3 \"Real(R)+Recip(G)\" '\ + self.namefile2 + ' | awk \'{print $2,$3,$4 }\' | sed \'/space primitive vectors,/d\' ') if len(np.array(temp.split(),dtype=float)) == 9: rprim = np.reshape(np.array(temp.split(),dtype=float),(3,3)) self.rprimd = np.array([[rprim,]*self.Nbtime]) else: self.rprimd = np.reshape(np.array(temp.split()[0:self.n_image*self.Nbtime*9],dtype=float),(self.Nbtime,self.n_image,3,3)) self.rprimd[:,:,0,:] = self.rprimd[:,:,0,:] * self.acell[0,:] self.rprimd[:,:,1,:] = self.rprimd[:,:,1,:] * self.acell[1,:] self.rprimd[:,:,2,:] = self.rprimd[:,:,2,:] * self.acell[2,:] #Reading velocity of particules: message = "Unable to read the velocity of particules" temp = commands.getoutput(' grep -A '+str(self.natom)+' \"Cartesian velocities\" '\ + self.namefile2 + ' | awk \'{print $1,$2,$3 }\'| sed \'/--/d\' | sed \'/Cartesian velocities (vel)/d\' ' ) self.vel = np.reshape(np.array(temp.split()[0:self.n_image*self.Nbtime*self.natom*3],dtype=float),(self.Nbtime,self.n_image,self.natom,3)) #Reading position of particules (Cartesian): message = "Unable to read the position of particules (cart)" temp = commands.getoutput(' grep -A '+str(self.natom)+' \"Cartesian coordinates\" '\ + self.namefile2 + ' | awk \'{print $1,$2,$3 }\'| sed \'/--/d\' | sed \'/Cartesian coordinates/d\' ' ) self.xcart = np.reshape(np.array(temp.split()[0:self.n_image*self.Nbtime*self.natom*3],dtype=float),(self.Nbtime,self.n_image,self.natom,3)) #Reading position of particules (Reduced): message = "Unable to read the position of particules (red)" temp = commands.getoutput(' grep -A '+str(self.natom)+' \"Reduced coordinates\" '\ + self.namefile2 + ' | awk \'{print $1,$2,$3 }\'| sed \'/--/d\' | sed \'/Reduced coordinates/d\' ' ) self.xred = np.reshape(np.array(temp.split()[0:self.n_image*self.Nbtime*self.natom*3],dtype=float),(self.Nbtime,self.n_image,self.natom,3)) #Reading stress of particules: message = "Unable to read stress of particules" if self.n_image==1: temp = commands.getoutput(' grep -A 3 \"Cartesian components of stress\" '+ self.namefile2\ + ' | awk \'{print $3,$6 }\'| sed \'/^\s/d\' | sed \'/of/d\' ') else: temp = commands.getoutput(' grep -A 3 \"Cartesian components of stress\" '+ self.namefile2\ + ' | awk \'{print $4,$7 }\'| sed \'/^\s/d\' | sed \'/of/d\'| sed \'/Pressure/d\'') self.temp = np.reshape(np.array(temp.split()[0:self.n_image*self.Nbtime*6],dtype=float),(self.Nbtime,self.n_image,6)) s = (self.Nbtime,self.n_image,6) self.stress = np.zeros(s) nb = len(temp) / 6 ni = self.Nbtime - nb self.stress[ni:] = self.temp self.stress[ni:,:,1]= self.temp[:,:,2] self.stress[ni:,:,2]= self.temp[:,:,4] self.stress[ni:,:,3]= self.temp[:,:,1] self.stress[ni:,:,4]= self.temp[:,:,3] del self.temp #Reading the type of particules: message = "Unable to read the type of particules" temp = commands.getoutput('grep -A '+str(int(self.natom/20))+ ' \" typat \" '+self.namefile2+' | sed \'s/typat//\'') self.typat = np.array(temp.split(), dtype=int) #Reading mass of particules: temp = commands.getoutput('grep -P -A' + str(int(max(self.typat))/3) +' \'amu\s(.*)\' '+self.namefile2+' | sed \'s/amu//\'') self.amu = np.array(temp.split(), dtype=float) #Reading znucl : message = "Unable to read znucl" temp = commands.getoutput('grep \"znucl \" '+self.namefile2+' | sed \'s/znucl//\'') self.znucl = np.array(temp.split(), dtype=float) #Reading dtion: try: #try to read dtion : temp = commands.getoutput('grep \"dtion \" '+self.namefile2+' | awk \'{print $2}\'') self.dtion = float(temp) except: #If dtion doesn't exist, it's set to 100 (default value) self.dtion = 100 self.goodFile = True #-----------Calculation of some quantities------------# #-----------not available in the output File----------# self.mass_calculation() self.volume_calculation() self.angles_calculation() self.E_kinCalculation() self.temperature_calculation() self.pressure_calculation() #-----------------------------------------------------# #-----------------------------------------------------# # except: # self.goodFile = False # print "Error, " +message+ ", please use netcdf file " # try: # QtGui.QMessageBox.critical(None,"Warning","Error, " +message+ ", please use netcdf file ")# # except : # pass; def read_netcdf(self): #Check the version of scipy for the netcdf library module = __import__('scipy') if float(module.__version__.split('.')[1]) >= 9: module = "data" else: module = "__array_data__" #if 1==1: try : #Set the two netcdf_file variable with the HOST and ncfile: self.hist_file = netcdf.netcdf_file(str(self.namefile1), 'r') self.nc_file = netcdf.netcdf_file(str(self.namefile2), 'r') #---------------READING DATA----------------------# self.ni = 1 #Set the initial step to 1 temp = getattr(self.hist_file.variables['etotal'],module).shape[0] self.nf = temp self.Nbtime = temp #Set the number of image (only 1 for netcdf) self.n_image = 1 #Reading the number of atom self.natom = self.hist_file.dimensions['natom'] #Reading the primitive vectors: self.rprimd = np.reshape(getattr(self.hist_file.variables['rprimd'],module),(self.Nbtime,self.n_image,3,3)) #Reading acell: self.acell = np.reshape(getattr(self.hist_file.variables['acell'],module),(self.Nbtime,self.n_image,3)) #Reading stress of particules: self.stress = np.reshape(getattr(self.hist_file.variables['strten'],module),(self.Nbtime,self.n_image,6)) #Reading the type of particules: self.typat = getattr(self.nc_file.variables['typat'],module) #Reading mass of particules: self.amu = getattr(self.nc_file.variables['amu'],module) #Reading potential energy: self.E_pot = getattr(self.hist_file.variables['etotal'],module) #Reading potential energy: self.E_kin_ion = getattr(self.hist_file.variables['ekin'],module) #Reading potential energy: self.mdtime = getattr(self.hist_file.variables['mdtime'],module) #Reading position of particules: self.xcart = np.reshape(getattr(self.hist_file.variables['xcart'],module),(self.Nbtime,self.n_image,self.natom,3)) self.xred = np.reshape(getattr(self.hist_file.variables['xred'] ,module),(self.Nbtime,self.n_image,self.natom,3)) #Reading forces of particules: self.fcart = np.reshape(getattr(self.hist_file.variables['fcart'],module),(self.Nbtime,self.n_image,self.natom,3)) self.fred = np.reshape(getattr(self.hist_file.variables['fred'] ,module),(self.Nbtime,self.n_image,self.natom,3)) #Reading velocity of particules: self.vel = np.reshape(getattr(self.hist_file.variables['vel'],module),(self.Nbtime,self.n_image,self.natom,3)) #Reading dtion: try: #try to read dtion : self.dtion = getattr(self.nc_file.variables['dtion'],module) except: # If dtion doesn't exist, it's set to default abinit value => 100 : self.dtion = 100 #Reading znucl of particules: self.znucl = getattr(self.nc_file.variables['znucl'],module) #---------------------------------------------------- self.goodFile = True #-----------Calculation of some quantities------------# #-----------not available in the output File----------# self.mass_calculation() self.volume_calculation() self.angles_calculation() self.E_kinCalculation() self.temperature_calculation() self.pressure_calculation() #-----------------------------------------------------# except: self.goodFile = False print "Error,unable to read this netcdf file " def E_kinCalculation(self): natom = self.getNatom() mass = self.getMass() vel = self.getVel() self.E_kin= np.zeros(self.Nbtime-1) self.E_kin[:] = np.sum((vel[:,:,0]**2+vel[:,:,1]**2+vel[:,:,2]**2)[:,:] * mass[:] * 0.5,axis=1) def mass_calculation(self): self.mass=[] typat = self.getTypat() amu = self.getAmu() for typ in typat: self.mass.extend([self.amu_emass*amu[int(typ-1)]]) def temperature_calculation(self): Nbstep = self.getNbStep()-1 E_kin = self.getE_kin() self.temperature = np.zeros(Nbstep) for itim in range(Nbstep): self.temperature[itim] = 2.*E_kin[itim]*self.Ha_eV/(3.*self.kb_eVK*self.getNatom()) def pressure_calculation(self): Nbstep = self.getNbStep()-1 TEMPER = self.getTemp() Vol = self.getVol() Stress = self.getStress() natom = self.getNatom() self.pressure= np.zeros(Nbstep) for itim in range(Nbstep): Pionique = (natom/Vol[itim])*self.kb_eVK*TEMPER[itim] * self.fact /self.Ha_eV PP=-(Stress[itim][0]+Stress[itim][1]+Stress[itim][2])/3. + Pionique self.pressure[itim] = PP def volume_calculation(self): Nbstep = self.getNbStep() rprim = self.getRPrim() self.volume = np.zeros(Nbstep) self.volume[:] = rprim[:,0,0]*(rprim[:,1,1]*rprim[:,2,2]-rprim[:,1,2]*rprim[:,2,1]) \ + rprim[:,0,1]*(rprim[:,1,2]*rprim[:,2,0]-rprim[:,1,0]*rprim[:,2,2]) \ + rprim[:,0,2]*(rprim[:,1,0]*rprim[:,2,1]-rprim[:,1,1]*rprim[:,2,0]) def angles_calculation(self): Nbstep = self.getNbStep() rprim = self.getRPrim() self.angles = np.zeros((Nbstep,3)) self.angles[:,0] = (np.arccos((rprim[:,0,0]*rprim[:,1,0]+rprim[:,0,1]*rprim[:,1,1]+rprim[:,0,2]*rprim[:,1,2])\ /np.sqrt((rprim[:,0,0]**2+rprim[:,0,1]**2+rprim[:,0,2]**2)*(rprim[:,1,0]**2+rprim[:,1,1]**2+rprim[:,1,2]**2)))) self.angles[:,1] = (np.arccos((rprim[:,0,0]*rprim[:,2,0]+rprim[:,0,1]*rprim[:,2,1]+rprim[:,0,2]*rprim[:,2,2])\ /np.sqrt((rprim[:,0,0]**2+rprim[:,0,1]**2+rprim[:,0,2]**2)*(rprim[:,2,0]**2+rprim[:,2,1]**2+rprim[:,2,2]**2)))) self.angles[:,2] = (np.arccos((rprim[:,2,0]*rprim[:,1,0]+rprim[:,2,1]*rprim[:,1,1]+rprim[:,2,2]*rprim[:,1,2])\ /np.sqrt((rprim[:,2,0]**2+rprim[:,2,1]**2+rprim[:,2,2]**2)*(rprim[:,1,0]**2+rprim[:,1,1]**2+rprim[:,1,2]**2)))) #-----------------------------# #-------ACCESSOR--------------# #-----------------------------# def isGoodFile(self): # return Boolean True/False according to the HIST and NC file: return self.goodFile def getNameFile(self,fullpath=False): # return string with the name of the file if self.goodFile: if fullpath : return self.namefile2 else: name = self.namefile2.split('/') return name[len(name)-1] else: return 0 def getNatom(self): # return integer with the number of atom if self.goodFile: return self.natom else: return 0 def getVol(self): # return double with the value of the volume (in bohr^3) if self.goodFile: return self.volume[self.ni:self.nf] # Temporarily. (ni-1) for start slicing at 0 else: return 0 def getTypat(self): # return 1D array with the type of particules if self.goodFile: return self.typat else: return 0 def getNTypat(self): # return the number of type of particules if self.goodFile: return np.max(self.typat) else: return 0 def getAmu(self): # return 1D array with the atomic mass if self.goodFile: return self.amu else: return 0 def getZnucl(self): # return 1D array with the number atomic of all the atoms if self.goodFile: return self.znucl else: return 0 def getMass(self): # return 1D array with the mass of all particules if self.goodFile: return self.mass else: return 0 def getNbTime(self): # return integer with the number of step in the simulation if self.goodFile : return self.Nbtime else: return 0 def getNImage(self): # return integer with the number of image in the simulation if self.goodFile : return self.n_image else: return 0 def getRPrim(self,image = 1): # return the primitive vector (Borh) if self.goodFile : return self.rprimd[self.ni-1:self.nf,image-1] else: return 0 def getNbStep(self): # return the number of step between ni and nf if self.goodFile: return self.nf-self.ni+1 else: return 0 def getDtion(self): # return ion time steps in atomic units of time if self.goodFile: return self.dtion else: return 0 def getVel(self,image = 1): # return 3D array with the velocity (vx,vy,z) for all particule at each time ( V[t][[at][vx,vy,vz]) for one image # between ni and nf if self.goodFile: if self.n_image == 1: return self.vel[self.ni:self.nf,image-1] # Temporarily because vel[t=0] = 0. put (ni-1) for start slicing at 0 elif self.n_image != 1 : return self.getVelCentroid() else: return 0 def getXCart(self,image = 1): # return 3D array with the position (x, y, z) for all particule at each time ( pos[t][at][x,y,z]) for one image in cartesian coordiantes # between ni and nf if self.goodFile: if self.n_image == 1: return self.xcart[self.ni:self.nf,image-1] # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 elif self.n_image != 1 : return self.getXCartCentroid() else: return 0 def getXRed(self,image = 1): # return 3D array with the position (x, y, z) for all particule at each time ( pos[t][at][x,y,z]) for one image in reduced coordinates # between ni and nf if self.goodFile: if self.n_image == 1: return self.xred[self.ni:self.nf,image-1] # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 elif self.n_image != 1 : return self.getXCartCentroid(xred=True) else: return 0 def getFCart(self,image = 1): # return 3D array with the forces (x, y, z) for all particule at each time ( pos[t][at][x,y,z]) for one image in cartesian coordinates # between ni and nf if self.goodFile: if self.n_image == 1: return self.fcart[self.ni:self.nf,image-1] # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 elif self.n_image != 1 : return self.getXCartCentroid(xred=True) else: return 0 def getFRed(self,image = 1): # return 3D array with the forces (x, y, z) for all particule at each time ( f[t][at][x,y,z]) for one image in reduced coordinates # between ni and nf if self.goodFile: if self.n_image == 1: return self.fred[self.ni:self.nf,image-1] # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 elif self.n_image != 1 : return self.getXCartCentroid(xred=True) else: return 0 def getVelCentroid(self): # return 3D array with the velocity (vx,vy,z) for all particule at each time ( V[t][[at][vx,vy,vz]) for the centroid # between ni and nf if self.goodFile: return np.sqrt(np.mean(self.vel[self.ni:self.nf,:]**2,axis=1)) # Temporarily because vel[t=0] = 0. put (ni-1) for start slicing at 0 else: return 0 def getXCartCentroid(self,xred=False): # return 3D array with the position (x, y, z) for all particule at each time ( pos[t][at][x,y,z]) for the centroid # between ni and nf if self.goodFile: if(xred): return np.mean(self.xred[self.ni:self.nf,:],axis=1) # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 else: return np.mean(self.xcart[self.ni:self.nf,:],axis=1) # Temporarily because pos[t=0] = 0. put (ni-1) for start slicing at 0 else: return 0 def getE_pot(self): # return 1D array with the potential energy (Ha) if self.goodFile : return self.E_pot[self.ni-1:self.nf-1] #Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return 0 def getMdtime(self): if self.goodFile : return self.mdtime[self.ni-1:self.nf-1] #Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return 0 def getE_kin_ion(self): # return 1D array with the Energy KINetic ionic (Ha) if self.goodFile : return self.E_kin_ion[self.ni-1:self.nf-1] #Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return 0 def getE_kin(self): # return 1D array with the kinetic energy (Ha) if self.goodFile : return self.E_kin[self.ni-1:self.nf-1]#Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return 0 def getE_Tot(self): # return 1D array with the total energy (Ha) between ni and nf if self.goodFile : E_kin = self.getE_kin() E_pot = self.getE_pot() Nbstep = self.getNbStep() ETOT = E_kin[0:Nbstep] + E_pot[0:Nbstep] return ETOT else: return 0 def getTemp(self): # return 1D array with the Temperature (K) if self.goodFile : return self.temperature[self.ni-1:self.nf-1]#Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return 0 def getIonMove(self): if self.goodFile : self.ionMove = [12] return self.ionMove else: return 0 def getStress(self,image = 1): # return 2D array with the Stess (GPa) if self.goodFile: if self.n_image==1: return self.stress[self.ni-1:self.nf-1,image-1] * 29421.033e0#Temporarily. put (nf) for finish slicing at 299, else: #Skip the last step to get the same array size (see vel and pos). return self.stress[self.ni-1:self.nf-1,image-1] #In PIMD, stress is already give in GPA# else: return 0 def getAcell(self,image = 1): if self.goodFile: acell = np.zeros(((self.nf-self.ni),image,3)) for i in range(3): acell[:,image-1,i]=np.sqrt(self.rprimd[self.ni:self.nf,image-1,i,0]**2+\ self.rprimd[self.ni:self.nf,image-1,i,1]**2+\ self.rprimd[self.ni:self.nf,image-1,i,2]**2) return acell[:,image-1,:] # Temporarily. (ni-1) for start slicing at 0 else: return 0 def getAngles(self): if self.goodFile: return self.angles[self.ni:self.nf] * 180 / np.pi # Temporarily. (ni-1) for start slicing at 0 else: return 0 def getPress(self): # return 1D array with the Pressure between ni and nf (GPa) if self.goodFile : return self.pressure[self.ni-1:self.nf-1] else: return 0 def getMoy(self,data): #return average of data if self.goodFile : nb = len(data) M=0.0 for i in range(nb): M = M + data[i] M = M / nb return M else: return 0 def getStandardDeviation(self,data ,averageData = 0): res = 0 if averageData == 0: averageData = self.getMoy(data) for i in range(len(data)): res += (data[i] - averageData)**2 res /= len(data) res = res**0.5 return res def getAtomName(self): #return array with the name of the atom PTOE = Analysis.PeriodicTableElement() znucl = self.getZnucl() name = [] for i in range(len(znucl)): name.append(PTOE.getName(znucl[i])) return name def getNi(self): #return ni (integer) return self.ni def getNf(self): #return nf (integer) return self.nf def setNi(self,pNi): #to set ni (integer) self.ni = pNi def setNf(self,pNf): #to set nf (integer) self.nf = pNf #-----------------------------# #-----------------------------# #-----------------------------# #######GROUND STATES OUTPUT (BETA)############# class outputFile: #-----------------------------# #-------CONSTRUCTOR-----------# #-----------------------------# Ha_eV = 27.2113834e0 # Conversion factor 1 Hartree = 27,2113834 eV def __init__(self, pnamefile): self.V = [] self.P = [] self.E = [] self.ecut = [] self.Nbk = [] self.ionMove = [] self.a = [] self.stress = {} self.begin = True self.date = "" self.hour = "" self.name = "" self.goodFile = False self.namefile ="" #-------READ THE LAST OUTPUT FILE OR THE CHOOSEN FILE--------# if pnamefile == "": OUT_list = [] filesdir = os.listdir(os.getcwd()) for files in filesdir : if os.path.isfile(files): stats = os.stat(files) fic_tuple = time.localtime(stats[8]), files lili = string.split(files, '.') if files.find('.out') == len(files)-4 and len(files)>4 : OUT_list.append(fic_tuple) if len(OUT_list) > 0 : OUT_list.sort() ; OUT_list.reverse() fic=OUT_list[0][1] if os.path.exists(fic): self.namefile = fic else: self.namefile = str(pnamefile) #-------------------------------------------------------------# if self.namefile !="": try : #if 1==1: self.ionMove = np.array(commands.getoutput('grep \"ionmov =\" '+ self.namefile + ' | awk \'{print $6}\'').split(), dtype=float) for i in range(len(self.ionMove)): if self.ionMove[i] == 12 : return self.date = commands.getoutput('grep \"Starting date\" ' + self.namefile + ' | awk \'{print $4,$5,$6,$7}\'') self.hour = commands.getoutput('grep \"( at\" '+ self.namefile + ' | awk \'{print $4}\'') self.Nbk = np.array(commands.getoutput('grep \"nkpt =\" '+ self.namefile + ' | awk \'{print $12}\'').split(), dtype=float) self.V = np.array(commands.getoutput('grep \"Unit cell volume ucvol\" '+ self.namefile + ' | awk \'{print $5}\'').split(), dtype=float) self.P = np.array(commands.getoutput('grep \"Pressure\" '+ self.namefile + ' | awk \'{print $8}\'').split(), dtype=float) self.E = np.array(commands.getoutput('grep \"etotal\" '+ self.namefile + '|sed /\"Total energy\"/d | awk \'{print $2}\'').split(), dtype=float) self.ecut = np.array(commands.getoutput('grep \"ecut(hartree)\" '+ self.namefile + ' | awk \'{print $2}\'').split(), dtype=float) self.a = np.array(commands.getoutput('grep -A '+ str(len(self.E)+1)+ ' \"after computation\" '+ self.namefile+ '| sed \'/acell/!d\' | awk \'{print $2}\'').split(), dtype=float) self.b = np.array(commands.getoutput('grep -A '+ str(len(self.E)+1)+ ' \"after computation\" '+ self.namefile+ '| sed \'/acell/!d\' | awk \'{print $3}\'').split(), dtype=float) self.c = np.array(commands.getoutput('grep -A '+ str(len(self.E)+1)+ ' \"after computation\" '+ self.namefile+ '| sed \'/acell/!d\' | awk \'{print $4}\'').split(), dtype=float) #Reading stress of particules: temp = commands.getoutput(' grep -A 3 \"stress tensor (GPa)\" '+ self.namefile + '|sed \'/Car/d\' |sed \'/--/d\' | awk \'{print $4,$7}\' ') temp = temp.split() s = (len(self.E),6) temp2 = np.zeros(s) for i in range(len(self.E)): temp2[i] = np.array( [ temp[(6*i)+0], temp[(6*i)+2], temp[(6*i)+4], temp[(6*i)+1], temp[(6*i)+3], temp[(6*i)+5] ], dtype=float) self.stress = temp2 #Reading the primitive vectors: temp = commands.getoutput(' grep -A 3 \"Real(R)+Recip(G)\" '+ self.namefile + ' | sed \'/space primitive vectors,/d\'|sed \'/--/d\' | awk \'{print $2,$3,$4 }\' ') temp = temp.split() s = (len(self.E),3,3) temp2 = np.zeros(s) for i in range(len(self.E)): a = np.array( [temp[(9*i)+0], temp[(9*i)+1], temp[(9*i)+2]], dtype=float) b = np.array( [temp[(9*i)+3], temp[(9*i)+4], temp[(9*i)+5]], dtype=float) c = np.array( [temp[(9*i)+6], temp[(9*i)+7], temp[(9*i)+8]], dtype=float) temp2[i] = np.array( [a,b,c], dtype=float) self.rprimd =temp2 if len(self.E) > 0: self.goodFile = True else : try : QtGui.QMessageBox.critical(None,"Warning","Your output file is not correct") except: print "Your output file is not correct" if len(self.a) == 1: self.a = np.zeros(len(self.E))+self.a[0] self.b = np.zeros(len(self.E))+self.b[0] self.c = np.zeros(len(self.E))+self.c[0] if len(self.ecut) == 1: self.ecut = np.zeros(len(self.E))+self.ecut[0] except: print "Unable to read file, please contact abinit Group" QtGui.QMessageBox.critical(None,"Warning","Unable to read file, please contact abinit Group") #-----------------------------# #-------ACCESSOR--------------# #-----------------------------# def getNameFile(self): # return string with the name of the file if self.goodFile: name = self.namefile.split('/') return name[len(name)-1] else: return 0 def getNbDataset(self): if self.goodFile: return len(self.E) else: return 0 def getE_Tot(self): if self.goodFile: return np.array(self.E) else: return 0 def getVol(self): if self.goodFile: return np.array(self.V) else: return 0 def getSigma(self): if self.goodFile: return self.sigma else: return 0 def getAcell(self): if self.goodFile: return self.acell else: return 0 def getPress(self): if self.goodFile: return np.array(self.P) else: return 0 def getNbK(self): if self.goodFile: return np.array(self.Nbk) else: return 0 def getIonMove(self): if self.goodFile: return np.array(self.ionMove) else: return 0 def getRPrim(self): if self.goodFile: return self.rprimd else: return 0 def getEcut(self): if self.goodFile: return np.array(self.ecut) else: return 0 def getA(self): if self.goodFile: return self.a else: return 0 def getB(self): if self.goodFile: return self.b else: return 0 def getC(self): if self.goodFile: return self.c else: return 0 def getDate(self): if self.goodFile: return self.date else: return 0 def getStress(self): if self.goodFile: return self.stress else: return 0 def getHour(self): if self.goodFile: return self.hour else: return 0 def isGoodFile(self): return self.goodFile