source: util/src/GetBondedPairs.py.in@ 4737be

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

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

# check for parameters
if len(sys.argv) < 5:
  wrerr("Usage: "+sys.argv[0]+" <datafile> <type1> <type2> <distance>\n")
  sys.exit(1)

#Get box boundaries (open parameter file)
input = open(sys.argv[1], "r")
type1 = sys.argv[2]
type2 = sys.argv[3]
distance = float(sys.argv[4])
filetype=sys.argv[1].split(".")
if (filetype[-1]=="data"):
  offset=1
  elementOffset=0
  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()
  offset=1
  elementOffset=0
elif (filetype[-1]=="pdb"):
  offset=5
  elementOffset=2

line=input.readline()
atoms=line.split()
max=[float(atoms[0+offset]), float(atoms[1+offset]), float(atoms[2+offset])]
min=[float(atoms[0+offset]), float(atoms[1+offset]), float(atoms[2+offset])]
for line in input:
  if "END" in line or "CONECT" in line or "MASTER" in line:
    continue
  atom=line.split()
  for i in range(3):
    if float(atom[i+offset]) < min[i]:
      min[i] = float(atom[i+offset])
    if float(atom[i+offset]) > max[i]:
      max[i] = float(atom[i+offset])
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))
cells_y=int(math.ceil(float(max[1])/CELLLENGTH))
cells_z=int(math.ceil(float(max[2])/CELLLENGTH))
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
line=input.readline()
line=input.readline()
for line in input:
  if "END" in line or "CONECT" in line or "MASTER" in line:
    continue
  # atom=[]
  fields=line.split()
  x=int(math.floor((float(fields[0+offset])-min[0])/CELLLENGTH))                # shift coordinates by box minimum
  y=int(math.floor((float(fields[1+offset])-min[1])/CELLLENGTH))
  z=int(math.floor((float(fields[2+offset])-min[2])/CELLLENGTH))
  typeA=fields[elementOffset]
  atom=[n, float(fields[0+offset]), float(fields[1+offset]), float(fields[2+offset]), typeA]
  cell[x][y][z][0] +=1
  cell[x][y][z].append(atom)
  n=n+1


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

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+1):
          #if on border continue periodic
          if e>cells_x-1:
            e=e-cells_x
          for r in range(j-1,j+1):
            #if on border continue periodic
            if r>cells_y-1:
              r=r-cells_y
            for t in range(k-1,k+2):
              #if on boarder continue periodic
              if t>cells_z-1:
                t=t-cells_z
              #efficiency in linked cell: ommit four boxes in middle layer
              if (i==e and r==j+1):
                continue
              if (i==e and r==j+1 and t==k+1):
                continue
              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])
                if (dist < distance*distance):
                  minid=atom1[0]
                  minid2=atom2[0]
                  min=dist
                  if (type1 in atom1[4] and type2 in atom2[4]) or (type2 in atom[4] and type1 in atom2[4]):
                    #print "%s%d\t%s%d\t%f" % (atom1[4],minid, atom2[4],minid2,min)
                    print "%d\t%d" % (minid, minid2)
                
#print v