#!@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]+" " 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"