Jump to content

Talk:Mod:Hunt Research Group/build freq

From ChemWiki

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()
#