#!@PYTHON@
#
# if a given block from a PDB was cut out, water molecules may have been cut through. This is bad for the neighbour numbers in the TIP3 potential.
# Hence, this script removes all dangling atoms.

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

# check for parameters
if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <datafile> <offset to SeqID>\n")
  sys.exit(1)

input = open(sys.argv[1], "r")
output = open(sys.argv[1]+".new", "w")
offsetID=int(sys.argv[2])

ids={}
for line in input:
  if "#" in line:
    continue
  if "END" in line:
    break
  entries=line.split()
  if "TIP3" in line:
    id=int(entries[offsetID])
    if ids.has_key(id):
      ids[id]=ids[id]+1
    else:
      ids[id]=1
    #wrout("for key %d: %d present\n" % (id, ids[id]))

wrout("The following ids will be removed: ")
for key in ids.keys():
  if ids[key] < 3:
    wrout("(%d = %d) " % (key, ids[key]))
wrout("\n")

input.seek(0)
origin=[0., 0., 0.]
id=1
for line in input:
  if "#" in line:
    output.write(line)
    continue
  if "END" in line:
    output.write(line)
    break
  entries=line.split()
  if "TIP3" in line:
    if (ids[int(entries[offsetID])] == 3):
      #wrout("%d present, keeping line: %s" % (ids[id], line))
      output.write("ATOM %6d %-4s%3s%2s %4s   %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s  \n" % (id, entries[2], entries[3], entries[4], entries[5], float(entries[6])-origin[0], float(entries[7])-origin[1], float(entries[8])-origin[2], float(entries[9]), float(entries[10]), entries[11], entries[2][0]))
      id+=1
    else:
     wrout("Only %d present, dropping line: %s" % (ids[int(entries[offsetID])], line))
  else:
    output.write("ATOM %6d %-4s %3s%3s%4s   %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s  \n" % (id, entries[2], entries[3], entries[4], entries[5], float(entries[6])-origin[0], float(entries[7])-origin[1], float(entries[8])-origin[2], float(entries[9]), float(entries[10]), entries[11], entries[2][0]))
    id+=1
input.close()
output.close()