466. morph xyz -- python script to morph .xyz files

Rather naively I was hoping that by comparing two  molecule .xyz files and generating an average of them I would be able to conveniently generate a half-decent transition state guess.

Turns out that it's not quite as simple. However, I've written the software, so I might as well share it.

Note that it's written in python 2.7 (i.e. not python 3)

Run the script without arguments for help. General usage is

morphxyz -i 1.xyz 2.xyz -o morph.xyz

So here it is:

morphxyz:

#!/usr/bin/python

import sys

def getvars(arguments):
 switches={}

 version='0.1'
 
 try:
  if "-i" in arguments:
   switches['in_one']=arguments[arguments.index('-i')+1]
   switches['in_two']=arguments[arguments.index('-i')+2]
   print 'Input: %s and %s'% (switches['in_one'],switches['in_two'])
  else:
   arguments="--help";
 except:
  arguments="--help";
  
 try:
  if "-o" in arguments:
   switches['o']=arguments[arguments.index('-o')+1].lower()
   print 'Output: %s'% switches['o']
  else:
   arguments="--help";
 except:
  arguments="--help";

 try:
  if "-w" in arguments:
   switches['w']=float(arguments[arguments.index('-w')+1])
   print 'Weighting: %i'% switches['w']
  else:
   print 'Assuming no weighting'
   switches['w']=1.0;
 except:
  switches['w']=1.0;

 doexit=0
 try:
  if ("-h" in arguments) or ("--help" in arguments):
   print '\t\t bytes2words version %s' % version
   print '\t-i\t two xyz files to morph'
   print '\t-o\t output file'
   print '\t-w\t weight one structure vs the other (1=average; 0=start; 2=end)'
   print 'Exiting'
   doexit=1
 except:
  a=0 # do nothing
 if doexit==1:
  sys.exit(0)

 return switches

def getcmpds(switches):
 
 cmpds={}
 
 g=open(switches['in_one'],'r') 
 n=0
 xyz=[]
 atoms=[]
 
 for line in g:
  n+=1
  line=line.rstrip('\n')
  if n==1:
   cmpds['atoms_one']=int(line)
  elif n==2:
   cmpds['title_one']=line
  else:
   line=line.split(' ')
   line=filter(None,line)
   xyz+=[[float(line[1]),float(line[2]),float(line[3])]]
   atoms+=[line[0].capitalize()]
 cmpds['coords_one']=xyz
 cmpds['elements_one']=atoms
 
 g.close
 
 g=open(switches['in_two'],'r') 
 n=0
 xyz=[]
 atoms=[]
 
 for line in g:
  n+=1
  line=line.rstrip('\n')
  if n==1:
   cmpds['atoms_two']=int(line)
  elif n==2:
   cmpds['title_two']=line
  else:
   line=line.split(' ')
   line=filter(None,line)
   xyz+=[[float(line[1]),float(line[2]),float(line[3])]]
   atoms+=[line[0].capitalize()]
 cmpds['coords_two']=xyz
 cmpds['elements_two']=atoms
 g.close
 
 cmpds['w']=switches['w']
 
 return cmpds

def morph(cmpds):
 coords_one=cmpds['coords_one']
 coords_two=cmpds['coords_two']
 
 coords_morph=[]
 coords_diff=[]
 for n in range(0,cmpds['atoms_one']):
  morph_x=coords_one[n][0]+cmpds['w']*(coords_two[n][0]-coords_one[n][0])/2.0
  morph_y=coords_one[n][1]+cmpds['w']*(coords_two[n][1]-coords_one[n][1])/2.0
  morph_z=coords_one[n][2]+cmpds['w']*(coords_two[n][2]-coords_one[n][2])/2.0
  diff_x=coords_two[n][0]-coords_one[n][0]
  diff_y=coords_two[n][1]-coords_one[n][1]
  diff_z=coords_two[n][2]-coords_one[n][2]
  coords_morph+=[[morph_x,morph_y,morph_z]]
  coords_diff+=[[diff_x,diff_y,diff_z]]
 cmpds['coords_morph']=coords_morph
 cmpds['coords_diff']=coords_diff
 return cmpds

def genxyzstring(coords,element):
 x_str='%10.5f'% coords[0]
 y_str='%10.5f'% coords[1]
 z_str='%10.5f'% coords[2]
 
 xyz_string=element+(3-len(element))*' '+10*' '+\
 (8-len(x_str))*' '+x_str+10*' '+(8-len(y_str))*' '+y_str+10*' '+(8-len(z_str))*' '+z_str+'\n'
 
 return xyz_string

def writemorph(cmpds,outfile):
 g=open(outfile,'w') 
 h=open('diff.xyz','w')
 g.write(str(cmpds['atoms_one'])+'\n'+'\n')
 h.write(str(cmpds['atoms_one'])+'\n'+'\n')
 
 for n in range(0,cmpds['atoms_one']):
  coords=cmpds['coords_morph'][n]
  diffcoords=cmpds['coords_diff'][n]
  
  g.write(genxyzstring(coords, cmpds['elements_one'][n]))
  h.write(genxyzstring(diffcoords, cmpds['elements_one'][n]))
    
 g.close
 h.close
 return 0

if __name__=="__main__":
 arguments=sys.argv[1:len(sys.argv)]
 switches=getvars(arguments)
 cmpds=getcmpds(switches)
 
 if cmpds['atoms_one']!=cmpds['atoms_two']:
  print 'The number of atoms differ. Exiting'
  sys.exit(1)
 elif cmpds['elements_one']!=cmpds['elements_two']:
  print 'The types of atoms differ. Exiting'
  sys.exit(1)
  
 cmpds=morph(cmpds)
 success=writemorph(cmpds,switches['o'])
 if success==0:
  print 'Conversion seems successful'