#!@PYTHON@
#
# converts mpqc hessian output files to pcp.hessian_{xx,xy,..,zz}.all for joiner/analyzer to work on

import sys, string, math
wrout = sys.stdout.write
wrerr = sys.stdout.write

if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <configbase> <prefix>\n")
  sys.exit(1)

#output = open(sys.argv[1]+".out", "r")
#config = open(sys.argv[1]+".conf", "r")

# 1. parse symmetrically stored data into list
input = open(sys.argv[1]+".hess", "r")
comment = input.readline()
atomline = input.readline()
atoms = string.split(atomline)
NoAtom = int(atoms[0])
for i in range(0,NoAtom):	# skip nuclear coordinates
  testline = input.readline()
hessianlist = []
NoValues=0
for line in input:
  if "End Hessian" in line:
    break
  else:
    values = string.split(line)
    for i in range(len(values)):
      NoValues = NoValues+1
      hessianlist.append(float(values[i]))
input.close()
print "I scanned "+str(NoValues)+" entries in the symmetric matrix for "+atoms[0]+" atoms. First value is "+str(hessianlist[0])+"."

# 2. Create full hessian matrix
hessianmatrix = []
NoValues=0
for i in range(3*NoAtom):	# row index
  for j in range(i+1):	# column index
    #print "["+str(3*NoAtom*i+j)+"] = ["+str(hessianlist[NoValues])+"]."
    hessianmatrix.append(hessianlist[NoValues])
    NoValues = NoValues+1
  for j in range(i+1,3*NoAtom):
    #print "["+str(3*NoAtom*i+j)+"] = [0]."
    hessianmatrix.append(0)
for i in range(3*NoAtom):	# row index
  for j in range(i+1,3*NoAtom):	# column index
    NoValues = NoValues+1
    index=j+3*NoAtom*i
    index2=i+3*NoAtom*j
    hessianmatrix[index]=hessianmatrix[index2]
print "I created the full hessian matrix from the symmetric one with a total of "+str(NoValues)+" entries."


# 3. store 9 (NoAtom x NoAtom)-matrices (xx,xy,...,zz)
textlist = [ "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" ]
for i in range(3):	# row index
  for j in range(3):	# column index
    wrout("Storing "+str(NoAtom)+"x"+str(NoAtom)+" matrix "+textlist[3*i+j]+" to file "+sys.argv[2]+".hessian_"+textlist[3*i+j]+".all ... ")
    output=open(sys.argv[2]+".hessian_"+textlist[3*i+j]+".all", "w")
    #print >>output,"# hessian matrix with "+str(NoAtom)+" by "+str(NoAtom)+" entries"
    output.write("# ")
    for m in range(NoAtom):	# row index
      output.write("Atom"+str(m+1)+"\t")
    output.write("\n")
    for m in range(NoAtom):	# row index
      for n in range(NoAtom): # column index
        #index = (3*NoAtom)*((NoAtom*i)+m) + (NoAtom*j)+n
        index = (3*NoAtom)*(i+(3*m)) + (j+(3*n))
        output.write("%f\t" % (hessianmatrix[index]))
      output.write("\n")
    wrout("done.\n")

print "Everything is converted."