#!@PYTHON@
#
# This code calculates the pair correlation between the water molecules and the non-convex surface of the cluster (i.e. the boundary defined by the Si and Ca atoms)

import sys, string, re, math, os
wrerr = sys.stderr.write
wrout = sys.stdout.write

if (len(sys.argv) < 4):
  print "Usage: "+sys.argv[0]+" <xyz> <dbond> <dat>"
  sys.exit(1)

# parse the pdb file and obtain maxid
print "Parsing the XYZ file "+sys.argv[1]+"."
xyzfile = open(sys.argv[1], "r")
atoms = {}
id=1
xyzfile.readline()
xyzfile.readline()
emptyline=re.compile('[ \t]*\n')
for line in xyzfile:
  #if "END" in line or "CONECT" in line or "MASTER" in line:
  if emptyline.match(line):
    continue
  fields = line.split()
  # id, element, x, y, z and ten neighbour ids
  atoms[str(id)] = [id, fields[0], float(fields[1]), float(fields[2]), float(fields[3]), -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
  id+=1
print "The biggest Id is %d" % (id)
xyzfile.close()

# Parsing the dbond file
print "Parsing the DBOND file "+sys.argv[2]+"."
dbondfile = open(sys.argv[2], "r")
test = dbondfile.readline()	# skip the first two lines
test = dbondfile.readline()
for line in dbondfile:
  fields = line.split()
  ids = [int(fields[0]), int(fields[1])]
  set = [0,0]
  if str(ids[0]) in atoms and str(ids[1]) in atoms:
    for i in range(5,15):
      for j in range(2):
        if set[j] == 0 and (atoms[ str(ids[j]) ])[i] == -1:
         (atoms[ str(ids[j]) ])[i] = ids [ (j+1) %2 ]
         set[j] = 1
      if set[0] == 1 and set[1] == 1:
        break
    if set[0] == 0 or set[1] == 0:
      print "ERROR: Either id already had ten neighbours, increase number in code!"
      sys.exit(1)
  else:
    print "ERROR: Either %d or %d are not listed in the PDB file!" % (ids[0], ids[1])
    sys.exit(1)

# parse the tecplot dat file (the first three lines are tecplot INFO)
datfile=open(sys.argv[3], "r")
line=datfile.readline()
line=datfile.readline()
line=datfile.readline()
points=[]
pcounter=0
for line in datfile:
  if not line.strip():	# end on first empty line, indicating end of points and beginning of triangles
    break
  else:
    fields=line.split()
    points.append([pcounter, float(fields[0]), float(fields[1]), float(fields[2])])
    pcounter=pcounter+1
print "There are %d endpoints of triangles." % (pcounter)

tcounter=0
triangles=[]
for line in datfile:
  fields=line.split()
  triangles.append([tcounter, int(fields[0]), int(fields[1]), int(fields[2])])
  corners=[points[ int(fields[0])-1 ], points[ int(fields[1])-1 ], points[ int(fields[2])-1 ]]
  # calculate middle point and add to points
  centroid=[0.,0.,0.]
  for i in range(3):
    for j in range(3):
      centroid[j] += corners[i][1+j]
  for i in range(3):
    centroid[i] = centroid[i]/3.
  points.append([pcounter+tcounter, centroid[0], centroid[1], centroid[2]])
  #print ("Centroid is at (%g,%g,%g) for triangle at (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)." % (centroid[0], centroid[1], centroid[2], corners[0][1],corners[0][2],corners[0][3], corners[1][1],corners[1][2],corners[1][3], corners[2][1],corners[2][2],corners[2][3]))
  
  tcounter=tcounter+1
print "There are %d triangles." % (tcounter)

# go through all Si and Ca and calculate CoG
CoG = [ 0, 0, 0]
key=atoms.keys()[0]
max=[atoms[key][2],atoms[key][3],atoms[key][4]]
min=[atoms[key][2],atoms[key][3],atoms[key][4]]
nr=0
for atom in atoms:
  #print atoms[atom]
  if "Si" in atoms[atom][1] or "Ca" in atoms[atom][1]:
    for i in range(3):
      CoG[i] += atoms[atom][2+i]
    nr = nr + 1
    for i in range(3):
      if min[i] > atoms[atom][2+i]:
        min[i] = atoms[atom][2+i]
      if max[i] < atoms[atom][2+i]:
        max[i] = atoms[atom][2+i]
for i in range(3):
  CoG[i] /= nr
print "The Center of Gravity is at %g,%g,%g." % (CoG[0], CoG[1], CoG[2])
print "Bounds of the cluster's atoms are [%f,%f,], [%f,%f], [%f,%f]." % (min[0],max[0], min[1],max[1], min[2],max[2])

# check bounds against bounds of cluster
max=[points[0][1],points[0][2],points[0][3]]
min=[points[0][1],points[0][2],points[0][3]]
for j in range(len(points)):
  for i in range(3):
    if min[i] > points[j][1+i]:
      min[i] = points[j][1+i]
    if max[i] < points[j][1+i]:
      max[i] = points[j][1+i]
print "Bounds of the cluster's tesselated surface are [%f,%f,], [%f,%f], [%f,%f]." % (min[0],max[0], min[1],max[1], min[2],max[2])

# find radius of cluster
bounds=[0.,0.,0.]
maximum=0
maxid=-1
for i in range(3):
  bounds[i] = max[i] - min[i]
  if maximum < bounds[i]:
    maxid=i
    maximum = bounds[i]
print "Cylinder direction is probably "+str(maxid)+"."
R=(bounds[(maxid+1)%3]+bounds[(maxid+2)%3])/4.

# search atom closest to CoG to allow visual check
mindist=1e+16
minid=-1
for atom in atoms:
  dist=0
  for i in range(3):
    temp = CoG[i] - atoms[atom][2+i]
    dist += temp*temp
  if dist < mindist:
    mindist = dist
    minid = atom
print "The atom closest to CoG is %d, element %s (%f,%f,%f) with distance %f." % (atoms[minid][0], atoms[minid][1], atoms[minid][2], atoms[minid][3], atoms[minid][4], math.sqrt(mindist))

# initialize the bin array
NrBins=100
BinWidth=0.5
bins = []
for i in range(2*NrBins+1):
  bins.append(0)

# get path
line=sys.argv[1].rpartition('/')
path=line[0]+line[1]
print "Path is "+path+"."

# open a number of xyz files
innerwaterfile=open(path+"innerwater","w")
outerwaterfile=open(path+"outerwater","w")
hyperdryclusterfile=open(path+"hyperdrycluster","w")
dryclusterfile=open(path+"drycluster","w")
waterdryclusterfile=open(path+"waterdrycluster","w")
layerclusterfile=open(path+"layercluster","w")
outerwaterfile=open(path+"outerwater","w")
nrcluster=0
for atom in atoms:
  if "Si" in atoms[atom][1] or "Ca" in atoms[atom][1]:
    #innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    #outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    hyperdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
    nrcluster+=1
nohyperdry=nrcluster
nodry=nrcluster
nowaterdry=nrcluster
nolayer=nrcluster

# go through every atom
elements=["H", "O"]
nrwaters = 0
nrothers = 0
inner=0
outer=0
for element in elements:
  print "Calculating for element "+element+"."
  for i in range(2*NrBins+1):
    bins[i] = 0.
  for atom in atoms:
    if element in atoms[atom][1]:
      water = 1
      # check whether its really water (i.e. not bound to Si)
      for i in range(5,15):
        if atoms[atom][i] == -1:
          break
        if str(atoms[atom][i]) in atoms:      
          if "Si" in atoms[ str(atoms[atom][i]) ][1] or "Ca" in atoms[ str(atoms[atom][i]) ][1]:
            water = 0
            break
        else:
          print "ERROR: key %d is not in list!" % (atoms[atom][i])
          sys.exit(1)
      # if it's water find closest boundary point and check whether atom is in- or outside of cluster
      if water == 1:
        mindist=1e+16
        vector=[0.,0.,0.]
        vector2=[0.,0.,0.]
        minptr=points[0]
        for ptr in points:
          dist = 0
          for i in range(3):
            vector[i] = ptr[1+i] - atoms[atom][2+i]
            dist += (vector[i])*(vector[i])
          if dist < mindist:
            minptr=ptr
            mindist = dist
            sign=1
            skp=0
            diff=0
            for i in range(3):
              vector2[i] = CoG[i] - atoms[atom][2+i]
              skp += vector[i]*vector2[i]
              diff += vector2[i]*vector2[i]
            if (skp < 0) or (diff < dist):
              sign=-1
        mindist = int(math.sqrt(mindist)/BinWidth)
        if mindist < NrBins:
          #if mindist < 8 and abs(sign) == 1:
            #for i in range(3):
              #vector[i] = minptr[1+i] - atoms[atom][2+i]
              #vector2[i] = CoG[i] - atoms[atom][2+i]
            #print "For atom %d vector to boundary is (%f,%f,%f) and vector to CoG (%f,%f,%f)." % (atoms[atom][0], vector[0], vector[1], vector[2], vector2[0], vector2[1], vector2[2])
            #print "atom %d vector is (%f,%f,%f) and boundary point vector is (%f,%f,%f)." % (atoms[atom][0], atoms[atom][2], atoms[atom][3], atoms[atom][4],minptr[1], minptr[2], minptr[3])
            #print "Atom %d goes into bin %d." % (atoms[atom][0], mindist)
          if sign == -1:
            bins[NrBins-mindist] = bins[NrBins-mindist] + 1
            innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
            dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
            waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
            layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
            nodry += 1
            nowaterdry += 1
            nolayer += 1
            inner += 1
          else:
            bins[NrBins+mindist] = bins[NrBins+mindist] + 1
            outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
            if mindist < 8:
              waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
              nowaterdry += 1
            if mindist < 11:
              layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
              nolayer += 1
            outer += 1
          nrwaters = nrwaters+1
        else:
          print "NOTE: Have to throw "+str(atoms[atom][0])+" away, as it would be put into Bin "+str(mindist)+" which is out of range 0,"+str(NrBins-1)+"."
      else:
        if element == elements[1]:
          nrothers = nrothers + 1
          nrcluster = nrcluster + 1
          #innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          #outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          hyperdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
          nodry += 1
          nohyperdry += 1
          nowaterdry += 1
          nolayer += 1
  
  #print "The Bins are as follows:"
  output = open(path+"histo_pair_correlation_"+element+".dat", "w")
  output.write("#Distance [A]\tnumber in bin\n")
  for i in range(NrBins):
    #print "%g\t%d" % (i*BinWidth, bins[i])
    output.write("%g\t%d\n" % (-(NrBins-i)*BinWidth, bins[i]))
  for i in range(NrBins,2*NrBins):
    #print "%g\t%d" % (i*BinWidth, bins[i])
    output.write("%g\t%d\n" % ((i-NrBins)*BinWidth, bins[i]))
  output.close()
  
  
  # calculate numerically the bins for homogeneous water distribution where we assume the cluster to be a cylinder (hence, just 2d and circle)
  # centered in the box of length 44 with a certain radius R
  a=22.0 # 44/2
  N=2200
  lastarea=0
  waterdensity=0.033456344	# (1g/18amu)/cm^3 -> 1/angstrom^3, i.e. number of hydrogen particles in Angstroem^3 volume
  if "H" in element:
    waterdensity = waterdensity*2
  output=open(path+"homogeneous_water_histo_"+element+".dat","w")
  for i in range(1,2*NrBins):
    area=0
    radius = float(i)*BinWidth
    h=(a)/float(N)
    for j in range(N):
      x = j*h
      if x < a and x < radius:
        y = (radius)*(radius) - x*x
        if y > a*a:
          y=a*a
        dx = h
        if (x+dx > a):
          dx -= a-x
        if (x+dx > radius):
          dx -= radius-x
        if (dx < 0):
          print "ERROR: dx is "+str(dx)+"."
        area += dx * math.sqrt(y)
    #print ("Area for radius %g is %g." % (radius, area-lastarea))
    volume = int((2.*a) * 4.* (area - lastarea) * waterdensity)	# 2a is length of cylinder, 4* due to only looking at a quadrant
    lastarea = area
    if radius <= R:
      volume = 0
    output.write("%g\t%d\n"  % ((i*BinWidth)-R, volume))
  output.close()

print "There were "+str(nrwaters)+" hydrogen and oxygen atoms and "+str(nrothers)+" connected to silica."
print "There were "+str(inner)+" inner and "+str(outer)+" outer hydrogen and oxygen atoms."

innerwaterfile.close()
outerwaterfile.close()
hyperdryclusterfile.close()
dryclusterfile.close()
waterdryclusterfile.close()
layerclusterfile.close()

# put number of atoms in front
innerwaterfile=open(path+"innerwater.xyz","w")
innerwaterfileold=open(path+"innerwater","r")
innerwaterfile.write(str(inner)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in innerwaterfileold:
  innerwaterfile.write(line)
innerwaterfile.close()
innerwaterfileold.close()
os.remove(path+"innerwater")

outerwaterfile=open(path+"outerwater.xyz","w")
outerwaterfileold=open(path+"outerwater","r")
outerwaterfile.write(str(outer)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in outerwaterfileold:
  outerwaterfile.write(line)
outerwaterfile.close()
outerwaterfileold.close()
os.remove(path+"outerwater")

hyperdryclusterfile=open(path+"hyperdrycluster.xyz","w")
hyperdryclusterfileold=open(path+"hyperdrycluster","r")
hyperdryclusterfile.write(str(nohyperdry)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in hyperdryclusterfileold:
  hyperdryclusterfile.write(line)
hyperdryclusterfile.close()
hyperdryclusterfileold.close()
os.remove(path+"hyperdrycluster")

dryclusterfile=open(path+"drycluster.xyz","w")
dryclusterfileold=open(path+"drycluster","r")
dryclusterfile.write(str(nodry)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in dryclusterfileold:
  dryclusterfile.write(line)
dryclusterfile.close()
dryclusterfileold.close()
os.remove(path+"drycluster")

waterdryclusterfile=open(path+"waterdrycluster.xyz","w")
waterdryclusterfileold=open(path+"waterdrycluster","r")
waterdryclusterfile.write(str(nowaterdry)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in waterdryclusterfileold:
  waterdryclusterfile.write(line)
waterdryclusterfile.close()
waterdryclusterfileold.close()
os.remove(path+"waterdrycluster")

layerclusterfile=open(path+"layercluster.xyz","w")
layerclusterfileold=open(path+"layercluster","r")
layerclusterfile.write(str(nolayer)+"\n\tCreated by "+sys.argv[0]+"\n")
for line in layerclusterfileold:
  layerclusterfile.write(line)
layerclusterfile.close()
layerclusterfileold.close()
os.remove(path+"layercluster")

print "finished."