Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #!/software/languages/python/2.5.2/gnu-4.2.4/bin/python
- import sys,re,numpy
- filename = sys.argv[1] #expects 1 argument, realative path to input file
- energy = [] #List of extracted energies
- rxcoord = [] #List of reaction coordinates
- cols = {} #Internal varaiable to track column names
- cols[1] = "ENERGY" #expect that energy and rxcoord are always first two colums (i think this is correct in all cases!)
- cols[2] = "RX.COORD"
- geoms = {} #Dictionary containing a list of all coordinates unsorted geoms{(step, X|Y|Z, atom#)} = coordinate
- noSteps = 0 # counters used to create numpy array at end
- noAtoms = 0
- try:
- inputfile = open(filename, "r")
- except:
- print ("file cannot be opened")
- for line in inputfile:
- if "Summary of reaction path following:" in line: #look for header line in input file
- waste = inputfile.next() # waste some lines
- hyphens = inputfile.next()
- line = inputfile.next()
- while re.search("(ENERGY\s+|RX\.COORD\s+|[XYZ]\d+\s+)+",line): #look for rows that contain column names
- numCol = 0
- match = re.findall("ENERGY\s+|RX\.COORD\s+|[XYZ]\d+\s+",line) #extract column names
- for col in match: #loop over column names
- numCol +=1
- col = col.strip(' \t\n\r')
- res = re.findall("([XYZ])(\d+)", col)
- for items in res:
- cols[numCol] = (items[0], items[1]) # store column name as tuple (X|Y|Z, atom#)
- if int(items[1]) > noAtoms: noAtoms = int(items[1]) #update the noAtoms counter if required
- line = inputfile.next()
- while re.search("\d+\s+([-|+]?\d+\.\d+\s+)+",line): #loop over the coordinate lines
- broken = line.split()
- if int(broken[0]) > noSteps: noSteps = int(broken[0]) #update noSteps counter if required
- for i in range(1,numCol+1,1):
- if i == 1 and cols[i] == "ENERGY": energy.append(broken[i]) #extract energies into list
- if i == 2 and cols[i] == "RX.COORD": rxcoord.append(broken[i]) #extract reaction coordinates into list
- if cols[i] != "ENERGY" and cols[i] != "RX.COORD":
- geoms[(broken[0],cols[i][0], cols[i][1])] = broken[i] #extract coordinates into dictionary
- try: #if this fails we have reached the end of the output
- line = inputfile.next()
- except:
- break # so we break out of this loop
- cols = {} #reset the colum names
- output=numpy.empty((noSteps,noAtoms,3)) #create empty array of rank 3. output[Step Number[Atom Number[X Y Z]]]
- for each in geoms: #loop over lines in geometry list
- if each[1] == "X": axes = 0 #number the axes for the array
- if each[1] == "Y": axes = 1
- if each[1] == "Z": axes = 2
- output[int(each[0])-1,int(each[2])-1,axes] = geoms[each] #fill array position with coordinate value
- for i in range(noSteps): #print each step
- print output[i,:,:]
Add Comment
Please, Sign In to add comment