| 1 | #!@PYTHON@
|
|---|
| 2 |
|
|---|
| 3 | # This programm takes as input an xyz file containing the border atoms of a cluster which is expected to consist of two touching clusters.
|
|---|
| 4 | # It will additionally take a .dat file containing the Convex Hull of the border atoms.
|
|---|
| 5 | # It will remove the Convex Hull and those atoms very close to it
|
|---|
| 6 | #
|
|---|
| 7 | # by Christian Neuen
|
|---|
| 8 |
|
|---|
| 9 |
|
|---|
| 10 | import sys, string, re, math
|
|---|
| 11 | wrerr = sys.stderr.write
|
|---|
| 12 | wrout = sys.stdout.write
|
|---|
| 13 | CUTOFF=20.0
|
|---|
| 14 | DirectNeighbourDistance=17.0
|
|---|
| 15 | eps=10**-2
|
|---|
| 16 |
|
|---|
| 17 |
|
|---|
| 18 | # check for parameters
|
|---|
| 19 | if len(sys.argv) < 3:
|
|---|
| 20 | wrerr("Usage: "+sys.argv[0]+" <datafile>.xyz <datafile>.dat <timestep> \n")
|
|---|
| 21 | sys.exit(1)
|
|---|
| 22 |
|
|---|
| 23 |
|
|---|
| 24 |
|
|---|
| 25 | def distance(a, b):
|
|---|
| 26 | """this functions will find and return the distance between a and b in 3 dim
|
|---|
| 27 | """
|
|---|
| 28 | dist =0.0
|
|---|
| 29 | for i in range(3):
|
|---|
| 30 | dist += (a[i]-b[i])*(a[i]-b[i])
|
|---|
| 31 | return(math.sqrt(dist))
|
|---|
| 32 |
|
|---|
| 33 |
|
|---|
| 34 |
|
|---|
| 35 | input=open(sys.argv[1], "r")
|
|---|
| 36 |
|
|---|
| 37 | line=input.readline() #here is number of atoms
|
|---|
| 38 | line=input.readline() #this is an empty line
|
|---|
| 39 | line=input.readline() #this is the first atom
|
|---|
| 40 | atom=line.split()
|
|---|
| 41 | maxco=[float(atom[1]), float(atom[2]), float(atom[3])]
|
|---|
| 42 | minco=[float(atom[1]), float(atom[2]), float(atom[3])]
|
|---|
| 43 |
|
|---|
| 44 | AllAtoms=[[atom[0], float(atom[1]), float(atom[2]), float(atom[3]), [] ]] # Atom id, Position 1-3, Storage for id and distance to ConvexHull
|
|---|
| 45 | # or Concavity Center)
|
|---|
| 46 |
|
|---|
| 47 |
|
|---|
| 48 | n=0
|
|---|
| 49 |
|
|---|
| 50 | for line in input:
|
|---|
| 51 | n=n+1
|
|---|
| 52 | atom=line.split()
|
|---|
| 53 | for i in range(3):
|
|---|
| 54 | if float(atom[i+1]) < minco[i]:
|
|---|
| 55 | minco[i] = float(atom[i+1])
|
|---|
| 56 | if float(atom[i+1]) > maxco[i]:
|
|---|
| 57 | maxco[i] = float(atom[i+1])
|
|---|
| 58 | AllAtoms.append([atom[0], float(atom[1]), float(atom[2]), float(atom[3]), [] ])
|
|---|
| 59 |
|
|---|
| 60 | BorderAtoms=AllAtoms[:]
|
|---|
| 61 |
|
|---|
| 62 | print len(AllAtoms), "=len of all atoms"
|
|---|
| 63 |
|
|---|
| 64 | print "Found Box bounds [%f,%f]x[%f,%f]x[%f,%f]." % (minco[0],maxco[0],minco[1],maxco[1],minco[2],maxco[2])
|
|---|
| 65 |
|
|---|
| 66 |
|
|---|
| 67 | input=open(sys.argv[2], "r")
|
|---|
| 68 |
|
|---|
| 69 | Temp=[]
|
|---|
| 70 | NumConBorderAtoms=[]
|
|---|
| 71 | NumConBorderFaces=[]
|
|---|
| 72 | ConBorderAtoms=[]
|
|---|
| 73 | ConBorderFaces=[]
|
|---|
| 74 |
|
|---|
| 75 | line = input.readline() #No relevant Info
|
|---|
| 76 | line = input.readline() #No relevant Info
|
|---|
| 77 | line = input.readline()
|
|---|
| 78 | Temp=line.split("=") #We get the # of Atoms in Convex Border and the # of Faces in Convex Hull
|
|---|
| 79 | NumConBorderAtoms=Temp[2].split(",")[0]
|
|---|
| 80 | NumConBorderFaces=Temp[3].split(",")[0]
|
|---|
| 81 |
|
|---|
| 82 | print NumConBorderAtoms, "Atoms mark the Convex Hull"
|
|---|
| 83 | print NumConBorderFaces, "Faces mark the Convex Hull"
|
|---|
| 84 |
|
|---|
| 85 | for i in range(int(NumConBorderAtoms)):
|
|---|
| 86 | line=input.readline()
|
|---|
| 87 | Temp=line.split()
|
|---|
| 88 | ConBorderAtoms.append([float(Temp[0])+67.01283, float(Temp[1])+19.49606, float(Temp[2])+36.74105]) ####This should be done differently!!!!!!!!!
|
|---|
| 89 |
|
|---|
| 90 | line =input.readline() #This is an empty line
|
|---|
| 91 |
|
|---|
| 92 | for i in range(int(NumConBorderFaces)):
|
|---|
| 93 | line=input.readline()
|
|---|
| 94 | Temp=line.split()
|
|---|
| 95 | ConBorderFaces.append([int(Temp[0])-1, int(Temp[1])-1, int(Temp[2])-1])
|
|---|
| 96 |
|
|---|
| 97 | #print len(ConBorderAtoms)
|
|---|
| 98 | #print len(ConBorderFaces)
|
|---|
| 99 |
|
|---|
| 100 | MinDist=max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])
|
|---|
| 101 | MinIndex=-1
|
|---|
| 102 | PlaneVector=[0.0, 0.0, 0.0]
|
|---|
| 103 |
|
|---|
| 104 | output=open("test2.xyz", "w") ##############################################Debug
|
|---|
| 105 | output.write("%d \n \n" %(len(AllAtoms)+13))
|
|---|
| 106 |
|
|---|
| 107 | BorderFace=[]
|
|---|
| 108 | SumDist=[]
|
|---|
| 109 |
|
|---|
| 110 | for j in range(int(NumConBorderFaces)): # this calculates an orthogonal vector to the border face. It will be used to get the
|
|---|
| 111 | # distance of an atom to the border
|
|---|
| 112 | betrag=0
|
|---|
| 113 | for l in range(3):
|
|---|
| 114 | PlaneVector[l]=(ConBorderAtoms[ConBorderFaces[j][1]][l-2]-ConBorderAtoms[ConBorderFaces[j][0]][l-2])*(ConBorderAtoms[ConBorderFaces[j][2]][l-1]-ConBorderAtoms[ConBorderFaces[j][0]][l-1]) - (ConBorderAtoms[ConBorderFaces[j][1]][l-1]-ConBorderAtoms[ConBorderFaces[j][0]][l-1])*(ConBorderAtoms[ConBorderFaces[j][2]][l-2]-ConBorderAtoms[ConBorderFaces[j][0]][l-2])
|
|---|
| 115 | betrag+=PlaneVector[l]*PlaneVector[l]
|
|---|
| 116 |
|
|---|
| 117 | betrag=math.sqrt(betrag) #the vector will be normalized.
|
|---|
| 118 | for l in range(3):
|
|---|
| 119 | PlaneVector[l]/=betrag
|
|---|
| 120 | BorderFace.append([ConBorderFaces[j][0], PlaneVector[:]])
|
|---|
| 121 | SumDist.append(0.0)
|
|---|
| 122 |
|
|---|
| 123 | """
|
|---|
| 124 | for j in range(len(BorderFace)):
|
|---|
| 125 | betrag=0
|
|---|
| 126 | for l in range(3):
|
|---|
| 127 | betrag+=BorderFace[j][1][l]*(ConBorderAtoms[ConBorderFaces[j][1]][l]-ConBorderAtoms[BorderFace[j][0]][l])
|
|---|
| 128 | # betrag=math.sqrt(betrag)
|
|---|
| 129 | print betrag
|
|---|
| 130 | """
|
|---|
| 131 |
|
|---|
| 132 | i=-1
|
|---|
| 133 | skalarp=0
|
|---|
| 134 | MaxDist=0
|
|---|
| 135 | MinDist=0
|
|---|
| 136 |
|
|---|
| 137 |
|
|---|
| 138 | while 1: # we go through all atoms and throw those away which are very close to the convex border
|
|---|
| 139 | # Min Dist is the minimum Distance of a single atom to all enclosing faces
|
|---|
| 140 | MinDist=max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])
|
|---|
| 141 | i+=1
|
|---|
| 142 | if i >= len(AllAtoms):
|
|---|
| 143 | break
|
|---|
| 144 | for j in range(len(BorderFace)):
|
|---|
| 145 | skalarp=0
|
|---|
| 146 | for l in range(3):
|
|---|
| 147 | skalarp+=(BorderFace[j][1][l])*(AllAtoms[i][l+1]-ConBorderAtoms[BorderFace[j][0]][l])
|
|---|
| 148 | if abs(skalarp)<MinDist:
|
|---|
| 149 | MinDist=abs(skalarp)
|
|---|
| 150 | AllAtoms[i][4]=[j, MinDist]
|
|---|
| 151 | if MinDist<3:
|
|---|
| 152 | # print "removing"
|
|---|
| 153 | output.write("H %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
|
|---|
| 154 | AllAtoms.remove(AllAtoms[i])
|
|---|
| 155 | i=i-1
|
|---|
| 156 | break
|
|---|
| 157 |
|
|---|
| 158 |
|
|---|
| 159 | MaxDist1 =-1.0
|
|---|
| 160 | MaxDist2 =-1.0
|
|---|
| 161 | MaxDist3 =-1.0
|
|---|
| 162 |
|
|---|
| 163 | for j in range(len(BorderFace)):
|
|---|
| 164 | MaxDist1 =-1.0
|
|---|
| 165 | MaxDist2 =-1.0
|
|---|
| 166 | MaxDist3 =-1.0
|
|---|
| 167 | for i in range(len(AllAtoms)):
|
|---|
| 168 | skalarp=0
|
|---|
| 169 | for l in range(3):
|
|---|
| 170 | skalarp+=(BorderFace[j][1][l])*(AllAtoms[i][l+1]-ConBorderAtoms[BorderFace[j][0]][l])
|
|---|
| 171 | SumDist[j]+=abs(skalarp)
|
|---|
| 172 |
|
|---|
| 173 | if AllAtoms[i][4][0]==j:
|
|---|
| 174 | if AllAtoms[i][4][1]>MaxDist1:
|
|---|
| 175 | MaxDist3=MaxDist2
|
|---|
| 176 | MaxDist2=MaxDist1
|
|---|
| 177 | MaxDist1=AllAtoms[i][4][1]
|
|---|
| 178 | elif AllAtoms[i][4][1]>MaxDist2:
|
|---|
| 179 | MaxDist3=MaxDist2
|
|---|
| 180 | MaxDist2=AllAtoms[i][4][1]
|
|---|
| 181 | elif AllAtoms[i][4][1]>MaxDist3:
|
|---|
| 182 | MaxDist3=AllAtoms[i][4][1]
|
|---|
| 183 | i=-1
|
|---|
| 184 | while 1:
|
|---|
| 185 | i+=1
|
|---|
| 186 | if i>=len(AllAtoms):
|
|---|
| 187 | break
|
|---|
| 188 | if AllAtoms[i][4][0]==j:
|
|---|
| 189 | if AllAtoms[i][4][1]<MaxDist3:
|
|---|
| 190 | output.write("H %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
|
|---|
| 191 | AllAtoms.remove(AllAtoms[i])
|
|---|
| 192 | i=i-1
|
|---|
| 193 | #print MaxDist1
|
|---|
| 194 | #print MaxDist2
|
|---|
| 195 | #print MaxDist3
|
|---|
| 196 |
|
|---|
| 197 |
|
|---|
| 198 |
|
|---|
| 199 |
|
|---|
| 200 |
|
|---|
| 201 | ########################################################
|
|---|
| 202 | #
|
|---|
| 203 | # Find direct and indirect neighbours
|
|---|
| 204 | # Remove sucht atoms which are isolated as a Concavity
|
|---|
| 205 | #
|
|---|
| 206 | #
|
|---|
| 207 | #
|
|---|
| 208 | ########################################################
|
|---|
| 209 |
|
|---|
| 210 | Neighbours=[]
|
|---|
| 211 | for i in range(len(AllAtoms)):
|
|---|
| 212 | Neighbours.append([])
|
|---|
| 213 | for i in range(len(AllAtoms)):
|
|---|
| 214 | for j in range(len(AllAtoms)):
|
|---|
| 215 | Neighbours[i].append(-1) # Initialized neighbour matrix. entry -1: not connected, entry 1 directly connected, entry 2 indirect connected
|
|---|
| 216 |
|
|---|
| 217 | for i in range(len(AllAtoms)):
|
|---|
| 218 | for j in range(len(AllAtoms)): # First we only look for direkt neighbours
|
|---|
| 219 | if distance(AllAtoms[i][1:4], AllAtoms[j][1:4])<DirectNeighbourDistance:
|
|---|
| 220 | Neighbours[i][j]=1
|
|---|
| 221 | Neighbours[j][i]=1
|
|---|
| 222 |
|
|---|
| 223 | for i in range(len(AllAtoms)):
|
|---|
| 224 | for j in range(len(AllAtoms)):
|
|---|
| 225 | for k in range(len(AllAtoms)):
|
|---|
| 226 | if Neighbours[j][k]<0:
|
|---|
| 227 | for l in range(len(AllAtoms)):
|
|---|
| 228 | if Neighbours[j][l]>0 and Neighbours[k][l]>0: #Now we mark indirekt Neigbours (Is N^4, as it is AllPairsShortestPath + Keeping the Information.
|
|---|
| 229 | Neighbours[j][k]=2
|
|---|
| 230 |
|
|---|
| 231 |
|
|---|
| 232 |
|
|---|
| 233 |
|
|---|
| 234 | i=-1 # Here we remove isolated concave atoms, we expect that a Concavity has more than one Atom.
|
|---|
| 235 | while 1:
|
|---|
| 236 | i+=1
|
|---|
| 237 | if i >=len(AllAtoms):
|
|---|
| 238 | break
|
|---|
| 239 | MaxDist=-0.5
|
|---|
| 240 | for j in range(len(AllAtoms)):
|
|---|
| 241 | if Neighbours[i][j]>0:
|
|---|
| 242 | MaxDist+=1
|
|---|
| 243 | if MaxDist<2:
|
|---|
| 244 | output.write("H %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
|
|---|
| 245 | for j in range(len(AllAtoms)):
|
|---|
| 246 | Neighbours[j][i]=5
|
|---|
| 247 | Neighbours[j].remove(Neighbours[j][i])
|
|---|
| 248 | Neighbours[i]=5
|
|---|
| 249 | Neighbours.remove(Neighbours[i])
|
|---|
| 250 | AllAtoms.remove(AllAtoms[i])
|
|---|
| 251 | i-=1
|
|---|
| 252 |
|
|---|
| 253 |
|
|---|
| 254 |
|
|---|
| 255 |
|
|---|
| 256 |
|
|---|
| 257 |
|
|---|
| 258 | #######################################################
|
|---|
| 259 | #
|
|---|
| 260 | # We have found all those indirectly connected.
|
|---|
| 261 | # We now sort those into Concavities
|
|---|
| 262 | # We calculate Centers of Concavities
|
|---|
| 263 | # We find distance to centers of concavities
|
|---|
| 264 | # Find 10th furthest distance.
|
|---|
| 265 | #
|
|---|
| 266 | ######################################################
|
|---|
| 267 |
|
|---|
| 268 |
|
|---|
| 269 | ConcaveClusters=[]
|
|---|
| 270 |
|
|---|
| 271 | for i in range(len(AllAtoms)): # Here we create a list of Concavities and the atoms of which it consists
|
|---|
| 272 | if Neighbours[i][i]>0:
|
|---|
| 273 | ConcaveClusters.append([i])
|
|---|
| 274 | for j in range(i+1, len(AllAtoms)):
|
|---|
| 275 | if Neighbours[i][j]>0:
|
|---|
| 276 | ConcaveClusters[-1].append(j)
|
|---|
| 277 | Neighbours[j][j]=-1
|
|---|
| 278 |
|
|---|
| 279 |
|
|---|
| 280 | print "Number of all Atoms", len(AllAtoms)
|
|---|
| 281 | print "ConcaveClusters",ConcaveClusters
|
|---|
| 282 | print "Number of ConcaveClusters", len(ConcaveClusters)
|
|---|
| 283 |
|
|---|
| 284 |
|
|---|
| 285 |
|
|---|
| 286 | ConcaveCenters=[] #this will be the center of the concave atoms
|
|---|
| 287 | MaxDist=0
|
|---|
| 288 | primal=[0.0, 0.0, 0.0]
|
|---|
| 289 |
|
|---|
| 290 | for i in range(len(ConcaveClusters)):
|
|---|
| 291 | ConcaveCenters.append([0.0, 0.0, 0.0])
|
|---|
| 292 | # j=0
|
|---|
| 293 | for q in ConcaveClusters[i]:
|
|---|
| 294 | for k in range(3):
|
|---|
| 295 | ConcaveCenters[i][k]+=AllAtoms[q][k+1]
|
|---|
| 296 | # ConcaveCenters[i][k]=(ConcaveCenters[i][k]*float(j)+AllAtoms[q][k+1]) / float(j+1)
|
|---|
| 297 | # j+=1
|
|---|
| 298 | for k in range(3):
|
|---|
| 299 | ConcaveCenters[i][k]/=float(len(ConcaveClusters[i]))
|
|---|
| 300 |
|
|---|
| 301 | print "ConcaveCenters", ConcaveCenters
|
|---|
| 302 |
|
|---|
| 303 |
|
|---|
| 304 | MaxDist=0
|
|---|
| 305 | j=0
|
|---|
| 306 | for i in range(len(AllAtoms)):
|
|---|
| 307 | if AllAtoms[i][4][1]>20:
|
|---|
| 308 | for q in range(3):
|
|---|
| 309 | primal[q]=(primal[q]*j+AllAtoms[i][q+1])/(j+1)
|
|---|
| 310 | j+=1
|
|---|
| 311 |
|
|---|
| 312 |
|
|---|
| 313 | output.write("Cs %f %f %f \n" %(primal[0], primal[1], primal[2]))
|
|---|
| 314 |
|
|---|
| 315 |
|
|---|
| 316 |
|
|---|
| 317 | for i in range(len(ConcaveClusters)):
|
|---|
| 318 | MinDist=[max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])]
|
|---|
| 319 | for j in ConcaveClusters[i]:
|
|---|
| 320 | dist=distance(AllAtoms[j][1:4], ConcaveCenters[i])
|
|---|
| 321 |
|
|---|
| 322 | AllAtoms[j][4]=[i, dist] #We overwrite the distance to the convex hull with distance to Concavity center. Same for respective id.
|
|---|
| 323 | for l in range(len(MinDist)): #We keep track of the 10 closest Atoms to Concavity Center.
|
|---|
| 324 | if dist<MinDist[l]:
|
|---|
| 325 | MinDist.insert(l, dist)
|
|---|
| 326 | if len(MinDist)>10:
|
|---|
| 327 | del MinDist[-1]
|
|---|
| 328 | break
|
|---|
| 329 | ConcaveClusters[i].append(MinDist[-1]) #respective 10th closest distance for center is saved here.
|
|---|
| 330 |
|
|---|
| 331 |
|
|---|
| 332 | ##############################################
|
|---|
| 333 | #
|
|---|
| 334 | # Debug
|
|---|
| 335 | #
|
|---|
| 336 | #
|
|---|
| 337 | #
|
|---|
| 338 | ##############################################
|
|---|
| 339 |
|
|---|
| 340 |
|
|---|
| 341 |
|
|---|
| 342 | for i in range(len(ConcaveClusters)):
|
|---|
| 343 | for j in range(len(ConcaveClusters[i])-1):
|
|---|
| 344 | if i ==0:
|
|---|
| 345 | output.write("O %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
|
|---|
| 346 | if i ==1:
|
|---|
| 347 | output.write("Li %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
|
|---|
| 348 | if i ==2:
|
|---|
| 349 | output.write("Fe %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
|
|---|
| 350 | if i ==3:
|
|---|
| 351 | output.write("He %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
|
|---|
| 352 |
|
|---|
| 353 |
|
|---|
| 354 | output.write("Ne %f %f %f \n" %(ConcaveCenters[i][0], ConcaveCenters[i][1], ConcaveCenters[i][2]))
|
|---|
| 355 |
|
|---|
| 356 |
|
|---|
| 357 | """
|
|---|
| 358 |
|
|---|
| 359 | #output.write("Si %f %f %f \n" %(primal[0], primal[1], primal[2]))
|
|---|
| 360 | #output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][0]][0], ConBorderAtoms[ConBorderFaces[Furthest][0]][1], ConBorderAtoms[ConBorderFaces[Furthest][0]][2]))
|
|---|
| 361 | #output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][1]][0], ConBorderAtoms[ConBorderFaces[Furthest][1]][1], ConBorderAtoms[ConBorderFaces[Furthest][1]][2]))
|
|---|
| 362 | #output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][2]][0], ConBorderAtoms[ConBorderFaces[Furthest][2]][1], ConBorderAtoms[ConBorderFaces[Furthest][2]][2]))
|
|---|
| 363 |
|
|---|
| 364 | #output.close()
|
|---|
| 365 | """
|
|---|
| 366 |
|
|---|
| 367 |
|
|---|
| 368 | ####################################################
|
|---|
| 369 | #
|
|---|
| 370 | # Until here we removed convex atoms. All remaining atoms
|
|---|
| 371 | # belong to some kind of concavity
|
|---|
| 372 | # We will now put them into the kartesian planes in polar coordinates and assign a degree of concavity:
|
|---|
| 373 | # max(len1*cos(alpha), len2*cos(beta)) / lenMid
|
|---|
| 374 | # Sort bei degree of concavity.
|
|---|
| 375 | # Test over best 1*10*20
|
|---|
| 376 | #
|
|---|
| 377 | ####################################################
|
|---|
| 378 | """
|
|---|
| 379 |
|
|---|
| 380 | def MergeSort(List, Typus):
|
|---|
| 381 | """# applies merge sort to list given
|
|---|
| 382 | """
|
|---|
| 383 | if len(List)<2:
|
|---|
| 384 | return List
|
|---|
| 385 | elif Typus==1:
|
|---|
| 386 | List1=List[:int(len(List)/2)]
|
|---|
| 387 | List2=List[int(len(List)/2):]
|
|---|
| 388 | return(MergeAngle(MergeSort(List1, 1), MergeSort(List2, 1)))
|
|---|
| 389 | elif Typus==2:
|
|---|
| 390 | List1=List[:int(len(List)/2)]
|
|---|
| 391 | List2=List[int(len(List)/2):]
|
|---|
| 392 | return(MergeKrit(MergeSort(List1, 2), MergeSort(List2, 2)))
|
|---|
| 393 |
|
|---|
| 394 |
|
|---|
| 395 | def MergeAngle(List1, List2):
|
|---|
| 396 | """# merges two lists in ordered way
|
|---|
| 397 | """
|
|---|
| 398 | List=[]
|
|---|
| 399 | while len(List1)>0 and len(List2)>0:
|
|---|
| 400 | if List1[0][4][0][1]<List2[0][4][0][1]:
|
|---|
| 401 | List.append(List1[0])
|
|---|
| 402 | List1.remove(List1[0])
|
|---|
| 403 | else:
|
|---|
| 404 | List.append(List2[0])
|
|---|
| 405 | List2.remove(List2[0])
|
|---|
| 406 | while len(List1)>0:
|
|---|
| 407 | List.append(List1[0])
|
|---|
| 408 | List1.remove(List1[0])
|
|---|
| 409 | while len(List2)>0:
|
|---|
| 410 | List.append(List2[0])
|
|---|
| 411 | List2.remove(List2[0])
|
|---|
| 412 | return(List)
|
|---|
| 413 |
|
|---|
| 414 |
|
|---|
| 415 |
|
|---|
| 416 |
|
|---|
| 417 | def MergeKrit(List1, List2):
|
|---|
| 418 | """# merges two lists in ordered way
|
|---|
| 419 | """
|
|---|
| 420 | List=[]
|
|---|
| 421 | while len(List1)>0 and len(List2)>0:
|
|---|
| 422 | if List1[0][4][1]>List2[0][4][1]:
|
|---|
| 423 | List.append(List1[0])
|
|---|
| 424 | List1.remove(List1[0])
|
|---|
| 425 | else:
|
|---|
| 426 | List.append(List2[0])
|
|---|
| 427 | List2.remove(List2[0])
|
|---|
| 428 | while len(List1)>0:
|
|---|
| 429 | List.append(List1[0])
|
|---|
| 430 | List1.remove(List1[0])
|
|---|
| 431 | while len(List2)>0:
|
|---|
| 432 | List.append(List2[0])
|
|---|
| 433 | List2.remove(List2[0])
|
|---|
| 434 | return(List)
|
|---|
| 435 |
|
|---|
| 436 |
|
|---|
| 437 | r=0.0
|
|---|
| 438 | phi=0.0
|
|---|
| 439 | sign=0.0
|
|---|
| 440 | Krit=0.0
|
|---|
| 441 | TempKrit=0.0
|
|---|
| 442 |
|
|---|
| 443 | for i in range(2):
|
|---|
| 444 | for j in range(i+1,3):
|
|---|
| 445 | for k in range(len(AllAtoms)):
|
|---|
| 446 | # print AllAtoms[k][i+1], primal[i]
|
|---|
| 447 | # print AllAtoms[k][j+1], primal[j]
|
|---|
| 448 | r=math.sqrt((AllAtoms[k][i+1]-primal[i])*(AllAtoms[k][i+1]-primal[i])+(AllAtoms[k][j+1]-primal[j])*(AllAtoms[k][j+1]-primal[j])) #r=sqrt(x^2+y^2)
|
|---|
| 449 | if AllAtoms[k][j+1]-primal[j]>0:
|
|---|
| 450 | phi=1
|
|---|
| 451 | else:
|
|---|
| 452 | phi=-1
|
|---|
| 453 | phi*=math.acos((AllAtoms[k][i+1]-primal[i])/r) # phi = sign(y)*arccos(x/r)
|
|---|
| 454 | AllAtoms[k][4]=[[r, phi], Krit]
|
|---|
| 455 |
|
|---|
| 456 | AllAtoms=MergeSort(AllAtoms, 1)
|
|---|
| 457 |
|
|---|
| 458 | for k in range(len(AllAtoms)):
|
|---|
| 459 | # print AllAtoms[k][4][0][1]
|
|---|
| 460 | Krit=((AllAtoms[k-2][4][0][0]*math.cos(AllAtoms[k-1][4][0][1]-AllAtoms[k-2][4][0][1])+ AllAtoms[k][4][0][0]*math.cos(AllAtoms[k][4][0][1]-AllAtoms[k-1][4][0][1]))/2)/AllAtoms[k-1][4][0][0]
|
|---|
| 461 | # Assign the ratio of the projected length of two outer atom legs to the inner leg. Krit > 1 is concavity
|
|---|
| 462 | if Krit>AllAtoms[k-1][4][1]:
|
|---|
| 463 | AllAtoms[k-1][4][1]=Krit
|
|---|
| 464 |
|
|---|
| 465 | AllAtoms=MergeSort(AllAtoms, 2)
|
|---|
| 466 |
|
|---|
| 467 | #for i in range(len(AllAtoms)):
|
|---|
| 468 | # print AllAtoms[i][4][1]
|
|---|
| 469 |
|
|---|
| 470 |
|
|---|
| 471 |
|
|---|
| 472 |
|
|---|
| 473 |
|
|---|
| 474 | """
|
|---|
| 475 |
|
|---|
| 476 |
|
|---|
| 477 | ####################################################
|
|---|
| 478 | #
|
|---|
| 479 | #
|
|---|
| 480 | #
|
|---|
| 481 | # We will now iterate over all triples and check the bonds
|
|---|
| 482 | # being cut by the respective plane
|
|---|
| 483 | #
|
|---|
| 484 | ####################################################
|
|---|
| 485 |
|
|---|
| 486 |
|
|---|
| 487 | CompleteAtoms =[]
|
|---|
| 488 | CompleteBonds =[]
|
|---|
| 489 |
|
|---|
| 490 | input = open("silica.grape.%.4d" %(int(sys.argv[3])), "r")
|
|---|
| 491 | line=input.readline()
|
|---|
| 492 |
|
|---|
| 493 | for line in input:
|
|---|
| 494 | CompleteAtoms.append([float(line.split()[1]), float(line.split()[2]), float(line.split()[3]) ])
|
|---|
| 495 | input.close()
|
|---|
| 496 |
|
|---|
| 497 | input=open("silica.dbond-SiOCaO.%.4d" %(int(sys.argv[3])), "r")
|
|---|
| 498 | line=input.readline()
|
|---|
| 499 |
|
|---|
| 500 | for line in input:
|
|---|
| 501 | CompleteBonds.append([int(line.split()[0]), int(line.split()[1])])
|
|---|
| 502 | input.close()
|
|---|
| 503 |
|
|---|
| 504 |
|
|---|
| 505 | CenterOfTriangle=[0.0, 0.0, 0.0]
|
|---|
| 506 | MinimumPlane=[0, 0, 0]
|
|---|
| 507 | SecondPlane=[0,0,0]
|
|---|
| 508 | ThirdPlane=[0,0,0]
|
|---|
| 509 | MinimumCutCount=len(CompleteBonds)
|
|---|
| 510 | SecondCutCount=len(CompleteBonds)
|
|---|
| 511 | ThirdCutCount=len(CompleteBonds)
|
|---|
| 512 | CutCount=0
|
|---|
| 513 | PlanePreCheck=0
|
|---|
| 514 |
|
|---|
| 515 | PlaneVector=[0.0, 0.0, 0.0]
|
|---|
| 516 | Atom1Vector=[0.0, 0.0, 0.0]
|
|---|
| 517 | Atom2Vector=[0.0, 0.0, 0.0]
|
|---|
| 518 | DirectionVector=[0.0, 0.0, 0.0]
|
|---|
| 519 |
|
|---|
| 520 |
|
|---|
| 521 | iteration=0
|
|---|
| 522 |
|
|---|
| 523 |
|
|---|
| 524 | for i in range(len(ConcaveCenters)):
|
|---|
| 525 | for j in range(i, len(ConcaveCenters)):
|
|---|
| 526 | for k in ConcaveClusters[i]+ConcaveClusters[j]:
|
|---|
| 527 | if k!=int(k) or AllAtoms[k][4][1]>ConcaveClusters[AllAtoms[k][4][0]][-1]:
|
|---|
| 528 | # if AllAtoms[k][4][1]>ConcaveClusters[AllAtoms[k][4][0]][-1] or (AllAtoms[i][4][0]==AllAtoms[j][4][0] and AllAtoms[i][4][0]==AllAtoms[k][4][0]):
|
|---|
| 529 |
|
|---|
| 530 | # (Neighbours[i][j]>0 and Neighbours[i][k]>0)
|
|---|
| 531 | # or not(AllAtoms[i][4][1]>20 or AllAtoms[j][4][1]>20 or AllAtoms[k][4][1]>20) or
|
|---|
| 532 | # (AllAtoms[i][4][1]>20 and AllAtoms[j][4][1]>20 and AllAtoms[k][4][1]>20):
|
|---|
| 533 | # or Neighbours[i][k]==1 or Neighbours[j][k]==1
|
|---|
| 534 | continue
|
|---|
| 535 |
|
|---|
| 536 | for q in range(3):
|
|---|
| 537 | CenterOfTriangle[q]= (ConcaveCenters[i][q] + ConcaveCenters[j][q] + AllAtoms[k][q+1])/3.0
|
|---|
| 538 |
|
|---|
| 539 |
|
|---|
| 540 | ### this calculates the crossproduct of the vectors in the plane, yielding an orthogonal vector on the plane
|
|---|
| 541 | PlaneVector[0] = (ConcaveCenters[j][1]-ConcaveCenters[i][1])*(AllAtoms[k][3]-ConcaveCenters[i][2]) - (ConcaveCenters[j][2]-ConcaveCenters[i][2])*(AllAtoms[k][2]-ConcaveCenters[i][1])
|
|---|
| 542 | PlaneVector[1] = (ConcaveCenters[j][2]-ConcaveCenters[i][2])*(AllAtoms[k][1]-ConcaveCenters[i][0]) - (ConcaveCenters[j][0]-ConcaveCenters[i][0])*(AllAtoms[k][3]-ConcaveCenters[i][2])
|
|---|
| 543 | PlaneVector[2] = (ConcaveCenters[j][0]-ConcaveCenters[i][0])*(AllAtoms[k][2]-ConcaveCenters[i][1]) - (ConcaveCenters[j][1]-ConcaveCenters[i][1])*(AllAtoms[k][1]-ConcaveCenters[i][0])
|
|---|
| 544 |
|
|---|
| 545 | ### cutting plane is done, next plane will ensure, that we only cut one leg of the cluster
|
|---|
| 546 |
|
|---|
| 547 | for q in range(3):
|
|---|
| 548 | DirectionVector[q] = CenterOfTriangle[q]-ConcaveCenters[i][q]
|
|---|
| 549 |
|
|---|
| 550 | """
|
|---|
| 551 | skalarpA=0.0
|
|---|
| 552 | skalarpB=0.0
|
|---|
| 553 | skalarpC=0.0
|
|---|
| 554 | for q in range(3):
|
|---|
| 555 | skalarpA += (ConcaveCenters[i][q]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[k][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
|
|---|
| 556 | skalarpB += (ConcaveCenters[i][hust][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[j][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
|
|---|
| 557 | skalarpC += (AllAtoms[j][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[k][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
|
|---|
| 558 | if skalarpA<0 or skalarpB<0 or skalarpC<0:
|
|---|
| 559 | continue #Concavity seperates triangle
|
|---|
| 560 | """
|
|---|
| 561 |
|
|---|
| 562 | PlanePreCheck=0
|
|---|
| 563 | for l in range(len(BorderAtoms)):
|
|---|
| 564 | skalarp=0
|
|---|
| 565 | for g in range(1,4):
|
|---|
| 566 | skalarp+=(BorderAtoms[l][g]-ConcaveCenters[i][g-1])*PlaneVector[g-1]
|
|---|
| 567 | if skalarp<0:
|
|---|
| 568 | PlanePreCheck+=1
|
|---|
| 569 | if PlanePreCheck<(4.0/5.0*len(BorderAtoms)) and PlanePreCheck>(len(BorderAtoms)/5.0):
|
|---|
| 570 | CutCount=0
|
|---|
| 571 | else:
|
|---|
| 572 | continue
|
|---|
| 573 | print "iteration", iteration
|
|---|
| 574 | iteration+=1
|
|---|
| 575 |
|
|---|
| 576 | for l in range(len(CompleteBonds)): ### Go through all the bonds, create vectors from first point in plane to both atoms
|
|---|
| 577 | if CutCount>=ThirdCutCount: ### This gives small speed enhancement as expensive planes will not be fully computed anymore.
|
|---|
| 578 | break
|
|---|
| 579 |
|
|---|
| 580 | Skalar1=0
|
|---|
| 581 | Skalar2=0
|
|---|
| 582 | # We again work with the sign of the skalar product
|
|---|
| 583 | """
|
|---|
| 584 | for q in range(3): #Calculate Vectors from Concavity to Atoms
|
|---|
| 585 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[i][q]
|
|---|
| 586 | Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[i][q]
|
|---|
| 587 |
|
|---|
| 588 | for q in range(3):
|
|---|
| 589 | Skalar1+=Atom1Vector[q]*DirectionVector[q]
|
|---|
| 590 | if Skalar1<0:
|
|---|
| 591 | continue
|
|---|
| 592 | Skalar1=0
|
|---|
| 593 | """
|
|---|
| 594 | for q in range(3): #Calculate Vectors from plane to Atoms
|
|---|
| 595 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[i][q]
|
|---|
| 596 | Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[i][q]
|
|---|
| 597 |
|
|---|
| 598 | for q in range(3):
|
|---|
| 599 | Skalar1+=Atom1Vector[q]*PlaneVector[q]
|
|---|
| 600 | Skalar2+=Atom2Vector[q]*PlaneVector[q]
|
|---|
| 601 | if Skalar1*Skalar2<0:
|
|---|
| 602 | # print Atom1Vector
|
|---|
| 603 | # print Atom2Vector
|
|---|
| 604 | CutCount+=1
|
|---|
| 605 |
|
|---|
| 606 |
|
|---|
| 607 | if CutCount<ThirdCutCount and CutCount>50:
|
|---|
| 608 | if CutCount<SecondCutCount:
|
|---|
| 609 | if CutCount<MinimumCutCount:
|
|---|
| 610 | ThirdPlane=SecondPlane[:]
|
|---|
| 611 | ThirdCutCount=SecondCutCount
|
|---|
| 612 | SecondPlane=MinimumPlane[:]
|
|---|
| 613 | SecondCutCount=MinimumCutCount
|
|---|
| 614 | MinimumCutCount=CutCount
|
|---|
| 615 | MinimumPlane=[i, j, k]
|
|---|
| 616 | print PlanePreCheck
|
|---|
| 617 | print MinimumCutCount
|
|---|
| 618 | print MinimumPlane
|
|---|
| 619 | else:
|
|---|
| 620 | ThirdPlane=SecondPlane[:]
|
|---|
| 621 | ThirdCutCount=SecondCutCount
|
|---|
| 622 | SecondPlane=[i, j, k]
|
|---|
| 623 | SecondCutCount=CutCount
|
|---|
| 624 | else:
|
|---|
| 625 | ThirdPlane=[i, j, k]
|
|---|
| 626 | ThirdCutCount=CutCount
|
|---|
| 627 |
|
|---|
| 628 |
|
|---|
| 629 |
|
|---|
| 630 | print MinimumCutCount
|
|---|
| 631 | print MinimumPlane
|
|---|
| 632 | print AllAtoms[MinimumPlane[0]]
|
|---|
| 633 | print AllAtoms[MinimumPlane[1]]
|
|---|
| 634 | print AllAtoms[MinimumPlane[2]]
|
|---|
| 635 |
|
|---|
| 636 |
|
|---|
| 637 | print SecondCutCount
|
|---|
| 638 | print SecondPlane
|
|---|
| 639 | print AllAtoms[SecondPlane[0]]
|
|---|
| 640 | print AllAtoms[SecondPlane[1]]
|
|---|
| 641 | print AllAtoms[SecondPlane[2]]
|
|---|
| 642 | print ThirdCutCount
|
|---|
| 643 | print ThirdPlane
|
|---|
| 644 | print AllAtoms[ThirdPlane[0]]
|
|---|
| 645 | print AllAtoms[ThirdPlane[1]]
|
|---|
| 646 | print AllAtoms[ThirdPlane[2]]
|
|---|
| 647 |
|
|---|
| 648 |
|
|---|
| 649 | #for i in range(len(Center)):
|
|---|
| 650 | # output.write("Ca %f %f %f \n" %(Center[i][0], Center[i][1], Center[i][2]))
|
|---|
| 651 |
|
|---|
| 652 | for i in range(2):
|
|---|
| 653 | output.write("Si %f %f %f \n" %(ConcaveCenters[MinimumPlane[i]][0], ConcaveCenters[MinimumPlane[i]][1], ConcaveCenters[MinimumPlane[i]][2]))
|
|---|
| 654 | output.write("Ar %f %f %f \n" %(ConcaveCenters[SecondPlane[i]][0], ConcaveCenters[SecondPlane[i]][1], ConcaveCenters[SecondPlane[i]][2]))
|
|---|
| 655 | output.write("Ca %f %f %f \n" %(ConcaveCenters[ThirdPlane[i]][0], ConcaveCenters[ThirdPlane[i]][1], ConcaveCenters[ThirdPlane[i]][2]))
|
|---|
| 656 |
|
|---|
| 657 | output.write("Si %f %f %f \n" %(AllAtoms[MinimumPlane[2]][1], AllAtoms[MinimumPlane[2]][2], AllAtoms[MinimumPlane[2]][3]))
|
|---|
| 658 | output.write("Ar %f %f %f \n" %(AllAtoms[SecondPlane[2]][1], AllAtoms[SecondPlane[2]][2], AllAtoms[SecondPlane[2]][3]))
|
|---|
| 659 | output.write("Ca %f %f %f \n" %(AllAtoms[ThirdPlane[2]][1], AllAtoms[ThirdPlane[2]][2], AllAtoms[ThirdPlane[2]][3]))
|
|---|
| 660 |
|
|---|
| 661 | output.close()
|
|---|
| 662 |
|
|---|
| 663 |
|
|---|
| 664 | ########################################################
|
|---|
| 665 | #
|
|---|
| 666 | # Having found a cutting plane we will now
|
|---|
| 667 | # seperate the Cluster, first cutting all bonds as
|
|---|
| 668 | # they are cut by the plane, the cutting the Atoms
|
|---|
| 669 | # which were the resting points of the plane
|
|---|
| 670 | #
|
|---|
| 671 | ########################################################
|
|---|
| 672 |
|
|---|
| 673 | output=open("silica.dbond-SiOCaO.Cut", "w")
|
|---|
| 674 |
|
|---|
| 675 | output.write("0 666666 \n")
|
|---|
| 676 |
|
|---|
| 677 |
|
|---|
| 678 | ### this calculates the crossproduct of the vectors in the plane, yielding an orthogonal vector on the plane
|
|---|
| 679 | PlaneVector[0] = (ConcaveCenters[MinimumPlane[1]][1]-ConcaveCenters[MinimumPlane[0]][1])*(AllAtoms[MinimumPlane[2]][3]-ConcaveCenters[MinimumPlane[0]][2]) - (ConcaveCenters[MinimumPlane[1]][2]-ConcaveCenters[MinimumPlane[0]][2])*(AllAtoms[MinimumPlane[2]][2]-ConcaveCenters[MinimumPlane[0]][1])
|
|---|
| 680 |
|
|---|
| 681 | PlaneVector[1] = (ConcaveCenters[MinimumPlane[1]][2]-ConcaveCenters[MinimumPlane[0]][2])*(AllAtoms[MinimumPlane[2]][1]-ConcaveCenters[MinimumPlane[0]][0]) - (ConcaveCenters[MinimumPlane[1]][0]-ConcaveCenters[MinimumPlane[0]][0])*(AllAtoms[MinimumPlane[2]][3]-ConcaveCenters[MinimumPlane[0]][2])
|
|---|
| 682 |
|
|---|
| 683 | PlaneVector[2] = (ConcaveCenters[MinimumPlane[1]][0]-ConcaveCenters[MinimumPlane[0]][0])*(AllAtoms[MinimumPlane[2]][2]-ConcaveCenters[MinimumPlane[0]][1]) - (ConcaveCenters[MinimumPlane[1]][1]-ConcaveCenters[MinimumPlane[0]][1])*(AllAtoms[MinimumPlane[2]][1]-ConcaveCenters[MinimumPlane[0]][0])
|
|---|
| 684 |
|
|---|
| 685 | ### cutting plane is done, next plane will ensure, that we only cut one leg of the cluster
|
|---|
| 686 | absolut=0
|
|---|
| 687 | for q in range(3):
|
|---|
| 688 | absolut+=PlaneVector[q]*PlaneVector[q]
|
|---|
| 689 | absolut=math.sqrt(absolut)
|
|---|
| 690 |
|
|---|
| 691 | for q in range(3):
|
|---|
| 692 | PlaneVector[q]/=absolut
|
|---|
| 693 |
|
|---|
| 694 | print PlaneVector
|
|---|
| 695 |
|
|---|
| 696 | for q in range(3):
|
|---|
| 697 | CenterOfTriangle[q]= (ConcaveCenters[MinimumPlane[0]][q] + ConcaveCenters[MinimumPlane[1]][q] + AllAtoms[MinimumPlane[2]][q+1])/3.0
|
|---|
| 698 | DirectionVector[q] = CenterOfTriangle[q]-ConcaveCenters[MinimumPlane[0]][q]
|
|---|
| 699 |
|
|---|
| 700 |
|
|---|
| 701 | l=-1
|
|---|
| 702 | #for l in range(len(CompleteBonds)): ### Go through all the bonds, create vectors from first point in plane to both atoms
|
|---|
| 703 | while 1:
|
|---|
| 704 | l+=1
|
|---|
| 705 | if l>=len(CompleteBonds):
|
|---|
| 706 | break
|
|---|
| 707 | Skalar1=0
|
|---|
| 708 | Skalar2=0
|
|---|
| 709 | # We again work with the sign of the skalar product
|
|---|
| 710 | """
|
|---|
| 711 | for q in range(3): #Calculate Vectors from Concavity to Atoms
|
|---|
| 712 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]
|
|---|
| 713 |
|
|---|
| 714 | for q in range(3):
|
|---|
| 715 | Skalar1+=Atom1Vector[q]*DirectionVector[q]
|
|---|
| 716 |
|
|---|
| 717 | if Skalar1<0:
|
|---|
| 718 | output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))
|
|---|
| 719 | del CompleteBonds[l]
|
|---|
| 720 | l-=1
|
|---|
| 721 | continue
|
|---|
| 722 |
|
|---|
| 723 | Skalar1=0
|
|---|
| 724 | """
|
|---|
| 725 | for q in range(3): #Calculate Vectors from plane to Atoms
|
|---|
| 726 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]
|
|---|
| 727 | Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]
|
|---|
| 728 |
|
|---|
| 729 | # if distance(Atom2Vector[:], [0, 0, 0])<1:
|
|---|
| 730 | # print Atom1Vector
|
|---|
| 731 | for q in range(3):
|
|---|
| 732 | Skalar1+=Atom1Vector[q]*PlaneVector[q]
|
|---|
| 733 | Skalar2+=Atom2Vector[q]*PlaneVector[q]
|
|---|
| 734 | if abs(Skalar1)>eps and abs(Skalar2)>eps:
|
|---|
| 735 | if Skalar1*Skalar2>0:
|
|---|
| 736 | output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))
|
|---|
| 737 | del CompleteBonds[l]
|
|---|
| 738 | l-=1
|
|---|
| 739 | else:
|
|---|
| 740 | del CompleteBonds[l]
|
|---|
| 741 | l-=1
|
|---|
| 742 | else:
|
|---|
| 743 | print Skalar1, Skalar2
|
|---|
| 744 |
|
|---|
| 745 |
|
|---|
| 746 | ######################
|
|---|
| 747 | #
|
|---|
| 748 | # Now we make small parallel shift and take care of the remaining bonds
|
|---|
| 749 | #
|
|---|
| 750 | #
|
|---|
| 751 | ######################
|
|---|
| 752 |
|
|---|
| 753 | ShiftCount=[0, 0]
|
|---|
| 754 | shift=10**-2
|
|---|
| 755 |
|
|---|
| 756 | for i in range(2):
|
|---|
| 757 | for l in range(len(CompleteBonds)):
|
|---|
| 758 | for q in range(3): #Calculate Vectors from plane to Atoms
|
|---|
| 759 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
|
|---|
| 760 | Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
|
|---|
| 761 | Skalar1=0
|
|---|
| 762 | Skalar2=0
|
|---|
| 763 | for q in range(3):
|
|---|
| 764 | Skalar1+=Atom1Vector[q]*PlaneVector[q]
|
|---|
| 765 | Skalar2+=Atom2Vector[q]*PlaneVector[q]
|
|---|
| 766 | if Skalar1*Skalar2<0:
|
|---|
| 767 | ShiftCount[i]+=1
|
|---|
| 768 |
|
|---|
| 769 | if ShiftCount[0]>ShiftCount[1]:
|
|---|
| 770 | i=1
|
|---|
| 771 | else:
|
|---|
| 772 | i=0
|
|---|
| 773 |
|
|---|
| 774 |
|
|---|
| 775 | for l in range(len(CompleteBonds)):
|
|---|
| 776 | for q in range(3): #Calculate Vectors from plane to Atoms
|
|---|
| 777 | Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
|
|---|
| 778 | Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
|
|---|
| 779 | Skalar1=0
|
|---|
| 780 | Skalar2=0
|
|---|
| 781 | for q in range(3):
|
|---|
| 782 | Skalar1+=Atom1Vector[q]*PlaneVector[q]
|
|---|
| 783 | Skalar2+=Atom2Vector[q]*PlaneVector[q]
|
|---|
| 784 | if Skalar1*Skalar2>0:
|
|---|
| 785 | output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))
|
|---|
| 786 |
|
|---|
| 787 |
|
|---|
| 788 |
|
|---|
| 789 | print len(CompleteBonds)
|
|---|
| 790 |
|
|---|
| 791 |
|
|---|
| 792 |
|
|---|
