#!/usr/bin/env python # -*- coding: utf-8 -*- """ Created on Wed Feb 3 12:24:35 2016 @author: tam10 """ import numpy as np import csv import argparse import os import sys xyz_file="" def get_measurements(input_file,atoms_to_measure,output_file='output.csv'): xyz_file=get_xyz(input_file) output_file=check_csv(output_file,input_file) with open(xyz_file,'r') as irc: irc_lines=[l.strip('\n').split() for l in irc.readlines()] mol_collection,reaction_coordinates_list,symbol_list,size=get_coordinates(irc_lines) measurement_list=measure(atoms_to_measure,mol_collection,symbol_list,size) if sys.version_info<=(3,0): output=open(output_file,'wb') else: output=open(output_file,'w',newline='') csv_writer=csv.writer(output) for index,line in enumerate(measurement_list[0]): if index==0: csv_writer.writerow(['Reaction Coordinate']+[[l for i,l in enumerate(m) if i==index][0] for m in measurement_list]) else: csv_writer.writerow([reaction_coordinates_list[index-1]]+[[l for i,l in enumerate(m) if i==index][0] for m in measurement_list]) output.close() def get_xyz(input_file): if input_file.split('.')[-1].lower()=='log': try: from irc2xyz import irc2xyz irc2xyz(input_file=input_file,write_reaction_coordinate=True) xyz_file='.'.join(input_file.split('.')[:-1]).lower()+'.xyz' except ImportError: raise ImportError("irc2xyz.py not found. Make sure irc2xyz.py is in the same directory as pairdistances.py") elif input_file.split('.')[-1].lower()=='xyz': xyz_file=input_file else: raise IOError('Input must be an xyz file or irc log file: \''+input_file+'\'') return xyz_file def check_csv(output_file,input_file): if not output_file: return ".".join(input_file.split('.')[:-1]) + '.csv' if output_file.split('.')[-1].lower()!='csv': raise IOError('Output must be a csv file: \''+output_file+'\'') if not os.path.isabs(output_file): output_file=os.sep.join(os.path.realpath(input_file).split(os.sep)[:-1]+[output_file]) return(output_file) def get_coordinates(irc_lines): mol_index,mol_collection,reaction_coordinates_list=-1,[],[] for l in irc_lines: if l[0].isdigit(): #First line represents number of atoms in molecule size,atom_index,symbol_list=int(l[0]),0,[] mol_collection.append([]) #This allows appending to a list within mol_collection mol_index+=1 else: if atom_index!=0: mol_collection[mol_index].append([float(c) for c in l[1:4]]) #Append coordinates symbol_list.append([l[0]]) #Append atom symbols assert size>=atom_index,'Error parsing file: \''+xyz_file+'\'' #Check there aren't more atoms than there should be elif atom_index==0: try: reaction_coordinates_list.append(l[6]) except IndexError: reaction_coordinates_list.append(mol_index+1) atom_index+=1 return mol_collection,reaction_coordinates_list,symbol_list,size def measure(atoms_to_measure,mol_collection,symbol_list,size): measurement_list=[] for atoms in atoms_to_measure: for atom_index in atoms: assert size>=atom_index,'Atom index \''+str(atom_index)+'\' greater than size of molecule ('+str(size)+'). Check input' ps=[] for mol in mol_collection: title_string=['-'.join([symbol_list[a-1][0]+str(a) for a in atoms])] #Title string in the format "C1-C2-H3" ps.append([np.array(mol[i-1]) for i in atoms]) if len(atoms)==2: distances=[round(np.linalg.norm(p[0]-p[1]),3) for p in ps] #Using numpy to subtract lists and calculate distances measurement_list.append(title_string+distances) elif len(atoms)==3: angles=[round(np.degrees(np.arccos(np.dot(p[0]-p[1], p[2]-p[1]) / (np.linalg.norm(p[0]-p[1]) * np.linalg.norm(p[2]-p[1])))),2) for p in ps] measurement_list.append(title_string+angles) elif len(atoms)==4: dihedrals=[round(np.degrees(np.arccos(np.dot(p[0]-p[1], p[3]-p[2]) / (np.linalg.norm(p[0]-p[1]) * np.linalg.norm(p[3]-p[2])))),2) for p in ps] measurement_list.append(title_string+dihedrals) return measurement_list if __name__=='__main__': parser=argparse.ArgumentParser(description="Generate a csv file (-o ) containing measurements between atoms (-l ) from an input xyz collection (-i )", epilog="Example:\n./pairdistances.py -i TS_IRC.xyz -l 2,4 1,12 2,4,5 1,2,4,5 -o TS_IRC.csv\n") parser.add_argument('-i', '--input', type=str, required=True, help="Input file to take geometries from. Must be an xyz collection or irc log file (irc2xyz.py required)", dest='input_name_arg') parser.add_argument('-l', '--list', nargs='+', required=True, help="List of atom measurements. Atoms within measurements must be comma separated, measurements must be separated with a space", dest='pair_list_arg') parser.add_argument('-o', '--output', type=str, help="Output .csv file name", dest='output_name_arg', default=None) args=parser.parse_args() try: input_name=args.input_name_arg except TypeError: parser.print_help() exit() try: pair_list=args.pair_list_arg pair_list=[[int(p) for p in pair.split(',')] for pair in pair_list] except TypeError: parser.print_help() exit() try: output_name=args.output_name_arg except TypeError: get_measurements(input_name,pair_list) exit() get_measurements(input_name,pair_list, output_name)