Talk:Mod:Hunt Research Group/Dbuild freq file

getdata script

• script to extract extract the last structure from a specific group of gaussian log files and build the corresponding frequency input files
• simply copy and paste this script into a file Dbuild_freq.py
• then run with "python3 Dbuild_freq.py no_files start_no input_file_name"

#!/usr/bin/env python3
#
# Dbuild_freq.py 
#
# python script to extract the last structure from a specific group of gaussian log files and 
# build the corresponding frequency input files
#
# to run the script type python3 Dbuild_freq.py no_files start_no input_file_name (without extension)
# no_files is the number of files to examine
# start_no is the starting file number
# use xx to replace the number of the file
#
import sys
import os
dir=os.getcwd()
 
if len(sys.argv) < 4:
  print('to run the script please type: python bfreq.py no_files start_no input_file_name \(without extension\)')
  sys.exit()
else:
  n_jobs=int(sys.argv[1])
  n_start=int(sys.argv[2])
  bjob=str(sys.argv[3])
  s='directory is '
  print('{0:}{1:}'.format(s,dir))
  ss='starting at '
  s=' cycle over '
  print('{0:}{1:}{2:}{3:}'.format(ss,n_start,s,n_jobs))
  s='base job expression: '
  print('{0:}{1:}'.format(s,bjob))
#close if
#
# set some defaults
nstore=[]
k=0
kk=0

# cycle over log files
while k < n_jobs :
  kk=n_start+k
  job=bjob.replace('xx',str(kk))
  joblog=dir+"/"+job+".log"
  s='searching file '
  print('{0:}{1:}'.format(s,joblog))
# set some defaults
  col_num=[]
  col_atom=[]
  col_type=[]
  col_x=[]
  col_y=[]
  col_z=[]
  go=False
  rep=0
# now set up the base file names
# 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=job.find('opt') 
  i_optf=job.find('optf') 
  if 'opt' in job:
    if 'optf' in job:
       job2=job.replace('optf','freq')
    else:
       job2=job.replace('opt','freq')
#   endif
  else:
    job2=job + '_freq'
#  endif
  out=job2+'.com'
  s='output file is '
  print('{0:}{1:}'.format(s,out))
#
# open file and read lines one by one
  f=open(joblog,"r")
  line =f.readline()
  while line :
    line =f.readline()

#==extract GEOMETRY ==
# check for standard orientation text set 'go' if the text appears
    if 'Standard orientation:' in line:
      go=True
      s='Structure in standard orientation found '
      if rep == 0 : print('{0:}'.format(s))
      n=1
      rep=rep+1
      count=0
# endif    
#
# if go is on then we need to extract and store the structure
    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))
        count=count+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

#now save to the store array nstore
#k is the count number of the logfile starting at 0
#we may have muliple geometries save only the last
#rep is the number of the last stored geometry
#count is the number of atoms in a geometry
  n=0
  if rep != 1 :
    n=len(col_num)-count
    s=' structures found only the last is saved'
    print('{0:}{1:}'.format(rep,s))
    if rep > 35 :
      s='** more than 35 structures opt may not be complete'
      print('{0:}'.format(s))
#   endif
# endif
  m=0
  while n < len(col_num):
#    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]))
    nstore.append([col_atom[n],col_x[n],col_y[n],col_z[n]])
#    print('{0:}'.format(nstore[k*count+m]))
    n=n+1
    m=m+1
#close while
#
#== build COM FILE ==
#
  o=open(out,"w")
#
# header information
  s="%nprocshared=32\n%mem=60000MB\n"
  o.write('{0:}'.format(s))
  s='%chk='+job2+'.chk\n'
  o.write('{0:}'.format(s))
#
# job information
  s="# freq b3lyp/6-31G(d,p) scrf=(smd,solvent=water) geom=cartesian int=ultrafine scf=conver=10 \n\n"
  o.write('{0:}'.format(s))
#
# geometry 
  s="Title Card Required\n\n-1 1\n"
  o.write('{0:}'.format(s))
  m=0
# structures are stored as atoms in molecule (total=count and counter=m) times the number of molecules (k)
  while m < count:
#  print('{0:}'.format(nstore[k*count+m]))
    o.write('{0:>4d} {1: 09.6f}  {2: 09.6f}  {3: 09.6f} \n'.format(nstore[k*count+m][0],nstore[k*count+m][1],nstore[k*count+m][2],nstore[k*count+m][3]))
    m=m+1
#close while
#
#  add the last blank lines
  s=" \n\n"
  o.write('{0:}'.format(s))
#
#close com file 
  f.close()
  o.close()
#
# update to next file
  k=k+1
#close while
#
sys.exit()