#!@PYTHON@
#
# This code calculates the pair correlation between the water molecules and the center of gravity of the cluster (defined by the Si and Ca atoms)

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

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

# parse the pdb file and obtain maxid
print "Parsing the PDB file "+sys.argv[1]+"."
pdbfile = open(sys.argv[1], "r")
atoms = {}
maxid=-1
for line in pdbfile:
  if "END" in line or "CONECT" in line or "MASTER" in line:
    continue
  fields = line.split()
  if maxid < int(fields[1]):
    maxid=int(fields[1])
  # id, element, x, y, z and ten neighbour ids
  last = fields[10].split('0')
  atoms[fields[1]] = [int(fields[1]), last[3], float(fields[5]), float(fields[6]), float(fields[7]), -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]

print "The biggest Id found was %d" % (maxid)

# 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)

# go through all Si and Ca and calculate CoG
CoG = [ 0, 0, 0]
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):
  CoG[i] /= nr
print "The Center of Gravity is at %g,%g,%g." % (CoG[0], CoG[1], CoG[2])

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

# go through every atom
nrwaters = 0
nrothers = 0
for atom in atoms:
  if "O" in atoms[atom][1]:
    water = 1
    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]:
          water = 0
          break
      else:
        print "ERROR: key %d is not in list!" % (atoms[atom][i])
        sys.exit(1)
    if water == 1:
      dist = 0
      for i in range(3):
        dist += (CoG[i] - atoms[atom][2+i])*(CoG[i] - atoms[atom][2+i])
      dist = int(math.sqrt(dist)/BinWidth)
      if dist < NrBins:
        bins[dist] = bins[dist] + 1
        nrwaters = nrwaters+1
    else:
      nrothers = nrothers + 1

print "There were %d water molecules and %d oxygens connected to silica." % (nrwaters, nrothers)

#print "The Bins are as follows:"
output = open("histo_pair_correlation.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%g\n" % (i*BinWidth, bins[i]))
output.close()

print "finished."