#!/usr/bin/env python
"""
Usage
efgparse_g09 gaussian.out efgs.out
where gaussian.out is the output from the G09 run and efgs.out is the destination file for the processed data.
"""
import sys

def periodictable(elementnumber):
	ptable={1:'H',2:'He',\
	3:'Li', 4:'Be',5:'B',6:'C',7:'N',8:'O',9:'F',10:'Ne',\
	11:'Na',12:'Mg',13:'Al',14:'Si',15:'P',16:'S',17:'Cl',18:'Ar',\
	19:'K',20:'Ca',\
	21:'Sc',22:'Ti',23:'V',24:'Cr',25:'Mn',26:'Fe',27:'Co',28:'Ni',29:'Cu',30:'Zn',\
	31:'Ga',32:'Ge',33:'As',34:'Se',35:'Br',36:'Kr',\
	37:'Rb',38:'Sr',\
	39:'Y',40:'Zr',41:'Nb',42:'Mo',43:'Tc',44:'Ru',45:'Rh',46:'Pd',47:'Ag',48:'Cd',\
	49:'In',50:'Sn',51:'Sb',52:'Te',53:'I',54:'Xe',\
	55:'Cs',56:'Ba',\
	57:'La',58:'Ce',59:'Pr',60:'Nd',61:'Pm',62:'Sm',63:'Eu',64:'Gd',65:'Tb',66:'Dy',67:'Ho',68:'Er',69:'Tm',70:'Yb',71:'Lu',\
	72:'Hf', 73:'Ta', 74:'W',75:'Re', 76:'Os', 77:'Ir',78:'Pt', 79:'Au', 80:'Hg',\
	81:'Tl', 82:'Pb', 83:'Bi',84:'Po',85:'At',86:'Rn',\
	87:'Fr',88:'Ra',\
	89:'Ac',90:'Th',91:'Pa',92:'U',93:'Np',94:'Pu',95:'Am',96:'Cm',97:'Bk',98:'Cf',99:'Es',100:'Fm',101:'Md',102:'No',\
	103:'Lr',104:'Rf',105:'Db',106:'Sg',107:'Bh',108:'Hs',109:'Mt',110:'Ds',111:'Rg',112:'Cn',\
	113:'Uut',114:'Fl',115:'Uup',116:'Lv',117:'Uus',118:'Uuo'}
	element=ptable[elementnumber]
	return element

def barns(elementnumber):
	# see periodictable() function above for quick conversion between atomic number and element
	#http://www.easyspin.org/documentation/isotopetable.html
	#Using the quadrupole moment of the most abundant stable isotope which possesses a quadrupole moment
	#I've used 2H, 7Li, 11B, 14N, 17O, 21Ne, 23Na, 25Mg, 33S, 35Cl, 39K, 47Ti, 51V, 53Cr, 63Cu
	# 69Ga, 73Ge, 79Br, 83Kr, 85Rb, 87Sr, 92Zr, 93Nb, 95Mo, 101Ru, 105Pd, 115In, 121Sb, 127I, 131Xe
	# 137Ba, 139La, 145Nd, 147Sm, 153Eu, 157Gd, 163Dy, 167Er, 173Yb, 177Hf, 181Ta, 187Re, 189Os
	# 193Ir, 201Hg, 235U
	btable={1:0.002860,2:0,\
	3:-0.0401, 4:0.05288,5:0.04059,6:0,7:0.02044,8:-0.02558,9:0,10:0.10155,\
	11:0.104,12:0.1994,13:0.1466,14:0,15:0,16:-0.0678,17:-0.08165,18:0,\
	19:0.0585,20:-0.0408,\
	21:-0.220,22:-0.302,23:-0.052,24:-0.150,25:0.330,26:0,27:0.420,28:0.162,29:-0.220,30:'Zn',\
	31:0.171,32:-0.196,33:0.314,34:0,35:0.313,36:0.259,\
	37:0.276,38:0.335,\
	39:0,40:-0.176,41:-0.320,42:-0.022,43:-0.129,44:0.457,45:0,46:0.660,47:0,48:0,\
	49:0.810,50:0,51:-0.360,52:0,53:-0.710,54:-0.114,\
	55:-0.00343,56:0.245,\
	57:0.200,58:0,59:-0.0589,60:-0.330,61:0.740,62:-0.259,63:2.412,64:1.350,65:1.432,66:2.648,67:3.580,68:3.565,69:0,70:2.800,71:3.56,\
	72:3.365, 73:3.170, 74:0,75:2.070, 76:0.856, 77:0.751,78:0, 79:0.547, 80:0.386,\
	81:0, 82:0, 83:-0.516,84:0,85:0,86:0,\
	87:0,88:0,\
	89:1.7,90:0,91:0,92:4.936,93:0,94:0,95:0,96:0,97:0,98:0,99:0,100:0,101:0,102:0,\
	103:0,104:0,105:0,106:0,107:0,108:0,109:0,110:0,111:0,112:0,\
	113:0,114:0,115:0,116:0,117:0,118:0}
	barn=btable[elementnumber]
	return barn	

def parseout(infile):
	f=open(infile,'r')

	#check for l202 and l602 in output
	l202=0 #l202 Reorients coordinates, calculates symmetry, and checks variables
	l602=0 #l602 1-electron properties (potential, field, and field gradient)
	efg=0
	collect=0
	structure=[]
	efgs=[]
	for line in f:
		if l202==1 and collect==2 and not ("---" in line):
			structure+=[line.rstrip()]
		if l602==1 and efg==1 and not  ("---" in line):
			efgs+=[line.rstrip()]
		if 'l202' in line:
			if l202==0:
				l202=1
				structure=[]
		if 'l602' in line:
			if l602==0:
				l602=1
				efgs=[]
		if l602==1 and 'tensor representation' in line:
			collect=1
		if l602==1 and collect==1 and 'Eigenvalues' in line:
			collect=2
		if l202==1 and "---------------------------------------------------------------------" in line:
			if collect<2: #found start ---, start collecting structure
				collect+=1
			elif collect==2: #found end ---, stop collecting structure
				collect=0
				l202=0
		if collect==2 and l602==1 and "-----------------------------------------------------" in line:
			if efg==0: #found start ---, start collecting structure
				efg=1
			elif efg==1: #found end ---, stop collecting structure
				efg=0
				l602=0
				collect=0
	return structure, efgs

def genxyzform(coords):
	x_str='%10.5f'% coords[0]
	y_str='%10.5f'% coords[1]
	z_str='%10.5f'% coords[2]
	xyz_form=[x_str, y_str,z_str]
 	return xyz_form

def genxyzstring(coords,elementnumber):
	x_str='%10.5f'% coords[0]
	y_str='%10.5f'% coords[1]
	z_str='%10.5f'% coords[2]
	element=periodictable(int(elementnumber))
	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 getstructures(rawdata,outputfile):
	n=1
	num="%03d" % (n,)
	g=open(outputfile+num+'.xyz','w')
	cartesian=[]
	xyz=[]
	for item in structure:
		coords=filter(None,item.split(' '))
		coordinates=[float(coords[3]),float(coords[4]),float(coords[5])]
		element=coords[1]
		xyz_form=genxyzform(coordinates)
		xyz+=[[periodictable(int(element)), xyz_form[0],xyz_form[1], xyz_form[2],barns(int(element))]]

		cartesian+=[genxyzstring(coordinates,element)]
	g.write(str(len(cartesian))+'\n')
	g.write('Structure '+str(n)+' '+'\n')
	for line in cartesian:
		g.write(line)
	g.close()
	return xyz

def genefg(coords):
	x_str='%10.5f'% coords[0]
	y_str='%10.5f'% coords[1]
	z_str='%10.5f'% coords[2]
	efg=[x_str, y_str,z_str]
	return efg


def formatefgs(efg):
	efgs=[]
	for item in efg:
		coords=filter(None,item.split(' '))
		coordinates=[float(coords[2]),float(coords[3]),float(coords[4])]
		efgs+=[coordinates]
	return efgs

def concatenate(xyz,efgs):
	catd=[]
	for n in range(0,len(xyz)):
		xyz[n].extend(efgs[n])
		catd+=[xyz[n]]
	return catd

def calculate(xyz):
	electron=1.5569E-25
	h=6.6260695729E-34
	
	for n in range(0,len(xyz)):
		v=[xyz[n][5],xyz[n][6],xyz[n][7]]
		tensors=map(abs,v)
		tensors.sort()
		if tensors[2] not in v:
			vzz=-tensors[2]
		else:
			vzz=tensors[2]
			
		Q=xyz[n][4]
		
		eta='%10.5f'% ((tensors[1]-tensors[0])/tensors[2])
		qcc='%10.5f'% (((electron*Q/h)*-vzz)/1.0E6)
		
		xyz[n].extend([vzz,eta,qcc])
	return xyz

def quickformat(floatvalue):
	stringvalue='%10.5f'% floatvalue
	return stringvalue

def formatoutput(xyz,outputfile):
	n=1
	num="%03d" % (n,)
	all_data=[]
	
	g=open(outputfile+num+'.data','w')

	xyz_formatted=""
	for n in range(0,len(xyz)):
		for m in range(4,len(xyz[n])-2):
			xyz[n][m]=quickformat(xyz[n][m])

		xyz_string=xyz[n][0]+(3-len(xyz[n][0]))*' '+10*' '+\
		(8-len(xyz[n][1]))*' '+xyz[n][1]+10*' '+(8-len(xyz[n][2]))*' '+xyz[n][2]+6*' '+(8-len(xyz[n][3]))*' '+xyz[n][3]+6*' '+\
		(8-len(xyz[n][4]))*' '+xyz[n][4]+6*' '++(8-len(xyz[n][5]))*' '+xyz[n][5]+6*' '++(8-len(xyz[n][6]))*' '+xyz[n][6]+6*' '+(8-len(xyz[n][7]))*' '+\
		xyz[n][7]+6*' '+(8-len(xyz[n][8]))*' '+xyz[n][8]+6*' '+(8-len(xyz[n][9]))*' '+xyz[n][9]+6*' '+(8-len(xyz[n][10]))*' '+str('%10.5f'% float(xyz[n][10]))+'\n'
		all_data+=[xyz_string]
	
	header='Atom'+(3-len('Atom'))*' '+10*' '+\
		(9-len('X'))*' '+'X'+10*' '+(10-len('Y'))*' '+'Y'+6*' '+(10-len('Z'))*' '+'Z'+\
		+6*' '+(10-len('Q'))*' '+'Q'+6*' '+(10-len('Tensor 1'))*' '+'Tensor 1'+6*' '+(10-len('Tensor 2'))*' '+'Tensor 2'+6*' '+(10-len('Tensor 3'))*' '+\
		'Tensor 3'+6*' '+(10-len('v_zz'))*' '+'v_zz'+6*' '+(10-len('eta'))*' '+'eta'+6*' '+(4-len('QCC'))*' '+'QCC (MHz)'+'\n'
	#g.write('Atom X Y Z Q XX YY ZZ eta QCC')
	print header,
	g.write(header)
	for line in all_data:
		g.write(line)
		print line,
	g.close()
	return all_data

if __name__ == "__main__":
	infile=sys.argv[1]
	outfile=sys.argv[2]

	structure,efg=parseout(infile)
	xyz=getstructures(structure,outfile)
	efgs=formatefgs(efg)
	xyz=concatenate(xyz,efgs)
	xyz=calculate(xyz)
	to_output=formatoutput(xyz,outfile)