source: util/src/GetBondsFromPDB.py.in@ b4fa915

#!@PYTHON@
#
# determines all distance between each atomic pair and prints minimum and respective ids

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

# check for parameters
if len(sys.argv) < 6:
  wrerr("Usage: "+sys.argv[0]+" <datafile> <bondfile> <offsetID> <offsetElement> <offsetXYZ> <output>\n")
  sys.exit(1)

#Get box boundaries (open parameter file)
input = open(sys.argv[1], "r")
inputbonds = open(sys.argv[2], "r")
offsetID=int(sys.argv[3])
offsetElement=int(sys.argv[4])
offsetXYZ=int(sys.argv[5])

# read in bond criteria
CELLLENGTH=0.
bonds=[]
for line in inputbonds:
  entries=line.split()
  if float(entries[2]) > CELLLENGTH:
    CELLLENGTH=float(entries[2])
  bonds.append([ entries[0], entries[1], float(entries[2])])
#  bonds.append([ entries[1], entries[0], float(entries[2])])
CELLLENGTH+=THRESHOLD

print str(CELLLENGTH)+" length of cell."

filetype=sys.argv[1].split(".")
if (filetype[-1]=="data"):
  while 1:
    line=input.readline()       # skip header (starting with #)
    if not re.compile("^#").match(line):
      break
elif (filetype[-1]=="xyz"):
  line=input.readline()
  line=input.readline()

line=input.readline()
atoms=line.split()
max=[float(atoms[offsetXYZ]), float(atoms[offsetXYZ+1]), float(atoms[offsetXYZ+2])]
min=[float(atoms[offsetXYZ]), float(atoms[offsetXYZ+1]), float(atoms[offsetXYZ+2])]
for line in input:
  if "#" in line:
    continue
  if "END" in line:
    break
  atom=line.split()
  for i in range(3):
    if float(atom[offsetXYZ+i]) < min[i]:
      min[i] = float(atom[offsetXYZ+i])
    if float(atom[offsetXYZ+i]) > max[i]:
      max[i] = float(atom[offsetXYZ+i])
print "Found Box bounds [%f,%f]x[%f,%f]x[%f,%f]." % (min[0],max[0],min[1],max[1],min[2],max[2])
for i in range(3):      # shift by minimum if below zero
  if min[i] < 0:
    max[i]-=min[i]
  else:
    min[i]=0
cells_x=int(math.ceil(float(max[0])/CELLLENGTH))+1
cells_y=int(math.ceil(float(max[1])/CELLLENGTH))+1
cells_z=int(math.ceil(float(max[2])/CELLLENGTH))+1
print "Number of cells in each axis direction (%f,%f,%f)." % (cells_x, cells_y, cells_z)
input.close()

#create 3 dim array to fit in box, each cell has size 5 Ang
cell=[]

for i in range(cells_x):
  cell.append([])
  for j in range(cells_y):
    cell[i].append([])
    for k in range(cells_z):
      cell[i][j].append([0])


# open files
input = open(sys.argv[1], "r")

# parse every atom into appropriate cell
n=0
filetype=sys.argv[1].split(".")
if (filetype[-1]=="data"):
  while 1:
    line=input.readline()       # skip header (starting with #)
    if not re.compile("^#").match(line):
      break
elif (filetype[-1]=="xyz"):
  line=input.readline()
  line=input.readline()
for line in input:
  if "#" in line:
    continue
  if "END" in line:
    break
  # atom=[]
  fields=line.split()
  x=int(math.floor((float(fields[offsetXYZ+0])-min[0])/CELLLENGTH))             # shift coordinates by box minimum
  y=int(math.floor((float(fields[offsetXYZ+1])-min[1])/CELLLENGTH))
  z=int(math.floor((float(fields[offsetXYZ+2])-min[2])/CELLLENGTH))
  atom=[int(fields[offsetID]), float(fields[offsetXYZ+0]), float(fields[offsetXYZ+1]), float(fields[offsetXYZ+2]), fields[offsetElement]]
  if int(fields[offsetID]) == 0:
    print "error"
  #wrout("Appending atom nr.%d (%f,%f,%f)in (%d,%d,%d)\n" % (n, float(fields[offsetXYZ+1]), float(fields[offsetXYZ+2]), float(fields[offsetXYZ+3]), x,y,z))
  cell[x][y][z][0] +=1
  cell[x][y][z].append(atom)
  n=n+1
print "I parsed "+str(n)+" atoms."

# go through every cell
min = 100000.
minid=-1
minid2=-1
atom1=[0,0.,0.,0.]
atom2=[0,0.,0.,0.]
v=0
nr=0

output  = open(sys.argv[6], "w")
for i in range(cells_x):
  for j in range(cells_y):
    for k in range(cells_z):
      v+=cell[i][j][k][0]

      #go through every atom in cell
      for l in range(1, cell[i][j][k][0]+1):
        atom1=cell[i][j][k][l]
  
        # go through cell and all lower neighbours
        for e in range(i-1,i+2):
          #if on border continue periodic
          if e<0 or e>cells_x-1:
            break
          for r in range(j-1,j+2):
            #if on border continue periodic
            if r<0 or r>cells_y-1:
              break
            for t in range(k-1,k+2):
              #if on boarder continue periodic
              if t<0 or t>cells_z-1:
                break
              w=cell[e][r][t][0]
              #go through all atoms in cell
              for m in range(1, w+1):
                atom2=cell[e][r][t][m]
                dist=0
                #when in same cell: ommit identical particles
                if (atom1[0]==atom2[0]):
                  continue
                for h in range(1,4):
                  dist += (atom1[h] - atom2[h])*(atom1[h] - atom2[h])
                for type in range(len(bonds)):
                  if (atom1[4] == bonds[type][0]) and (atom2[4] == bonds[type][1]):
                    if ((dist < (bonds[type][2]+THRESHOLD)*(bonds[type][2]+THRESHOLD)) and (dist > ((bonds[type][2]-THRESHOLD))*((bonds[type][2]-THRESHOLD)))):
                      output.write("%d\t%d\n" % (atom1[0], atom2[0]))
                      nr+=1

output.close()                
#print v
print str(nr)+" bonds were found."