python script to build frequency file from optimisation output
- copy the code into a script called "energy.py"
- to execute the script type "python build_freq.py file_name step
- if you have your execute directory in your path, you can alter the first line to your python executable and then you can just type "build_freq.py filename step"
- if step is not given the last structure from the optimisation file will be used
- if step is given the geometry from that step will be used
- NOTE: you will need to alter the string which prints the job key words into the gaussian input file
#!/opt/local/bin/python3
#
# build_freq.py
#
# python script to extract a given step structure, or the last structure from a gaussian log file and build a frequency input file
# to run the script type:
#
# python build_freq.py input_file_name step_no
#
# if no step_no is given the last geometry will be used
#
import sys
import os
dir=os.getcwd()
#==set up FILES ==
# if there are 2 arguments to use the last step (of the optimisation)
# if there are 3 arguments then use a specific step_no from the optimisation
if len(sys.argv) == 2:
infile=str(sys.argv[1])
step_no=int(0)
elif len(sys.argv) == 3:
infile=str(sys.argv[1])
step_no=int(sys.argv[2])
#since counting starts from zero
step_no=step_no-1
#endif
# if there is less or more arguments something is wrong
else:
print('to run the script please type: python build_freq.py input_file_name step_no')
sys.exit()
#endif
# check to see if the user has used the .log extenion or not
# find returns -1 if the sub-string is not found, or the index of where the sub-string starts
# we cut out the .log part if .log is found
i_infile=infile.find('.log')
if (i_infile == -1):
base=infile
else:
base=infile[:i_infile]
#endif
# now set up the base file names
log_file=base+'.log'
geom_file=base+'.xyz'
# if opt is not in the filename add freq to the filename
# if opt is already in the file change this to freq
i_opt=base.find('opt')
i_optf=base.find('optf')
if 'opt' in base:
print ('yes')
if 'optf' in base:
base2=base.replace('optf','freq')
else:
base2=base.replace('opt','freq')
# endif
else:
base2=base + '_freq'
#endif
com_file=base2+'.com'
# print out some reference information
s='directory is '
print('{0:}{1:}'.format(s,dir))
s='input file is '
print('{0:}{1:}'.format(s,infile))
s='output file is '
print('{0:}{1:}'.format(s,com_file))
#sys.exit()
#
#==set some DEFAULTS ==
col_num=[]
col_atom=[]
col_type=[]
col_x=[]
col_y=[]
col_z=[]
rep=0
n=0
go=False
atom_total=0
geom_count=0
#
#==start READING THE FILE ==
# open file and read lines one by one
f=open(log_file,"r")
# line is a variable which contains the line
line =f.readline()
while line :
line =f.readline()
# print('{0:}'.format(line))
#==extract GEOMETRY ==
# check for standard orientation text set 'go' if the text appears
# geom_count is the number of the geometry stored
# atom_total will be the number of atoms
# n will be stepping through the lines in std orientation section so set to zero for each time through
# if this is the first geometry print that a geometry is found
if 'Standard orientation:' in line:
go=True
s='Structure in standard orientation found '
n=1
geom_count=geom_count+1
if geom_count == 1 : print('{0:}'.format(s))
s='geometry count= '
# print('{0:}{1:}'.format(s,geom_count))
# endif
#
# if go is on then we need to extract and store the structure
# skip the next 4 lines, on the 5th start collecting structure
# m is our switch to stop, we stop once we get to the Rotational constants line
# m is also our switch to not print the lines of dashes
# we need to know how many atoms there are in total as well
# we simply append each data list into the arrays col which has components num,atom,type,x,y,z
while go:
n = n+1
m = 0
line =f.readline()
if 'Rotational constants (GHZ):' in line: m=-1
if '---------' in line: m=1
if n > 4 and m == 0 :
# print('line: {0:}'.format(line))
if geom_count == 1 : atom_total=atom_total+1
a1,b1,c1,d1,e1,f1=line.rstrip().split()
col_num.append(int(a1))
col_atom.append(int(b1))
col_type.append(int(c1))
col_x.append(float(d1))
col_y.append(float(e1))
col_z.append(float(f1))
if m == -1 : go=False
# close while
#close while
#== build COM FILE ==
# first give all the comp details
c=open(com_file,"w")
s="%nprocshared=48\n%mem=122000MB\n"
c.write('{0:}'.format(s))
s='%chk='+base2+'.chk\n'
c.write('{0:}'.format(s))
s="# freq m06/6-311G(d,p) geom=cartesian int=ultrafine scf=conver=10 \n\n"
c.write('{0:}'.format(s))
s="Title Card Required\n\n0 1\n"
c.write('{0:}'.format(s))
#other options are
#s="# freq b3lyp/gen pseudo=cards empiricaldispersion=gd3bj geom=cartesian \n "
#s="gfinput int=ultrafine scf=conver=10 test\n\n"
# now put in the coordinates
# if step_no is zero we only use the last geometry
# if step_no is given, we use that geometry
# geom_count is the number of geometries
# atom_total is the numer of atoms
# if taking the last geometry count backwards
# n is the internal counter
# n_start defines where to start printing
# n_stop defines wher to stop printing
#
n=0
n_start=0
n_stop=0
#print ('{0:}{1:}'.format('step_no=',step_no))
#print ('{0:}{1:}'.format('atom_total=',atom_total))
#print ('{0:}{1:}'.format('geom_count=',geom_count))
if step_no == 0 :
n_start=len(col_num)-atom_total
n_stop=len(col_num)
s=' structures found only the last is used'
print('{0:}{1:}'.format(geom_count,s))
else:
n_start=atom_total*step_no
n_stop=n_start+atom_total
print('{0:}{1:}{2:}{3:}'.format('using structure ',step_no+1,' of ',geom_count))
# endif
#print ('{0:}{1:}'.format('n_start=',n_start))
#print ('{0:}{1:}'.format('n_stop=',n_stop))
n=n_start
while n >= n_start and n < n_stop:
c.write('{0:>4d} {1: 09.6f} {2: 09.6f} {3: 09.6f} \n'.format(col_atom[n],col_x[n],col_y[n],col_z[n]))
# print('{0:>4d} {1:>4d} {2: 09.6f} {3: 09.6f} {4: 09.6f}'.format(col_num[n],col_atom[n],col_x[n],col_y[n],col_z[n]))
n=n + 1
#close while
# now put in the gen information
#s="\nCl 0\n6-311+G(d,p)\n****\n"
#t="In 0\nLANL2DZ\n****\n\n"
#u="In 0\nLANL2DZ\n\n\n"
#c.write('{0:}{1:}{2:}'.format(s,t,u))
# and add the last blank lines
s=" \n\n"
c.write('{0:}'.format(s))
# close all open files
f.close()
c.close()
#
