#!/usr/bin/env python # -*- coding: utf-8 -*- """ Created on Wed Dec 9 18:21:43 2015 @author: tam10 """ import sys def irc2xyz(input_file,print_energies=False,csv_name="",kj=False,reverse_order=False, write_reaction_coordinate=False): """ Converts IRC log file to an xyz collection input_file: (Required) Path to file to convert print_energies: (Optional) Print energies to stdout csv_name: (Optional) Path for csv output (Energies and dipole moments) kj: (Optional) False: (Default) csv energies in kcal/mol; True: csv energies in kj/mol reverse_order: (Optional) Reverse order of xyz frames """ coordinates=[] coordinates_list=[] forces=[] forces_list=[] energies=[] energies_list=[] charges_symbols=[] charges_symbols_list=[] dipoles=[] dipoles_list=[] reaction_coordinates=[] reaction_coordinates_list=[] expect_coordinates=False coordinates_mode=False expect_forces=False forces_mode=False expect_mulliken=False mulliken_mode=False dipole_mode=False point_number="0" path_number="0" au_kcal=627.50960803059 au_kj=2625.5002 print_energies=print_energies in (True,"True") kj=kj in (True,"True") reverse_order=reverse_order in (True,"True") with open(input_file) as irc_file: irc_lines = irc_file.readlines() for line_number,irc_line in enumerate(irc_lines): split_line=irc_line.split() if coordinates_mode: #Using booleans to determine where the parser is in the file e.g. is the next section a list of coordinates? if split_line[0].isdigit(): coordinates.append([split_line[0]]+split_line[-3:]) else: coordinates_mode=False if forces_mode: if split_line[0].isdigit(): forces.append([split_line[0]]+split_line[-3:]) else: forces_mode=False #This line might fail if Mulliken charges aren't present. Would need to find another way to get symbol mappings if mulliken_mode: #Using Mulliken charges section to get symbol mapping for xyz file (Gaussian deals with atom number) if split_line[0].isdigit(): charges_symbols.append(split_line) else: mulliken_mode=False if dipole_mode: try: dipoles.append([str(float(split_line[-1]))]) except: dipole_mode=False if "SCF Done:" in irc_line: #Energy of each step is taken from here energies.append([str(float(split_line[4]))]) if expect_coordinates: if "----" in irc_line: #This should always be true if expect_coordinates is true, but there's a chance this word might be used in titlecard coordinates_mode=True expect_coordinates=False if expect_forces: if "----" in irc_line: forces_mode=True expect_forces=False if expect_mulliken: expect_mulliken=False mulliken_mode=True if "Point Number:" in irc_line: #Here we get the path and point number for each geometry point_number=split_line[2] path_number=split_line[5] try: reaction_coordinates_list.append([path_number,point_number,irc_lines[line_number+2].split()[8]]) except IndexError: reaction_coordinates_list.append([path_number,point_number,'0.00000']) [coordinates_list.append([path_number, point_number]+c) for c in coordinates] [forces_list.append([path_number, point_number]+f) for f in forces] if len(charges_symbols)==max([int(e[2]) for e in coordinates_list]): charges_symbols_list=charges_symbols [energies_list.append([path_number, point_number]+e) for e in energies] [dipoles_list.append([path_number, point_number]+d) for d in dipoles] coordinates=[] forces=[] charges_symbols=[] energies=[] dipoles=[] if len(split_line)>2 and "Coordinates" in split_line: expect_coordinates=True if len(split_line)>2 and "Forces" in split_line and not "Forces:" in split_line: expect_forces=True if "Mulliken charges:" in irc_line: expect_mulliken=True if "Dipole moment" in irc_line: dipole_mode=True number_of_atoms=max([int(e[2]) for e in coordinates_list]) paths=list(set([int(e[0]) for e in energies_list])) path_point_list=[] for path in paths: #Create a list of points for each path. If using both directions, the other direction is reversed. if reverse_order: if path==1: points=list(set([int(e[1]) for e in energies_list if int(e[0])==path])) points.reverse() else: points=list(set([int(e[1]) for e in energies_list if int(e[0])==path])) path_point_list.append(points) else: if path==2: points=list(set([int(e[1]) for e in energies_list if int(e[0])==path])) points.reverse() else: points=list(set([int(e[1]) for e in energies_list if int(e[0])==path])) path_point_list.append(points) ordered_paths=[2,1][-len(path_point_list):] #odd syntax so path_point_list[2] isn't used if the irc is forward only if reverse_order: ordered_paths.reverse() with open(input_file.split(".")[0]+".xyz","w") as xyz_output: for path in ordered_paths: for point in path_point_list[path-1]: energy=[e[2] for e in energies_list if e[0]==str(path) and e[1]==str(point)] reaction_coordinate=[r[2] for r in reaction_coordinates_list if r[0]==str(path) and r[1]==str(point)] if path==2: sign_str='-' else: sign_str='' xyz_output.write(str(number_of_atoms)+'\n') if write_reaction_coordinate: xyz_output.write(" Energy = "+"{0: >20}".format(energy[0])+" Reaction Coordinate = "+"{0: >20}".format(sign_str+reaction_coordinate[0])+'\n') else: xyz_output.write(" Energy = "+"{0: >20}".format(energy[0])+'\n') if print_energies: print(" Energy = "+"{0: >20}".format(energy[0])+" Path = "+"{0: <2}".format(str(path))+" Point = "+str(point)) for atom in range(number_of_atoms): symbol=[cs[1] for cs in charges_symbols_list if cs[0]==str(atom+1)][0] pos_list=[c[-3:] for c in coordinates_list if c[0]==str(path) and c[1]==str(point) and c[2]==str(atom+1)][0] force_list=[f[-3:] for f in forces_list if f[0]==str(path) and f[1]==str(point) and f[2]==str(atom+1)][0] xyz_output.write("{0: >2}".format(symbol)+"".join(["{0: >14}".format(p) for p in pos_list])+" 0.00"+"".join(["{0: >14}".format(f) for f in force_list])+'\n') if len(csv_name)>0: normalise_to=True normalise_au=float(0) with open(csv_name,"w") as csv_output: if kj: csv_output.write("Reaction Coordinate,Energy (kj/mol),Dipole Moment\n") else: csv_output.write("Reaction Coordinate,Energy (kcal/mol),Dipole Moment\n") for path in ordered_paths: for point in path_point_list[path-1]: au=float([e[2] for e in energies_list if e[0]==str(path) and e[1]==str(point)][0])-normalise_au if normalise_to: normalise_au=au au=float(0) normalise_to=False try: dipole=[d[2] for d in dipoles_list if d[0]==str(path) and d[1]==str(point)][0] except: dipole="Unkown" if kj: energy=au*au_kj else: energy=au*au_kcal if path==2: csv_output.write(",".join(["-"+str(point),str(energy),dipole,'\n'])) else: csv_output.write(",".join([str(point),str(energy),dipole,'\n'])) if __name__ == '__main__': args = sys.argv[1:] if len(args) == 0: print(" \ Converts IRC log file to an xyz collection\n \ Arguments:\n \ 0 (str) : Path to file to convert\n \ 1 (bool): (Optional) Print energies to stdout\n \ 2 (str) : (Optional) Path for csv output (Energies and dipole moments)\n \ 3 (bool): (Optional) False: (Default) csv energies in kcal/mol; True: csv energies in kj/mol\n \ 4 (bool): (Optional) Reverse frame order\n \ 5 (bool): (Optional) Print reaction coordinate alongside energy") if len(args) == 1: irc2xyz(args[0]) elif len(args) == 2: irc2xyz(args[0],args[1]) elif len(args) == 3: irc2xyz(args[0],args[1],args[2]) elif len(args) == 4: irc2xyz(args[0],args[1],args[2],args[3]) elif len(args) == 5: irc2xyz(args[0],args[1],args[2],args[3],args[4]) elif len(args) == 6: irc2xyz(args[0],args[1],args[2],args[3],args[4],args[5])