UsefullRandomScripts/extract_Stationary_structures.py at master · MarieTJ/UsefullRandomScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import os
import sys


if len(sys.argv) != 2 : sys.exit(' Usage : ExtractStationaryStructures.py prefix_  . With your files being prefix_angle.log')

prefix=sys.argv[1]
atomdict={'1':'H', '6':'C' ,'7':'N' , '8':'O', '16':'S' }
outputenergy=open('%s-energySP.dat' %(prefix) ,'w')
outputlist=open('%s-listfiles' %(prefix),'w')
for filestem in range(-180,180,1):
    file=prefix+str(filestem)+'.log'
    if os.path.isfile( file ):
        i=0
        filepointer=open(file, 'r')
        filecontent=filepointer.readlines()
        print(file)
        for linenumber in range(len(filecontent)):
          if '\\HF='  in  filecontent[linenumber]:

              lineenergy= (filecontent[linenumber]+filecontent[linenumber+1]).replace("\n", "").replace(" ", "").split('\\')
              k=0
              while lineenergy[k][:3]!= 'HF=': k+=1
              else :
                  if lineenergy[k][-2:]!='\n' : hfenergy=lineenergy[k]
                  else : hfenergy=filecontent[linenumber+1].split('\\')[0]
                  print( hfenergy)
        SP=False
        for linenumber in range(len(filecontent)):
            if 'Stationary point found' in  filecontent[linenumber]:
                SP=True
                print('Stationary point found in file %s' %(file))
                linecoord= linenumber
                while 'Center     Atomic      Atomic             Coordinates (Angstroms)' not in filecontent[linecoord]:
                    linecoord+=-1

                linecoord+=2
                fileout=open('SP_%s_%s.pdb' %(file[:-4],i),'w')
                #i+=1 ###Uncomment to save more than one stationanry point it to
                fileout.write('# Energy = %s \n' %(hfenergy))
                j=1

                while '----------------------------------------------------' not in filecontent[linecoord+j]:
                        linesplited=filecontent[linecoord+j].split()
                        #   print (filecontent[linecoord+j].split())
                        fileout.write('ATOM     %s   %s  AAA     1    %s%s%s  1.00  0.00 \n' %(linesplited[0].rjust(2),\
                                atomdict[linesplited[1]] , \
                                str("{0:.3f}".format(float(linesplited[3]))).rjust(8), \
                                str("{0:.3f}".format(float(linesplited[4]))).rjust(8), \
                                str("{0:.3f}".format(float(linesplited[5]))).rjust(8)))
                        j+=1
    #now write t  he energy of this structure
        if SP==True:
            print(hfenergy, hfenergy.split("HF=")[1], hfenergy.split("HF=")[0])
            outputenergy.write("%.8f\n" % float(hfenergy.split("HF=")[1]))
            outputlist.write('SP_%s_%s.pdb \n' %(file[:-4],i))