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'
Given the energies could be put in the second line of the xyz input files. Any thoughts regarding the Hammond postulate, with the notion of giving you a weighted guess?
ReplyDeleteThis was a very silly script which ultimately didn't give me anything useful -- I had much more luck using QST2 in G09.
DeleteAlthough I don't show it in the 'usage', you can weight the output by using the -w switch, and I guess you can then use Hammond's postulate to weight it towards the products for endothermic reactions, and vice versa. Again, I ultimately had little use for the morphed structures and abandoned this approach.
Anyway, thanks for the feedback!