Talk:Mod:Hunt Research Group/extract crest conf

#!/opt/local/bin/python3
#
# extract_crest_conf.py
#
# python script to extract and print single structure from a trajectory of structures given by CREST
# to run the script type python extract_conf.py conf_no 
# the input file is crest_conformers.xyz
# conf_no is the structure to extract
# output file will be frame_x.xyz
#
import sys
import os
dir=os.getcwd()
#
# get input
#
if len(sys.argv) == 1:
  print('to run the script please type: python extract_crest_conf.py conf_no')
  sys.exit()
else:
  conf=str(sys.argv[1])
  conf_no=int(sys.argv[1])
  base='conf_'
  file='crest_rotamers_sorted.xyz'
  outfile=base+conf+'.xyz'
  s='directory is '
  print('{0:}{1:}'.format(s,dir))
  s='input file is '
  print('{0:}{1:}'.format(s,file))
  s='output file is '
  print('{0:}{1:}'.format(s,outfile))
#close if
#
# set some defaults
#
col_num=[]
col_atom=[]
col_energy=[]
col_x=[]
col_y=[]
col_z=[]
#
# determine file length
#
cmd='wc crest_rotamers_sorted.xyz'
out=os.popen(cmd).read()
# s='system call '
# print('{0:}'.format(out))
s1,s2,s3,s4=out.rstrip().split()
# s='returns '
# print('{0:} {1:} {2:} {3:} {4:}'.format(s,s1,s2,s3,s4))
file_length=int(s1)
s='file length '
#print('{0:}{1:}'.format(s,file_length))
#
# find start and stop points
#
f=open(file,"r")
line =f.readline()
no_atoms=int(line.rstrip())
s='number of atoms '
print('{0:}{1:}'.format(s,no_atoms))
block_length=no_atoms+2
start=1+(conf_no-1)*block_length
finish=conf_no*block_length
tot_conf_no=int(file_length/block_length)
s='number of conformers '
print('{0:}{1:}'.format(s,tot_conf_no))
s='printing conformer '
print('{0:}{1:}'.format(s,conf_no))
s1='start at line '
s2='  finish at line '
#print('{0:}{1:}{2:}{3:}'.format(s1,start,s2,finish))
#
# read lines and store data           
#
# i=conformer_pointer
# j=atom_pointer
j=1
while line:
#  print(j)
  line =f.readline()
#  print(line)
  store=[]
  store=line.rstrip().split()
  e1=float(store[0])
  col_energy.append(e1)
#  print(e1)
  i=0
  while i < no_atoms:
    line =f.readline()
#    print(line)
    a1,b1,c1,d1=line.rstrip().split()
    col_atom.append(a1)
    col_x.append(float(b1))
    col_y.append(float(c1))
    col_z.append(float(d1))
    i=i+1
# close while
  line =f.readline()
  j=j+1
#close while

#
#now write to the .xyz file
#
c=open(outfile,"w")
print('{0:}'.format(no_atoms))
c.write('{0:} \n'.format(no_atoms))
print('{0:}'.format(col_energy[conf_no]))
c.write('{0:} \n'.format(col_energy[conf_no]))
#
start=(conf_no-1)*no_atoms
finish=conf_no*no_atoms
s1='start at line '
s2='  finish at line '
#print('{0:}{1:}{2:}{3:}'.format(s1,start,s2,finish))
n=start
while n < finish:
  c.write('{0:3s}  {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:3s}  {1: 09.6f}  {2: 09.6f}  {3: 09.6f}'.format(col_atom[n],col_x[n],col_y[n],col_z[n]))
  n=n + 1
#close while

# close all open files
f.close()
c.close()
exit(0)
#