source: util/src/CheckMap.py.in@ d97af9

#!@PYTHON@
#
# checks whether a given bond file actually makes sense in terms of a max distance

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

# check arguments
if len(sys.argv) < 3:
        print "Usage: "+sys.argv[0]+" <srcBondFile> <srcPositionFile> <max_dist>"
        sys.exit(1)


EPSILON=1e-3
CUTOFF=2.
inputBonds = open(sys.argv[1], "r")
inputPos = open(sys.argv[2], "r")
maxdist = (float)(sys.argv[3])

def distance(atom1, atom2):
  """# Here we calculate the distance between atom1 and atom2
  """
  dist = 0
  for i in range(len(atom1)):
    dist+=(atom1[i]-atom2[i])*(atom1[i]-atom2[i])
  return(  math.sqrt(dist))



print "Scanning XYZ file"+sys.argv[2]+"."
Atoms = []
Bonds = []
x = [0.0, 0.0, 0.0]

for line in inputPos:
        entries = line.split()
        if (len(entries)<1 or (entries[0]!="O" and entries[0]!="Si" and entries[0]!="Ca")):
                continue
        for n in range(3):
                x[n] = float(entries[1+n])
        Atoms.append([x[0], x[1], x[2]])
inputPos.close()

for line in inputBonds:
  entries = line.split()
  Bonds.append([int(entries[0]), int(entries[1])])
inputBonds.close()

NumViolators =0

for i in range(len(Bonds)):
  try:
    dist =distance(Atoms[Bonds[i][0]], Atoms[Bonds[i][1]])
  except IndexError:
    print "Bonds["+str(i)+"][0] = "+str(Bonds[i][0])+" or Bonds["+str(i)+"][1] = "+str(Bonds[i][1])+" does not exist."
    dist=89
        #if Bonds[i][0]==6096:
  #  print "Atoms ", Bonds[i], "have distance: ", dist
  #  print Atoms[Bonds[i][0]]
  #  print Atoms[Bonds[i][1]]
  if dist>maxdist and dist<88:
    NumViolators+=1
    #print "Atoms ", Bonds[i], Atoms[Bonds[i][0]], Atoms[Bonds[i][1]], "have distance: ", dist

print NumViolators, "Atoms have unreasonable bonds"