#!@PYTHON@

# This programm takes as input an xyz file containing the border atoms of a cluster which is expected to consist of two touching clusters.
# It will additionally take a .dat file containing the Convex Hull of the border atoms.
# It will remove the Convex Hull and those atoms very close to it
# 
# by Christian Neuen


import sys, string, re, math
wrerr = sys.stderr.write
wrout = sys.stdout.write
CUTOFF=20.0
DirectNeighbourDistance=17.0
eps=10**-2


# check for parameters
if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <datafile>.xyz <datafile>.dat <timestep> \n")
  sys.exit(1)



def distance(a, b): 
    """this functions will find and return the distance between a and b in 3 dim
    """
    dist =0.0
    for i in range(3):
      dist += (a[i]-b[i])*(a[i]-b[i])
    return(math.sqrt(dist))



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

line=input.readline() #here is number of atoms
line=input.readline() #this is an empty line
line=input.readline() #this is the first atom
atom=line.split()
maxco=[float(atom[1]), float(atom[2]), float(atom[3])]
minco=[float(atom[1]), float(atom[2]), float(atom[3])]

AllAtoms=[[atom[0], float(atom[1]), float(atom[2]), float(atom[3]), [] ]] # Atom id, Position 1-3, Storage for id and distance to ConvexHull 
                                                                            # or Concavity Center)


n=0

for line in input:
  n=n+1  
  atom=line.split()
  for i in range(3):
    if float(atom[i+1]) < minco[i]:
      minco[i] = float(atom[i+1])
    if float(atom[i+1]) > maxco[i]:
      maxco[i] = float(atom[i+1])
  AllAtoms.append([atom[0], float(atom[1]), float(atom[2]), float(atom[3]), [] ])

BorderAtoms=AllAtoms[:]

print len(AllAtoms), "=len of all atoms"  

print "Found Box bounds [%f,%f]x[%f,%f]x[%f,%f]." % (minco[0],maxco[0],minco[1],maxco[1],minco[2],maxco[2])


input=open(sys.argv[2], "r")

Temp=[]
NumConBorderAtoms=[]
NumConBorderFaces=[]
ConBorderAtoms=[]
ConBorderFaces=[]

line = input.readline() #No relevant Info
line = input.readline() #No relevant Info
line = input.readline() 
Temp=line.split("=")    #We get the # of Atoms in Convex Border and the # of Faces in Convex Hull
NumConBorderAtoms=Temp[2].split(",")[0]
NumConBorderFaces=Temp[3].split(",")[0]

print NumConBorderAtoms, "Atoms mark the Convex Hull"
print NumConBorderFaces, "Faces mark the Convex Hull"

for i in range(int(NumConBorderAtoms)):
  line=input.readline()
  Temp=line.split()
  ConBorderAtoms.append([float(Temp[0])+67.01283, float(Temp[1])+19.49606, float(Temp[2])+36.74105])  ####This should be done differently!!!!!!!!!

line =input.readline() #This is an empty line

for i in range(int(NumConBorderFaces)):
  line=input.readline()
  Temp=line.split()
  ConBorderFaces.append([int(Temp[0])-1, int(Temp[1])-1, int(Temp[2])-1])

#print len(ConBorderAtoms)
#print len(ConBorderFaces)

MinDist=max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])
MinIndex=-1
PlaneVector=[0.0, 0.0, 0.0]

output=open("test2.xyz", "w") ##############################################Debug
output.write("%d \n \n" %(len(AllAtoms)+13))

BorderFace=[]
SumDist=[]

for j in range(int(NumConBorderFaces)):  # this calculates an orthogonal vector to the border face. It will be used to get the 
                                          # distance of an atom to the border
  betrag=0
  for l in range(3):
    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])
    betrag+=PlaneVector[l]*PlaneVector[l]

  betrag=math.sqrt(betrag)   #the vector will be normalized.
  for l in range(3):
    PlaneVector[l]/=betrag
  BorderFace.append([ConBorderFaces[j][0], PlaneVector[:]])
  SumDist.append(0.0)

"""
for j in range(len(BorderFace)):
  betrag=0
  for l in range(3):
    betrag+=BorderFace[j][1][l]*(ConBorderAtoms[ConBorderFaces[j][1]][l]-ConBorderAtoms[BorderFace[j][0]][l])
#  betrag=math.sqrt(betrag)
  print betrag
"""

i=-1
skalarp=0
MaxDist=0
MinDist=0


while 1:            # we go through all atoms and throw those away which are very close to the convex border
#  Min Dist is the minimum Distance of a single atom to all enclosing faces
  MinDist=max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])
  i+=1
  if i >= len(AllAtoms):
    break
  for j in range(len(BorderFace)):
    skalarp=0
    for l in range(3):
      skalarp+=(BorderFace[j][1][l])*(AllAtoms[i][l+1]-ConBorderAtoms[BorderFace[j][0]][l])
    if abs(skalarp)<MinDist:
      MinDist=abs(skalarp)
      AllAtoms[i][4]=[j, MinDist]
    if MinDist<3:
#      print "removing"
      output.write("H  %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
      AllAtoms.remove(AllAtoms[i])
      i=i-1
      break


MaxDist1 =-1.0
MaxDist2 =-1.0
MaxDist3 =-1.0

for j in range(len(BorderFace)):
  MaxDist1 =-1.0
  MaxDist2 =-1.0
  MaxDist3 =-1.0
  for i in range(len(AllAtoms)):
    skalarp=0
    for l in range(3):
      skalarp+=(BorderFace[j][1][l])*(AllAtoms[i][l+1]-ConBorderAtoms[BorderFace[j][0]][l])
    SumDist[j]+=abs(skalarp)

    if AllAtoms[i][4][0]==j:
      if AllAtoms[i][4][1]>MaxDist1:
        MaxDist3=MaxDist2
        MaxDist2=MaxDist1
        MaxDist1=AllAtoms[i][4][1]
      elif AllAtoms[i][4][1]>MaxDist2:
        MaxDist3=MaxDist2
        MaxDist2=AllAtoms[i][4][1]
      elif AllAtoms[i][4][1]>MaxDist3:
        MaxDist3=AllAtoms[i][4][1]
  i=-1
  while 1:
    i+=1
    if i>=len(AllAtoms):
      break
    if AllAtoms[i][4][0]==j:
      if AllAtoms[i][4][1]<MaxDist3:
        output.write("H  %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
        AllAtoms.remove(AllAtoms[i])
        i=i-1
  #print MaxDist1
  #print MaxDist2
  #print MaxDist3



  

########################################################
#
# Find direct and indirect neighbours
# Remove sucht atoms which are isolated as a Concavity
#
#
#
########################################################

Neighbours=[]
for i in range(len(AllAtoms)):
  Neighbours.append([])
for i in range(len(AllAtoms)):
  for j in range(len(AllAtoms)):
    Neighbours[i].append(-1)   # Initialized neighbour matrix. entry -1: not connected, entry 1 directly connected, entry 2 indirect connected

for i in range(len(AllAtoms)):
  for j in range(len(AllAtoms)):  # First we only look for direkt neighbours
    if distance(AllAtoms[i][1:4], AllAtoms[j][1:4])<DirectNeighbourDistance:
      Neighbours[i][j]=1
      Neighbours[j][i]=1
    
for i in range(len(AllAtoms)):
  for j in range(len(AllAtoms)):
    for k in range(len(AllAtoms)):
      if Neighbours[j][k]<0:
        for l in range(len(AllAtoms)):
          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.
            Neighbours[j][k]=2
        



i=-1     # Here we remove isolated concave atoms, we expect that a Concavity has more than one Atom.
while 1:
  i+=1
  if i >=len(AllAtoms):
    break
  MaxDist=-0.5
  for j in range(len(AllAtoms)):
    if Neighbours[i][j]>0:
      MaxDist+=1
  if MaxDist<2:
    output.write("H  %f %f %f \n" %(AllAtoms[i][1], AllAtoms[i][2], AllAtoms[i][3] ))
    for j in range(len(AllAtoms)):
      Neighbours[j][i]=5
      Neighbours[j].remove(Neighbours[j][i])
    Neighbours[i]=5
    Neighbours.remove(Neighbours[i])
    AllAtoms.remove(AllAtoms[i])
    i-=1






#######################################################
#
# We have found all those indirectly connected.
# We now sort those into Concavities
# We calculate Centers of Concavities
# We find distance to centers of concavities
# Find 10th furthest distance.
#
######################################################


ConcaveClusters=[]

for i in range(len(AllAtoms)):   # Here we create a list of Concavities and the atoms of which it consists
  if Neighbours[i][i]>0:
    ConcaveClusters.append([i])
    for j in range(i+1, len(AllAtoms)):
      if Neighbours[i][j]>0:
        ConcaveClusters[-1].append(j)
        Neighbours[j][j]=-1


print "Number of all Atoms", len(AllAtoms)
print "ConcaveClusters",ConcaveClusters
print "Number of ConcaveClusters", len(ConcaveClusters)



ConcaveCenters=[]  #this will be the center of the concave atoms
MaxDist=0
primal=[0.0, 0.0, 0.0]

for i in range(len(ConcaveClusters)):
  ConcaveCenters.append([0.0, 0.0, 0.0])
#  j=0
  for q in ConcaveClusters[i]:
    for k in range(3):
      ConcaveCenters[i][k]+=AllAtoms[q][k+1]
#      ConcaveCenters[i][k]=(ConcaveCenters[i][k]*float(j)+AllAtoms[q][k+1]) / float(j+1)
#    j+=1
  for k in range(3):
    ConcaveCenters[i][k]/=float(len(ConcaveClusters[i]))

print "ConcaveCenters", ConcaveCenters
    

MaxDist=0
j=0
for i in range(len(AllAtoms)):
  if AllAtoms[i][4][1]>20:
    for q in range(3):
      primal[q]=(primal[q]*j+AllAtoms[i][q+1])/(j+1)
    j+=1


output.write("Cs  %f %f %f \n" %(primal[0], primal[1], primal[2]))



for i in range(len(ConcaveClusters)):
  MinDist=[max(maxco[0]-minco[0], maxco[1]-minco[1], maxco[2]-minco[2])]
  for j in ConcaveClusters[i]:
    dist=distance(AllAtoms[j][1:4], ConcaveCenters[i])

    AllAtoms[j][4]=[i, dist]     #We overwrite the distance to the convex hull with distance to Concavity center. Same for respective id.
    for l in range(len(MinDist)): #We keep track of the 10 closest Atoms to Concavity Center.
      if dist<MinDist[l]:
        MinDist.insert(l, dist)
        if len(MinDist)>10:
          del MinDist[-1]
        break
  ConcaveClusters[i].append(MinDist[-1])  #respective 10th closest distance for center is saved here.


##############################################
#
# Debug
#
#
#
##############################################



for i in range(len(ConcaveClusters)):
  for j in range(len(ConcaveClusters[i])-1):
    if i ==0:
      output.write("O  %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
    if i ==1:
      output.write("Li  %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
    if i ==2:
      output.write("Fe  %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))
    if i ==3:
      output.write("He  %f %f %f \n" %(AllAtoms[ConcaveClusters[i][j]][1], AllAtoms[ConcaveClusters[i][j]][2], AllAtoms[ConcaveClusters[i][j]][3] ))


  output.write("Ne  %f %f %f \n" %(ConcaveCenters[i][0], ConcaveCenters[i][1], ConcaveCenters[i][2]))


"""

#output.write("Si %f %f %f \n" %(primal[0], primal[1], primal[2]))
#output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][0]][0], ConBorderAtoms[ConBorderFaces[Furthest][0]][1], ConBorderAtoms[ConBorderFaces[Furthest][0]][2]))
#output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][1]][0], ConBorderAtoms[ConBorderFaces[Furthest][1]][1], ConBorderAtoms[ConBorderFaces[Furthest][1]][2]))
#output.write("Ca %f %f %f \n" %(ConBorderAtoms[ConBorderFaces[Furthest][2]][0], ConBorderAtoms[ConBorderFaces[Furthest][2]][1], ConBorderAtoms[ConBorderFaces[Furthest][2]][2]))

#output.close()
"""


####################################################
#
# Until here we removed convex atoms. All remaining atoms 
# belong to some kind of concavity
# We will now put them into the kartesian planes in polar coordinates and assign a degree of concavity:
# max(len1*cos(alpha), len2*cos(beta)) / lenMid
# Sort bei degree of concavity.
# Test over best 1*10*20
#
####################################################
"""

def MergeSort(List, Typus):
  """#  applies merge sort to list given
"""
  if len(List)<2:
    return List
  elif Typus==1:
    List1=List[:int(len(List)/2)]
    List2=List[int(len(List)/2):]
    return(MergeAngle(MergeSort(List1, 1), MergeSort(List2, 1)))
  elif Typus==2:
    List1=List[:int(len(List)/2)]
    List2=List[int(len(List)/2):]
    return(MergeKrit(MergeSort(List1, 2), MergeSort(List2, 2)))
    

def MergeAngle(List1, List2):
  """# merges two lists in ordered way
"""
  List=[]
  while len(List1)>0 and len(List2)>0:
    if List1[0][4][0][1]<List2[0][4][0][1]:
      List.append(List1[0])
      List1.remove(List1[0])
    else:
      List.append(List2[0])
      List2.remove(List2[0])
  while len(List1)>0:
    List.append(List1[0])
    List1.remove(List1[0])
  while len(List2)>0:
    List.append(List2[0])
    List2.remove(List2[0])
  return(List)




def MergeKrit(List1, List2):
  """# merges two lists in ordered way
"""
  List=[]
  while len(List1)>0 and len(List2)>0:
    if List1[0][4][1]>List2[0][4][1]:
      List.append(List1[0])
      List1.remove(List1[0])
    else:
      List.append(List2[0])
      List2.remove(List2[0])
  while len(List1)>0:
    List.append(List1[0])
    List1.remove(List1[0])
  while len(List2)>0:
    List.append(List2[0])
    List2.remove(List2[0])
  return(List)
    

r=0.0
phi=0.0
sign=0.0
Krit=0.0
TempKrit=0.0

for i in range(2):
  for j in range(i+1,3):
    for k in range(len(AllAtoms)):
#      print AllAtoms[k][i+1], primal[i]
#      print AllAtoms[k][j+1], primal[j]
      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)
      if AllAtoms[k][j+1]-primal[j]>0:
        phi=1
      else:
        phi=-1
      phi*=math.acos((AllAtoms[k][i+1]-primal[i])/r)  # phi = sign(y)*arccos(x/r)
      AllAtoms[k][4]=[[r, phi], Krit]
    
    AllAtoms=MergeSort(AllAtoms, 1)

    for k in range(len(AllAtoms)):
#      print AllAtoms[k][4][0][1]
      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]
      # Assign the ratio of the projected length of two outer atom legs to the inner leg. Krit > 1 is concavity
      if Krit>AllAtoms[k-1][4][1]:
        AllAtoms[k-1][4][1]=Krit

AllAtoms=MergeSort(AllAtoms, 2)

#for i in range(len(AllAtoms)):
#  print AllAtoms[i][4][1]





"""


####################################################
#
# 
# 
# We will now iterate over all triples and check the bonds
# being cut by the respective plane
#
####################################################


CompleteAtoms =[]
CompleteBonds =[]

input = open("silica.grape.%.4d" %(int(sys.argv[3])), "r")
line=input.readline()

for line in input:
  CompleteAtoms.append([float(line.split()[1]), float(line.split()[2]), float(line.split()[3]) ])
input.close()

input=open("silica.dbond-SiOCaO.%.4d" %(int(sys.argv[3])), "r")
line=input.readline()

for line in input:
  CompleteBonds.append([int(line.split()[0]), int(line.split()[1])])
input.close()


CenterOfTriangle=[0.0, 0.0, 0.0]
MinimumPlane=[0, 0, 0]
SecondPlane=[0,0,0]
ThirdPlane=[0,0,0]
MinimumCutCount=len(CompleteBonds)
SecondCutCount=len(CompleteBonds)
ThirdCutCount=len(CompleteBonds)
CutCount=0
PlanePreCheck=0

PlaneVector=[0.0, 0.0, 0.0]
Atom1Vector=[0.0, 0.0, 0.0]
Atom2Vector=[0.0, 0.0, 0.0]
DirectionVector=[0.0, 0.0, 0.0]


iteration=0


for i in range(len(ConcaveCenters)):
  for j in range(i, len(ConcaveCenters)):
    for k in ConcaveClusters[i]+ConcaveClusters[j]:
      if k!=int(k) or AllAtoms[k][4][1]>ConcaveClusters[AllAtoms[k][4][0]][-1]:
#      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]):

# (Neighbours[i][j]>0 and Neighbours[i][k]>0)
# or  not(AllAtoms[i][4][1]>20 or AllAtoms[j][4][1]>20 or AllAtoms[k][4][1]>20) or
# (AllAtoms[i][4][1]>20 and AllAtoms[j][4][1]>20 and AllAtoms[k][4][1]>20): 
# or Neighbours[i][k]==1 or Neighbours[j][k]==1
        continue

      for q in range(3):
        CenterOfTriangle[q]= (ConcaveCenters[i][q] + ConcaveCenters[j][q] + AllAtoms[k][q+1])/3.0
      
      
      ### this calculates the crossproduct of the vectors in the plane, yielding an orthogonal vector on the plane
      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])
      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])
      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])
      
      ### cutting plane is done, next plane will ensure, that we only cut one leg of the cluster

      for q in range(3):
        DirectionVector[q] = CenterOfTriangle[q]-ConcaveCenters[i][q]

      """
      skalarpA=0.0
      skalarpB=0.0
      skalarpC=0.0
      for q in range(3):
        skalarpA += (ConcaveCenters[i][q]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[k][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
        skalarpB += (ConcaveCenters[i][hust][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[j][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
        skalarpC += (AllAtoms[j][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q]) * (AllAtoms[k][q+1]-ConcaveCenters[ConcaveCenters[i][hust][4][0]][q])
      if skalarpA<0 or skalarpB<0 or skalarpC<0:
        continue  #Concavity seperates triangle
      """
      
      PlanePreCheck=0
      for l in range(len(BorderAtoms)):
        skalarp=0
        for g in range(1,4):
          skalarp+=(BorderAtoms[l][g]-ConcaveCenters[i][g-1])*PlaneVector[g-1]
        if skalarp<0:
          PlanePreCheck+=1
      if PlanePreCheck<(4.0/5.0*len(BorderAtoms)) and PlanePreCheck>(len(BorderAtoms)/5.0):
        CutCount=0
      else:
        continue
      print "iteration", iteration
      iteration+=1

      for l in range(len(CompleteBonds)):  ### Go through all the bonds, create vectors from first point in plane to both atoms
        if CutCount>=ThirdCutCount:       ### This gives small speed enhancement as expensive planes will not be fully computed anymore.
          break

        Skalar1=0
        Skalar2=0
        # We again work with the sign of the skalar product
        """
        for q in range(3):                   #Calculate Vectors from Concavity to Atoms
          Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[i][q]
          Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[i][q]

        for q in range(3):
          Skalar1+=Atom1Vector[q]*DirectionVector[q]
        if Skalar1<0:
          continue
        Skalar1=0
        """
        for q in range(3):                   #Calculate Vectors from plane to Atoms
          Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[i][q]
          Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[i][q]

        for q in range(3):
          Skalar1+=Atom1Vector[q]*PlaneVector[q]
          Skalar2+=Atom2Vector[q]*PlaneVector[q]
        if Skalar1*Skalar2<0:
#          print Atom1Vector
#          print Atom2Vector
          CutCount+=1
        

      if CutCount<ThirdCutCount and CutCount>50:
        if CutCount<SecondCutCount:
          if CutCount<MinimumCutCount:
            ThirdPlane=SecondPlane[:]
            ThirdCutCount=SecondCutCount
            SecondPlane=MinimumPlane[:]
            SecondCutCount=MinimumCutCount
            MinimumCutCount=CutCount
            MinimumPlane=[i, j, k]
            print PlanePreCheck
            print MinimumCutCount
            print MinimumPlane
          else:
            ThirdPlane=SecondPlane[:]
            ThirdCutCount=SecondCutCount
            SecondPlane=[i, j, k]
            SecondCutCount=CutCount
        else:
          ThirdPlane=[i, j, k]
          ThirdCutCount=CutCount
          


print MinimumCutCount
print MinimumPlane
print AllAtoms[MinimumPlane[0]]
print AllAtoms[MinimumPlane[1]]
print AllAtoms[MinimumPlane[2]]


print SecondCutCount
print SecondPlane
print AllAtoms[SecondPlane[0]]
print AllAtoms[SecondPlane[1]]
print AllAtoms[SecondPlane[2]]
print ThirdCutCount
print ThirdPlane
print AllAtoms[ThirdPlane[0]]
print AllAtoms[ThirdPlane[1]]
print AllAtoms[ThirdPlane[2]]


#for i in range(len(Center)):
#  output.write("Ca %f %f %f \n" %(Center[i][0], Center[i][1], Center[i][2]))

for i in range(2):
  output.write("Si %f %f %f \n" %(ConcaveCenters[MinimumPlane[i]][0], ConcaveCenters[MinimumPlane[i]][1], ConcaveCenters[MinimumPlane[i]][2]))
  output.write("Ar %f %f %f \n" %(ConcaveCenters[SecondPlane[i]][0], ConcaveCenters[SecondPlane[i]][1], ConcaveCenters[SecondPlane[i]][2]))
  output.write("Ca %f %f %f \n" %(ConcaveCenters[ThirdPlane[i]][0], ConcaveCenters[ThirdPlane[i]][1], ConcaveCenters[ThirdPlane[i]][2]))

output.write("Si %f %f %f \n" %(AllAtoms[MinimumPlane[2]][1], AllAtoms[MinimumPlane[2]][2], AllAtoms[MinimumPlane[2]][3]))
output.write("Ar %f %f %f \n" %(AllAtoms[SecondPlane[2]][1], AllAtoms[SecondPlane[2]][2], AllAtoms[SecondPlane[2]][3]))
output.write("Ca %f %f %f \n" %(AllAtoms[ThirdPlane[2]][1], AllAtoms[ThirdPlane[2]][2], AllAtoms[ThirdPlane[2]][3]))

output.close()


########################################################
#
# Having found a cutting plane we will now
# seperate the Cluster, first cutting all bonds as 
# they are cut by the plane, the cutting the Atoms 
# which were the resting points of the plane
#
########################################################

output=open("silica.dbond-SiOCaO.Cut", "w")

output.write("0 666666 \n")


      ### this calculates the crossproduct of the vectors in the plane, yielding an orthogonal vector on the plane
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])

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])

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])
      
      ### cutting plane is done, next plane will ensure, that we only cut one leg of the cluster
absolut=0
for q in range(3):
  absolut+=PlaneVector[q]*PlaneVector[q]
absolut=math.sqrt(absolut)

for q in range(3):
  PlaneVector[q]/=absolut

print PlaneVector

for q in range(3):
  CenterOfTriangle[q]= (ConcaveCenters[MinimumPlane[0]][q] + ConcaveCenters[MinimumPlane[1]][q] + AllAtoms[MinimumPlane[2]][q+1])/3.0
  DirectionVector[q] = CenterOfTriangle[q]-ConcaveCenters[MinimumPlane[0]][q]


l=-1
  #for l in range(len(CompleteBonds)):  ### Go through all the bonds, create vectors from first point in plane to both atoms
while 1:
    l+=1
    if l>=len(CompleteBonds):
      break
    Skalar1=0
    Skalar2=0
    # We again work with the sign of the skalar product
    """
    for q in range(3):                   #Calculate Vectors from Concavity to Atoms
      Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]

    for q in range(3):
      Skalar1+=Atom1Vector[q]*DirectionVector[q]
    
    if Skalar1<0:
      output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))
      del CompleteBonds[l]
      l-=1
      continue

    Skalar1=0
    """
    for q in range(3):                   #Calculate Vectors from plane to Atoms
      Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]
      Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]
    
   # if distance(Atom2Vector[:], [0, 0, 0])<1:
   #   print Atom1Vector
    for q in range(3):
      Skalar1+=Atom1Vector[q]*PlaneVector[q]
      Skalar2+=Atom2Vector[q]*PlaneVector[q]
    if abs(Skalar1)>eps and abs(Skalar2)>eps:
      if Skalar1*Skalar2>0:
        output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))
        del CompleteBonds[l]
        l-=1
      else:
        del CompleteBonds[l]
        l-=1
    else:
      print Skalar1, Skalar2


######################
#
# Now we make small parallel shift and take care of the remaining bonds
#
#
######################

ShiftCount=[0, 0]
shift=10**-2

for i in range(2):
  for l in range(len(CompleteBonds)):
    for q in range(3):                   #Calculate Vectors from plane to Atoms
      Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
      Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
    Skalar1=0
    Skalar2=0
    for q in range(3):
      Skalar1+=Atom1Vector[q]*PlaneVector[q]
      Skalar2+=Atom2Vector[q]*PlaneVector[q]
    if Skalar1*Skalar2<0:
      ShiftCount[i]+=1

if ShiftCount[0]>ShiftCount[1]:
  i=1
else:
  i=0


for l in range(len(CompleteBonds)):
    for q in range(3):                   #Calculate Vectors from plane to Atoms
      Atom1Vector[q]=CompleteAtoms[CompleteBonds[l][0]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
      Atom2Vector[q]=CompleteAtoms[CompleteBonds[l][1]-1][q]-ConcaveCenters[MinimumPlane[0]][q]+(-1**i)*shift*PlaneVector[q]
    Skalar1=0
    Skalar2=0
    for q in range(3):
      Skalar1+=Atom1Vector[q]*PlaneVector[q]
      Skalar2+=Atom2Vector[q]*PlaneVector[q]
    if Skalar1*Skalar2>0:
      output.write("%d %d \n" %(CompleteBonds[l][0], CompleteBonds[l][1]))



print len(CompleteBonds)