#Python script for determining proton affinities from 2 gamess output files by Alistair King
#Cartesian coordinates assumed for input files
#linear molecules with three or more atoms should be calculated by hand

import sys

#check for arguements

if len(sys.argv) != 5:
	print 'Usage: python ProtonAffinity.py <acidfile.out> <basefile.out> <Acid Geometry> <Base Geometry>\n'
	print 'Acid & Base Geometries M (Monoatomic), L (Linear), N (Non-linear)\n'
	print len(sys.argv)
	sys.exit(1)
	

Acidfilename=sys.argv[1]
Basefilename=sys.argv[2]

AcidGeom=str(sys.argv[3])
BaseGeom=str(sys.argv[4])

print AcidGeom
print BaseGeom


#some constants

H2kcal = 627.53 # for conversion from Hartrees to kcals
T = 298.15 # Kelvin
R = 0.0019872 # kcal.mol-1.K-1


#################acid#####################
#open the file, copy as a string and close

Acidfile=open(Acidfilename)

Afilestring=Acidfile.read()

Acidfile.close()


#rfind to locate the last ENERGY COMPONENTS in the string

index = Afilestring.rfind('ENERGY COMPONENTS', 1)

index = Afilestring.find('TOTAL ENERGY', index)

AElec = H2kcal*float(Afilestring[index+14:index+35])

print "\nThe electronic energy for the protonated species is " + str(AElec) + " kcal/mol"


#locate the THERMOCHEMISTRY section in the string

index = Afilestring.rfind('THE HARMONIC ZERO POINT ENERGY IS', 1)

if index == -1:
	AHZPE = 0

else:
	index = Afilestring.find('CM**-1/MOLECULE', index)
	AHZPE = float(Afilestring[index+20:index+33])

print "The harmonic zero-point energy for the protonated species is " + str(AHZPE) + " kcal/mol\n"


if AcidGeom == 'M':
	AErot = 0
	AEtrans = R*T*3/2
	print "Note: The Acid is assumed to be monoatomic!"
elif AcidGeom == 'N':
	AErot = R*T*3/2
	AEtrans = R*T*3/2
	print "Note: The Acid is assumed to be non-linear!"
elif AcidGeom == 'L':
	AErot = R*T
	AEtrans = R*T*3/2
	print "Note: The Acid is assumed to be linear!"
else:
    print "Please Check the input parameters for the Acid Geometry\n"

print "The Acids rotational energy is " + str(AErot) + " kcal/mol"
print "The Acids translational energy is " + str(AEtrans) + " kcal/mol\n\n"

#################base######################
#open the file, copy as a string and close

Basefile=open(Basefilename)

Bfilestring=Basefile.read()

Basefile.close()


#rfind to locate the last ENERGY COMPONENTS in the string

index = Bfilestring.rfind('ENERGY COMPONENTS', 1)

index = Bfilestring.find('TOTAL ENERGY', index)

BElec = H2kcal*float(Bfilestring[index+14:index+35])

print "\nThe electronic energy for the unprotonated species is " + str(BElec) + " kcal/mol"


#locate the THERMOCHEMISTRY section in the string

index = Bfilestring.rfind('THE HARMONIC ZERO POINT', 1)


if index == -1:
	BHZPE = 0

else:
	index = Bfilestring.find('CM**-1/MOLECULE', index)
	BHZPE = float(Bfilestring[index+20:index+33])

print "The harmonic zero-point energy for the unprotonated species is " + str(BHZPE) + " kcal/mol\n"

	
if BaseGeom == 'M':
	BErot = 0
	BEtrans = R*T*3/2
	print "Note: The Base is assumed to be monoatomic!"
elif BaseGeom == 'N':
	BErot = R*T*3/2
	BEtrans = R*T*3/2
	print "Note: The Base is assumed to be non-linear!"
elif BaseGeom == 'L':
	BErot = R*T
	BEtrans = R*T*3/2
	print "Note: The Base is assumed to be linear!"
else:
    print "Please Check the input parameters for the Base Geometry\n"
	

print "The Bases rotational energy is " + str(BErot) + " kcal/mol"
print "The Bases translational energy is " + str(BEtrans) + " kcal/mol\n\n\n"


###############final calculation################

AE = AErot+AEtrans+AElec+AHZPE

print "The acid total energy is " + str(AE) + " kcal/mol\n"

BE = BErot+BEtrans+BElec+BHZPE

print "The base total energy is " + str(BE) + " kcal/mol\n"

deltaE = AE - BE

PA =  R*T*5/2 - deltaE

PA*=-1

print "The proton affinity is " + str(PA) + " kcal/mol\n\n"