#Script for outputting energies for a bondlengthscan in Terachem by Alistair King

import sys
import matplotlib.pyplot as plt

#check for arguements

if len(sys.argv) != 2:
	print 'Usage: python bondlengthenergies.py <geometry.xyz> \n'
	print len(sys.argv)
	sys.exit(1)
	
#Open the file and read each line to a list, then close file

filename=sys.argv[1]

file=open(filename)

list = []

for line in file: 
        line = line.strip()
        list.append(line)
        
file.close()

#create two lists for the x and y values (bond length and energy)

xlist = []

ylist = []

#create a file for the comma separated values

output = open('scanenergiesxy.xy','w')

#search and extract the energies and bond lengths from the list
#then write them to file and the new lists..then close file

for item in list:
	index = item.find('Converged')
	while index == 0:
		index = item.find('Energy')
		print index
		text = item[index+34:index+42]
		output.write(text)
		output.write(',')
		ylist.append(text)
		print str(text)
		text = item[index+7:index+30]
		output.write(text)
		output.write('\n')
		xlist.append(text)
		print str(text)	
			
output.close()

#plot the x & y values...then save an image

plt.plot(ylist, xlist, 'ro-')
plt.ylabel('Energy (Hartree)')
plt.xlabel('Bond Length (Angstroms)')
plt.tight_layout()

plt.savefig('scanenergies.png')
plt.show()