'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' ACD CHEMBASIC DEMO PROGRAM                                          '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'                                                                     '
' Molecular 3D Editor//POLYHEDR.BAS                                   '
'                                                                     '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'                                                                     '
' The utility measures Porai-Koshits/Muetterties parameters of        '
' the molecular polyhedron                                            '
'                                                                     '
' NOTE: a special attention is paid to 8-coordination                 '
'                                                                     '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''




CONST TITLE="ChemBasic Molecular Editor // Polyhedron"
CONST TITLETBOXCONTENT="Porai-Koshits/Muetterties analysis of polyhedral shapes"
CONST RAD_TO_DEG = 57.29577951
CONST MAXDIHEDR = 200           ' max number of dihedrals
CONST NVERT8A = 4               ' number of type A vertices
CONST NVERT8B = 4               ' number of type B vertices
CONST L0 = 200                  ' left
CONST U0 = 200                  ' upper corner
CONST DOUBLE_ZERO = 1E-300      ' min value at double arithmetics


Dim nvert As Integer             ' number of vertices
Dim Pi As Double
Dim one As Double
Dim zero As Double
Dim rspage As Object            ' results page
Dim CR As String                ' <CR>







'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Main As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' POLYHEDR.BAS                                                        '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim naround As Integer, OK As Boolean
Dim page,diag,asm,struc,atcenter,ataround() As Object
  MAIN="Failed or nothing to do!"

  CR = Chr(13) & Chr(10)

  ' Get 1st structure from the current page
  page=ActiveDocument.ActivePage
  If page.Diagrams.Count<1 Then Exit Function
  asm=Assemblies.AddFromCS(page.Diagrams.Item(1))
  If asm=NULL Then Exit Function
  struc=asm.Structures.Item(1)
  If struc=NULL Then Exit Function
  diag = CreateACopy(struc)
  If diag=NULL Then Exit Function

  ' Do the job
  OK=LabelDiagramWithNumbers(diag)
  If Not OK Then Exit Function
  OK=ParsePolyhedron(struc,atcenter,naround,ataround)
  If Not OK Then Exit Function
  Call DihedralAngles(diag,struc,naround,ataround)

  Main="Completed."
End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function CreateACopy(struc As Object) As Object
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Create a copy of diagram on a newly added page                      '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim l,u,w,h As Integer, normsz As Double, diag As Object
  CreateACopy=NULL

  rspage=ActiveDocument.AddEmpty
  diag=rspage.Diagrams.AddEmpty
  If diag=NULL Then Exit Function
  With diag
    .Depict (struc)
    .GetBound(l,u,w,h)
    normsz = 1
    If (w > 1200 - L0) Then normsz = (1200 - L0) / w
    w = Fix(w * normsz)
    h = Fix(h * normsz)
    If (h > 1400 - (U0 + 150)) Then normsz = (1400 - (U0 + 150)) / h
    w = Fix(w * normsz)
    h = Fix(h * normsz)
    .SetBound(L0,U0+150,w,h)
  End With

  CreateACopy=diag
End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function ParsePolyhedron(struc As Object,atcenter As Object,naround As Integer,ataround() As Object) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Find polyhedron's centre and connected atoms                        '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim icenter,nat,k,nn,maxnn As Integer
Dim sresp,TitleTBoxCont As String, OK As Boolean
Dim asm,mol,at,candid,around,TitleTBox As Object
  ParsePolyhedron=FALSE

  asm=struc.Assembly
  nat=asm.Count
  mol=struc.Molecule

  ' Select candidate for central atom of polyhedron
  maxnn = 0
  For Each at In asm
    nn=mol.AssocAtoms(at).Count
    If (maxnn < nn) Then
      maxnn = nn
      candid = at
    End If
  Next at
  icenter=asm.Index(candid)

  ' Ask for the centre
  sresp = UserIOBox("Please supply a central atom of polyhedron", TITLE, Str(icenter))
  icenter = Fix(val(sresp))

  ' Check and parse the results
  If (icenter<1) Or (icenter>nat) Then Exit Function
  atcenter = mol.Assembly.Item(icenter)
  around = mol.AssocAtoms(atcenter)
  naround = around.Count
  ReDim ataround(naround) As Object
  For k = 1 To naround
    ataround(k) = around.Item(k)
  Next k
  ' Printouts
  TitleTBox = rspage.TextBoxes.AddEmpty
  TitleTBox.SetContent (TITLETBOXCONTENT)
  TitleTBox.SetBound(L0,U0,800,100)

  ParsePolyhedron=TRUE
End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub DihedralAngles(diag As Object, struc As Object,nat As Integer, ats() As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,j,k,l, u, w, h, atnum(),nedges,nsteps,iedge() As Integer
Dim dihedr(MAXDIHEDR),x(),y(),z(),normsz As Double
Dim asm,mol,at,at1,diagP,TitleTBox As Object
Dim IsInPoly As Boolean

  ' Prepare to work
  CR = Chr(13) & Chr(10)
  one = 1 : zero = 0 : Pi = 4 * Atan(1.0)
  asm=struc.Assembly
  mol=struc.Molecule
  For i = 1 To asm.Count
    Struc.Assembly.Item(i).SetName(Str(i))
  Next i
  nvert = nat
  ReDim x(nvert),y(nvert),z(nvert) As Double      ' atomic(vertices) coordinates
  ReDim atnum(nvert) As Integer                   ' remember atomic numbers
  ReDim iedge(nvert, nvert) As Integer            ' contains numbers of dihedral angles
                                                  ' in the order they were found
  For i = 1 To nvert
    at = ats(i)
    struc.GetAtomXYZ(at,x(i),y(i),z(i))
    atnum(i) = Fix(Val(at.GetName))
  Next i

  ' Find a convex polyhedron which consists of nvert vertices
  Print "Finding convex polyhedron..."
  Print "The following dihedral angles found:"
  nedges = FindConvexPolyhedron(nsteps, iedge, dihedr, ats, atnum, x, y, z, struc)
  Print "Considered: ", nsteps, " possible dihedrals;  Found: ", nedges, " present"
  Print "Done."

  ' Prepare the polyhedron for drawing
  If (nvert = 0) Then Exit Sub
  If (nvert = 2) Or (nvert = 3) Then
    iedge(1, 2) = 1 : iedge(2, 1) = 1
  End If
  If (nvert = 3) Then
    iedge(1, 3) = 2 : iedge(3, 1) = 2 : iedge(2, 3) = 3 : iedge(3, 2) = 3
  End If
  For k = 1 To asm.Count
    at = asm.Item(k)
    IsInPoly = False
    For i = 1 To nvert
    If (Str(atnum(i)) = at.GetName) Then
      IsInPoly = True
      Exit For
    End If
    Next i
    If (IsInPoly) Then
      ' Add bonds
      For j = 1 To nvert
        at1 = ats(j)
        If (iedge(i, j) > 0) Then mol.AddBond(at, at1, 1)
      Next j
    Else
      Kill (at)
      k = k - 1
    End If
  Next k

  ' Draw the polyhedron
  diagP = rspage.Diagrams.AddEmpty
  diagP.Depict(struc)
  diag.GetBound(l,u,w,h)
  normsz=1.0
  If (w > 1600 - 1200) Then normsz = (1600 - 1200) / w
  w = Fix(w * normsz)
  h = Fix(h * normsz)
  If (h > 700 - 200) Then normsz = (700 - 200) / h
  w = Fix(w * normsz)
  h = Fix(h * normsz)
  diagP.SetBound(1200,200,w,h)

  ' Analyze the polyhedron
  Select Case nvert
    Case 0
      Print "Number of atoms = 0"
    Case 1
      Print "There is no need in definition of the proximate polyhedron for the 1-atom family"
    Case 2
      Print "There is no need in definition of the proximate polyhedron for the 2-atoms family"
    Case 3
      Print "There is no need in definition of the proximate polyhedron for the 3-atoms family"
    Case 4
      Print "Going to define the proximate polyhedron"
      Print "for the 4-atoms tetrahedron family"
      Call Proceed4Coord(iedge, dihedr, atnum, x, y, z)
    Case 5
      Print "Going to define the proximate polyhedron"
      Print "for the 5-atoms trigonal bipyramid family"
      Call Proceed5Coord(iedge, dihedr, atnum, x, y, z)
    Case 6
      Print "Going to define the proximate polyhedron"
      Print "for the 6-atoms octahedron family"
      Call Proceed6Coord(iedge, dihedr, atnum, x, y, z)
    Case 7
      Print "Going to define the proximate polyhedron"
      Print "for the 7-atoms capped octahedron family"
      Call Proceed7Coord(iedge, dihedr, atnum, x, y, z)
    Case 8
      Print "Going to define the proximate polyhedron"
      Print "for the 8-atoms dodecahedron family"
      Call Proceed8Coord(iedge, dihedr, atnum, x, y, z)
    Case 9
      Print "Definition of the proximate polyhedron is not realized yet for 9-atoms family"
    Case 10
      Print "Definition of the proximate polyhedron is not realized yet for 10-atoms family"
    Case Else
      Print "Definition of the proximate polyhedron is not realized yet for more than 10-atoms family"
  End Select 'nvert


End Sub 'DihedralAngles





'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function FindConvexPolyhedron(nsteps As Integer, iedge() As Integer, dihedr() As Double, ats() As Object, atnum() As Integer, x() As Double, y() As Double, z() As Double, struc As Object) As Integer
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Find a convex polyhedron for given atoms                            '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim j,i,i1,i2,i3,i4,i5,CurDh,DhAngBoxHt As Integer
Dim A1,B1,C1,d1,A2,B2,C2,d2,xi,yi,zi,ang,infsmall As Double
Dim messg, DhAngBoxCont As String, IsConvex,ReDefX As Boolean
Dim at1,at2,at3,at4,DhAngBox As Object

  infsmall = 0.0001
  ReDefX = FALSE
  For i = 1 To nvert
    For j = 1 To nvert
      iedge(i, j) = 0
    Next j
  Next i
  CurDh = 0
  nsteps = 0
  DhAngBoxCont = "Dihedral angles list between atoms:"
  DhAngBox = rspage.TextBoxes.AddEmpty

  ' Analyze now
  For i1 = 1 To nvert - 1
    For i2 = i1 + 1 To nvert
      For i3 = 1 To nvert - 1
        If (i3 <> i1) And (i3 <> i2) Then
          For i4 = i3 + 1 To nvert
            If (i4 <> i1) And (i4 <> i2) Then
              nsteps = nsteps + 1
              IsConvex = TRUE
              Call GetPlane(A1, B1, C1, d1, x(i1), x(i2), x(i3), y(i1), y(i2), y(i3), z(i1), z(i2), z(i3))
              Call GetPlane(A2, B2, C2, d2, x(i1), x(i2), x(i4), y(i1), y(i2), y(i4), z(i1), z(i2), z(i4))
              ang = (PI - Abs(struc.GetTAngle(ats(i3), ats(i1), ats(i2), ats(i4)))) * 180 / Pi
              ' If an edge between i3, i4 already exists and dihedral angle is 0.0
              ' then edge between i1, i2 isn't possible (crossed edges in planar dihedron)
              If (ang <= zero) And (iedge(i3, i4) > 0) Then IsConvex = False
              If (ang <= zero) Or (ang >= Pi * 180 / Pi) Then
                x(i4) = x(i4) + infsmall
                y(i4) = y(i4) + infsmall
                z(i4) = z(i4) + infsmall
                ' Note: if all vertices lie on the one side with NEW i4 from the plane i1-i2-i3
                ' then those vertices lie on the one side with i3 from ihe plane i1-i2-i4
                ReDefX = TRUE
              End If
              For i5 = 1 To nvert
                If (i5 <> i1) And (i5 <> i2) And (i5 <> i3) And (i5 <> i4) Then
                  Call GetIntersPoint(xi, yi, zi, A1, B1, C1, d1, x, y, z, i4, i5)
                  If (mini(x(i4), x(i5)) <= xi) And (xi <= maxi(x(i4), x(i5))) And (mini(y(i4), y(i5)) <= yi) And (yi <= maxi(y(i4), y(i5))) And (mini(z(i4), z(i5)) <= zi) And (zi <= maxi(z(i4), z(i5))) Then IsConvex = False
                  Call GetIntersPoint(xi, yi, zi, A2, B2, C2, d2, x, y, z, i3, i5)
                  If (mini(x(i3), x(i5)) <= xi) And (xi <= maxi(x(i3), x(i5))) And (mini(y(i3), y(i5)) <= yi) And (yi <= maxi(y(i3), y(i5))) And (mini(z(i3), z(i5)) <= zi) And (zi <= maxi(z(i3), z(i5))) Then IsConvex = False
                End If 'i5
                If (IsConvex = False) Then Exit For 'i5
              Next i5
              If (ReDefX) Then
                ' This is a weak point when i1-i2-i3-i4 lie in one plane
                ' it is possible not to find the dihedral in that plane
                ' (the probability is very small though)
                x(i4) = x(i4) - infsmall
                y(i4) = y(i4) - infsmall
                z(i4) = z(i4) - infsmall
                ReDefX = FALSE
              End If
              If (IsConvex) Then
                CurDh = CurDh + 1
                dihedr(CurDh) = ang
                iedge(i1, i2) = CurDh
                iedge(i2, i1) = CurDh
                messg = Str(atnum(i1)) & "," & Str(atnum(i2)) & "," & Str(atnum(i3)) & " and " & Str(atnum(i1)) & "," & Str(atnum(i2)) & "," & Str(atnum(i4)) & " = " & FStr(dihedr(CurDh), 2, 1) & " deg."
                DhAngBoxCont = DhAngBoxCont & CR & messg
                DhAngBoxHt = CurDh * 50 + 50
                DhAngBox.SetContent (DhAngBoxCont)
                DhAngBox.SetBound(L0,1400,800,DhAngBoxHt)
                Print CurDh, ": between atoms ", messg
              End If 'IsConvex
            End If 'i4
          Next i4
        End If 'i3
      Next i3
    Next i2
  Next i1

  FindConvexPolyhedron = CurDh
End Function 'FindConvexPolyhedron




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Proceed8Coord(iedge() As Integer, dihedr() As Double, atnum() As Integer, x() As Double, y() As Double, z() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Treat 8-coordination shapes                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i, j As Integer
Dim A(NVERT8A) As Integer       ' vertices of type A
Dim b(NVERT8B) As Integer       ' vertices of type B
Dim PhiBox As Object            ' TextBox containing phi angles
Dim PhiBoxCont As String        ' TextBox contents with phi angles
Dim AllPhiBox As Object         ' TextBox containing phi angles for all polyhedra
Dim AllPhiBoxCont As String     ' TextBox contents with phi angles for all polyhedra
Dim phi1, phi2 As Double        ' phi1, phi2 angles
Dim IsDodFam As Boolean         ' is from dodecahedron family
Dim dt(4) As Double             ' delta angles
Dim Dev(3) As Double            ' square deviation from ideal polyhedra
Dim Poly(3) As String           ' possible polyhedron
Dim StdAng(3, 5) As Double      ' standart phi, delta1-4 for dod, bctp, sap
                                ' the foolowing is neseccary in any  ProceedNCoord
Dim DeltaBox As Object          ' TextBox containing delta angles
Dim DeltaBoxCont As String      ' TextBox contents with delta angles
Dim PolyhBox As Object          ' TextBox containing polyhedron information
Dim PolyhBoxCont As String      ' TextBox contents with polyhedron information
Dim AllDeltaBox As Object       ' TextBox containing delta angles for all polyhedra
Dim AllDeltaBoxCont As String   ' TextBox contents with delta angles for all polyhedra
Const PhiWt = 1                 ' phi weight to the deviation

  ' Standard shapes
  StdAng(1, 1) = 0
  StdAng(1, 2) = 29.5
  StdAng(1, 3) = 29.5
  StdAng(1, 4) = 29.5
  StdAng(1, 5) = 29.5
  Poly(1) = "D2d-DODECAHEDRON"
  StdAng(2, 1) = 16.1
  StdAng(2, 2) = 0
  StdAng(2, 3) = 21.7
  StdAng(2, 4) = 48.2
  StdAng(2, 5) = 48.2
  Poly(2) = "C2v-BICAPPED TRIGONAL PRISM"
  StdAng(3, 1) = 24.5
  StdAng(3, 2) = 0
  StdAng(3, 3) = 0
  StdAng(3, 4) = 52.5
  StdAng(3, 5) = 52.5
  Poly(3) = "D4d-TETRAGONAL ANTIPRISM"

  ' Separating all vertices to A and B types following Porai-Koshits
  Print "Separating vertices by A and B type..."
  IsDodFam = SeparateVertToAB(A, b, iedge)
  Print "Done."

  'For dodecahedral shapes
  If (IsDodFam) Then

    ' Reorder A, B according to pic.3 in of Porai-Koshits' paper
    Print "Reordering vertices A and B..."
    Call ReorderAB(A, b, iedge, dihedr)
    Print "Done."

    ' Show deltas
    dt(1) = dihedr(iedge(b(1), b(2)))
    dt(2) = dihedr(iedge(b(3), b(4)))
    dt(3) = dihedr(iedge(b(2), b(3)))
    dt(4) = dihedr(iedge(b(4), b(1)))
    DeltaBoxCont = "Dihedral angles delta belonging to edges of type b:"
    DeltaBoxCont = DeltaBoxCont & CR & "delta1 (belonging to the edge " & Str(atnum(b(1))) & "-" & Str(atnum(b(2))) & ") = " & FStr(dt(1), 2, 1) & " deg."
    DeltaBoxCont = DeltaBoxCont & CR & "delta2 (belonging to the edge " & Str(atnum(b(3))) & "-" & Str(atnum(b(4))) & ") = " & FStr(dt(2), 2, 1) & " deg."
    DeltaBoxCont = DeltaBoxCont & CR & "delta3 (belonging to the edge " & Str(atnum(b(2))) & "-" & Str(atnum(b(3))) & ") = " & FStr(dt(3), 2, 1) & " deg."
    DeltaBoxCont = DeltaBoxCont & CR & "delta4 (belonging to the edge " & Str(atnum(b(4))) & "-" & Str(atnum(b(1))) & ") = " & FStr(dt(4), 2, 1) & " deg."
    DeltaBox = rspage.TextBoxes.AddEmpty
    DeltaBox.SetContent (DeltaBoxCont)
    DeltaBox.SetBound(1200,700,1000,250)

    ' Get phi between A2-A4--B4-B2 and A1-A3--B3-B1
    Print "Calculating phi angles..."
    phi1 = GetPhi(A(2), A(4), b(4), b(2), x, y, z) * 180 / Pi
    phi2 = GetPhi(A(1), A(3), b(3), b(1), x, y, z) * 180 / Pi
    Print "Done."

    PhiBoxCont = "Phi angles in diagonal trapezoids:"
    PhiBoxCont = PhiBoxCont & CR & "phi1 (in the trapezoid " & Str(atnum(b(1))) & "-" & Str(atnum(A(1))) & "-" & Str(atnum(A(3))) & "-" & Str(atnum(b(3))) & ") = " & FStr(phi1, 2, 1) & " deg."
    PhiBoxCont = PhiBoxCont & CR & "phi2 (in the trapezoid " & Str(atnum(b(2))) & "-" & Str(atnum(A(2))) & "-" & Str(atnum(A(4))) & "-" & Str(atnum(b(4))) & ") = " & FStr(phi2, 2, 1) & " deg."
    PhiBox = rspage.TextBoxes.AddEmpty
    PhiBox.SetContent (PhiBoxCont)
    PhiBox.SetBound(1200,1000,1000,150)

    AllDeltaBoxCont = "Delta angles in the ideal polyhedra:" & CR
    For i = 1 To 3
      AllDeltaBoxCont = AllDeltaBoxCont & Poly(i) & ":" & CR
      For j = 2 To 5
        AllDeltaBoxCont = AllDeltaBoxCont & FStr(StdAng(i, j), 2, 1)
        If (j < 5) Then AllDeltaBoxCont = AllDeltaBoxCont & ", " Else AllDeltaBoxCont = AllDeltaBoxCont & " deg." & CR
      Next j
    Next i
    AllDeltaBox = rspage.TextBoxes.AddEmpty
    AllDeltaBox.SetContent (AllDeltaBoxCont)
    AllDeltaBox.SetBound(1200,1400,1000,350)

    AllPhiBoxCont = "Phi angles in the ideal polyhedra:" & CR
    For i = 1 To 3
      AllPhiBoxCont = AllPhiBoxCont & Poly(i) & ": " & FStr(StdAng(i, 1), 2, 1) & " deg." & CR
    Next i
    AllPhiBox = rspage.TextBoxes.AddEmpty
    AllPhiBox.SetContent (AllPhiBoxCont)
    AllPhiBox.SetBound(1200,1850,1000,350)
    For i = 1 To 3
      Dev(i) = PhiWt * (Square(phi1 - StdAng(i, 1)) + Square(phi2 - StdAng(i, 1))) / (4 + 2 - 1)
      For j = 2 To 5
        Dev(i) = Dev(i) + Square(dt(j - 1) - StdAng(i, j)) / (4 + 2 - 1)
      Next j
      Print Dev(i)
    Next i
    PolyhBoxCont = "The proximate polyhedron is" & CR
    If (Dev(1) <= Dev(2)) And (Dev(1) <= Dev(3)) Then PolyhBoxCont = PolyhBoxCont & Poly(1)
    If (Dev(2) <= Dev(1)) And (Dev(2) <= Dev(3)) Then PolyhBoxCont = PolyhBoxCont & Poly(2)
    If (Dev(3) <= Dev(1)) And (Dev(3) <= Dev(2)) Then PolyhBoxCont = PolyhBoxCont & Poly(3)


  Else   'IsDodFam


    PolyhBoxCont = "This polyhedron does not belong" & CR & "to the dodecahedron family"


  End If 'IsDodFam


  PolyhBox = rspage.TextBoxes.AddEmpty
  PolyhBox.SetContent (PolyhBoxCont)
  PolyhBox.SetBound(1600,500,700,100)

End Sub 'Proceed8Coord



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function SeparateVertToAB(A() As Integer, b() As Integer, iedge() As Integer) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Classify 8-coordination shape vertces to A or B type                '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,j,n_edges,nB,nA As Integer
  nA = 0
  nB = 0
  SeparateVertToAB = TRUE

  For i = 1 To nvert
    n_edges = 0
    For j = 1 To nvert
      If (iedge(i, j) > 0) Then n_edges = n_edges + 1
    Next j
    If (n_edges = 5) Then
      nB = nB + 1
      If (nB > 4) Then
        SeparateVertToAB = FALSE
        Exit Function
      End If
    b(nB) = i
    Else
      nA = nA + 1
      If (nA > 4) Then
        SeparateVertToAB = FALSE
        Exit Function
      End If
      A(nA) = i
    End If 'n_edges
  Next i

  If (nA <> 4) Or (nB <> 4) Then SeparateVertToAB = FALSE
End Function 'SeparateVertToAB


'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub ReorderAB(A() As Integer, b() As Integer, iedge() As Integer, dihedr() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,j,k,l As Integer, min_dih As Double
Dim AA(NVERT8A) As Integer
  ' Reordering B to a chain (B1-B2, B2-B3, B3-B4, B4-B1)
  ' finding B1, B2 as having minimal dihedral angle;
  ' B3 as connected to B2; B4 as left one
  min_dih = 2 * Pi * 180 / PI
  For i = 1 To NVERT8B - 1
    For j = i + 1 To NVERT8B
      If (iedge(b(i), b(j)) > 0) Then
        If (dihedr(iedge(b(i), b(j))) < min_dih) Then
          k = i
          l = j
          min_dih = dihedr(iedge(b(i), b(j)))
        End If ' dihedr
      End If ' iedge
    Next j
  Next i
  Call swap(b(k), b(1))
  Call swap(b(l), b(2))
  If (iedge(b(2), b(3)) = 0) Then Call swap(b(3), b(4)) 'swap B3, B4

  ' Find A1 as connected to B1, B2, B4;
  '         A2 as connected to B1, B2, B3;
  '         A3 as connected to B2, B3, B4;
  '         A4 as connected to B1, B3, B4.
  For i = 1 To NVERT8A
    AA(i) = A(i)
  Next i
  For i = 1 To NVERT8A
    If (iedge(AA(i), b(1)) > 0) And (iedge(AA(i), b(2)) > 0) And (iedge(AA(i), b(4)) > 0) Then A(1) = AA(i)
    If (iedge(AA(i), b(1)) > 0) And (iedge(AA(i), b(2)) > 0) And (iedge(AA(i), b(3)) > 0) Then A(2) = AA(i)
    If (iedge(AA(i), b(2)) > 0) And (iedge(AA(i), b(3)) > 0) And (iedge(AA(i), b(4)) > 0) Then A(3) = AA(i)
    If (iedge(AA(i), b(1)) > 0) And (iedge(AA(i), b(3)) > 0) And (iedge(AA(i), b(4)) > 0) Then A(4) = AA(i)
  Next i

End Sub 'ReorderAB



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Proceed4Coord(iedge() As Integer, dihedr() As Double, atnum() As Integer, x() As Double, y() As Double, z() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Treat 4-coordination shapes                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i, j As Integer
Dim k1, l1 As Integer           ' edges of dt1
Dim k2, l2 As Integer           ' edges of dt2
Dim MinDt As Double
Dim Dev(2) As Double            ' squared deviation from ideal polyhedra
Dim dt(2) As Double             ' delta angles
Dim StdAng(2, 2) As Double      ' standard delta1-2 for tetrahedron and tetragon
Dim Poly(2) As String           ' possible polyhedron
Dim DeltaBox As Object          ' TextBox containing delta angles
Dim DeltaBoxCont As String      ' TextBox contents with delta angles
Dim PolyhBox As Object          ' TextBox containing polyhedron information
Dim PolyhBoxCont As String      ' TextBox contents with polyhedron information
Dim AllDeltaBox As Object       ' TextBox containing delta angles for all polyhedra
Dim AllDeltaBoxCont As String   ' TextBox contents with delta angles for all polyhedra

  ' Standard shapes
  StdAng(1, 1) = 109.467
  StdAng(1, 2) = 109.467
  Poly(1) = "Td-TETRAHEDRON"
  StdAng(2, 1) = 0
  StdAng(2, 2) = 0
  Poly(2) = "D4h-TETRAGON"

  MinDt = 2 * Pi * 180 / Pi
  For i = 1 To nvert
    For j = 1 To nvert
      If (iedge(i, j) > 0) Then
        If (dihedr(iedge(i, j)) < MinDt) Then
          MinDt = dihedr(iedge(i, j))
          k1 = i
          l1 = j
        End If
      End If
    Next j
  Next i

  ' Find the second dihedral
  For i = 1 To nvert
    If (i <> k1) And (i <> l1) Then
      k2 = i
      Exit For
    End If
  Next i
  For i = 1 To nvert
    If (i <> k1) And (i <> l1) And (i <> k2) Then
      l2 = i
      Exit For
    End If
  Next i

  If (MinDt = zero) Then
    iedge(k2, l2) = 6
    dihedr(6) = 0
  End If

  dt(1) = dihedr(iedge(k1, l1))
  dt(2) = dihedr(iedge(k2, l2))
  DeltaBoxCont = "Dihedral angles delta which characterize polyhedron:"
  DeltaBoxCont = DeltaBoxCont & CR & "delta1 (belonging to the edge " & Str(atnum(k1)) & "-" & Str(atnum(l1)) & ") = " & FStr(dt(1), 2, 1) & " deg."
  DeltaBoxCont = DeltaBoxCont & CR & "delta2 (belonging to the edge " & Str(atnum(k2)) & "-" & Str(atnum(l2)) & ") = " & FStr(dt(2), 2, 1) & " deg."
  DeltaBox = rspage.TextBoxes.AddEmpty
  DeltaBox.SetContent (DeltaBoxCont)
  DeltaBox.SetBound(1200,700,1000,150)

  AllDeltaBoxCont = "Delta angles in the ideal polyhedra:" & CR
  For i = 1 To 2
    AllDeltaBoxCont = AllDeltaBoxCont & Poly(i) & ":" & CR
      For j = 1 To 2
        AllDeltaBoxCont = AllDeltaBoxCont & FStr(StdAng(i, j), 2, 1)
          If (j < 2) Then AllDeltaBoxCont = AllDeltaBoxCont & ", " Else AllDeltaBoxCont = AllDeltaBoxCont & " deg." & CR
      Next j
  Next i
  AllDeltaBox = rspage.TextBoxes.AddEmpty
  AllDeltaBox.SetContent (AllDeltaBoxCont)
  AllDeltaBox.SetBound(1200,1400,1000,350)

  For i = 1 To 2
    Dev(i) = 0
    For j = 1 To 2
      Dev(i) = Dev(i) + Square(dt(j) - StdAng(i, j)) / (2 - 1)
    Next j
    Print Dev(i)
  Next i
  PolyhBoxCont = "The proximate polyhedron is" & CR
  If (Dev(1) <= Dev(2)) Then PolyhBoxCont = PolyhBoxCont & Poly(1) Else PolyhBoxCont = PolyhBoxCont & Poly(2)
  PolyhBox = rspage.TextBoxes.AddEmpty
  PolyhBox.SetContent (PolyhBoxCont)
  PolyhBox.SetBound(1600,500,700,100)

End Sub 'Proceed4Coord



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Proceed6Coord(iedge() As Integer, dihedr() As Double, atnum() As Integer, x() As Double, y() As Double, z() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Treat 6-coordination shapes                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i, j, k, l As Integer       ' ONLY BY B1-criterion
Dim b(2, 3) As Integer          ' 3 pairs of atoms forming edges b1
Dim MinDt As Double
Dim IsMin As Boolean
Dim count(nvert) As Integer     ' number of eges at each vertex
Dim IsOctFam As Boolean         ' does polyhedron belongs to the octahedron family
Dim IsInReacPth As Boolean      ' does polyhedron shape lie on the reaction path octahedron-trigonal prism
Dim Dev(2) As Double            ' squared deviation from ideal polyhedra
Dim dt1(3) As Double            ' delta angles
Dim StdAng(2, 3) As Double      ' standart delta1-2 for trig. prism and octahedron
Dim Poly(2) As String           ' possible polyhedron
Dim DeltaBox As Object          ' TextBox containing delta angles
Dim DeltaBoxCont As String      ' TextBox contents with delta angles
Dim PolyhBox As Object          ' TextBox containing polyhedron information
Dim PolyhBoxCont As String      ' TextBox contents with polyhedron information
Dim AllDeltaBox As Object       ' TextBox containing delta angles for all polyhedra
Dim AllDeltaBoxCont As String   ' TextBox contents with delta angles for all polyhedra
Dim foundb1 As Integer

  ' Standard shapes
  StdAng(1, 1) = 0
  StdAng(1, 2) = 0
  StdAng(1, 3) = 0
  Poly(1) = "D3h-TRIGONAL PRISM"
  StdAng(2, 1) = 70.5
  StdAng(2, 2) = 70.5
  StdAng(2, 3) = 70.5
  Poly(2) = "D3d-OCTAHEDRON"


  IsOctFam = TRUE
  For i = 1 To nvert
    count(i) = 0
    For j = 1 To nvert
      If (iedge(i, j) > 0) Then
        count(i) = count(i) + 1
      End If
    Next j
    If (count(i) <> 4) Then IsOctFam = FALSE
  Next i

  ' For octahedron family
  If (IsOctFam) Then
    ' Find b1-dihedrals
    foundb1=0
    For k = 1 To 3
      MinDt = 2 * Pi * 180 / Pi
      For i = 1 To nvert
        For j = 1 To nvert
        If i<>j Then
          IsMin = True
            If (iedge(i, j) > 0) Then
              If (dihedr(iedge(i, j)) <= MinDt) Then
                If (k > 1) Then
                  For l = 1 To k - 1
                    If (b(1, l) = i) Or (b(2, l) = j) Or (b(2, l) = i) Or (b(1, l) = j) Then IsMin = False
                  Next l
                End If 'k
                If (IsMin) Then
                  MinDt = dihedr(iedge(i, j))
                  b(1, k) = i
                  b(2, k) = j
                  foundb1=k
                End If 'IsMin
              End If 'dihedr
            End If 'iedge
          End If
          Next j
        Next i
      dt1(k) = MinDt 'dihedr(iedge(b1(1, k), b1(2, k)))
    Next k
    IsInReacPth = True
    For i = 1 To nvert
      For j = 1 To nvert
        If (iedge(i, j) > 0) Then
          If (dihedr(iedge(i, j)) < dt1(3)) And (iedge(i, j) <> iedge(b(1, 1), b(2, 1))) And (iedge(i, j) <> iedge(b(1, 2), b(2, 2))) Then IsInReacPth = False
        End If
      Next j
    Next i

    DeltaBoxCont = "Dihedral angles delta belonging to edges of the type b1:"
    If foundb1>=3 Then
    For i = 1 To 3
      DeltaBoxCont = DeltaBoxCont & CR & "delta" & Str(i) & " (belonging to the edge " & Str(atnum(b(1, i))) & "-" & Str(atnum(b(2, i))) & ") = " & FStr(dt1(i), 2, 1) & " deg."
    Next i
    Else
      DeltaBoxCont = DeltaBoxCont & CR & " * UNDEFINED * Could not find three b1-type dihedrals"
    End If
    If (IsInReacPth = False) Then DeltaBoxCont = DeltaBoxCont & CR & CR & "Warning! The shape of the polyhedron" & CR & " deviates from the reaction pathway" & CR & Poly(2) & " <-> " & CR & Poly(1)
    Set DeltaBox = rspage.TextBoxes.AddEmpty
    DeltaBox.SetContent (DeltaBoxCont)
    DeltaBox.SetBound(1200,700,1100,450)

    AllDeltaBoxCont = "Delta angles in the ideal polyhedra:" & CR
    For i = 1 To 2
      AllDeltaBoxCont = AllDeltaBoxCont & Poly(i) & ":" & CR
      For j = 1 To 3
        AllDeltaBoxCont = AllDeltaBoxCont & FStr(StdAng(i, j), 2, 1)
        If (j < 3) Then AllDeltaBoxCont = AllDeltaBoxCont & ", " Else AllDeltaBoxCont = AllDeltaBoxCont & " deg." & CR
      Next j
    Next i
    AllDeltaBox = rspage.TextBoxes.AddEmpty
    AllDeltaBox.SetContent (AllDeltaBoxCont)
    AllDeltaBox.SetBound(1200,1400,1000,350)

    For i = 1 To 2
      Dev(i) = 0
      For j = 1 To 3
        Dev(i) = Dev(i) + Square(dt1(j) - StdAng(i, j)) / (3 - 1)
      Next j
      Print Dev(i)
    Next i
    PolyhBoxCont = "The proximate polyhedron is" & CR
    If (Dev(1) <= Dev(2)) Then PolyhBoxCont = PolyhBoxCont & Poly(1) Else PolyhBoxCont = PolyhBoxCont & Poly(2)


  Else    'If IsOctFam


    PolyhBoxCont = "This polyhedron does not belong" & CR & "to the octahedron family"


  End If  'If IsOctFam


  PolyhBox = rspage.TextBoxes.AddEmpty
  PolyhBox.SetContent (PolyhBoxCont)
  PolyhBox.SetBound(1600,500,700,200)
End Sub 'Proceed6Coord



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Proceed5Coord(iedge() As Integer, dihedr() As Double, atnum() As Integer, x() As Double, y() As Double, z() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Treat 5-coordination shapes                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i, j As Integer
Dim E(3) As Integer             ' 3 E vertices
Dim nE As Integer               ' number of E vertices
Dim Dev(2) As Double            ' squared deviation from ideal polyhedra
Dim dt(3) As Double             ' delta angles
Dim b(2, 3) As Integer          ' 3 pairs of atoms formative edges e
Dim StdAng(2, 3) As Double      ' standart delta1-2 for trig. prism and octahedron
Dim IsTrBypFam As Boolean       ' does polyhedron belongs to the trig. bipyr. family
Dim count(nvert) As Integer     ' number of eges at each vertex
Dim Poly(2) As String           ' possible polyhedron
Dim DeltaBox As Object          ' TextBox containing delta angles
Dim DeltaBoxCont As String      ' TextBox contents with delta angles
Dim PolyhBox As Object          ' TextBox containing polyhedron information
Dim PolyhBoxCont As String      ' TextBox contents with polyhedron information
Dim AllDeltaBox As Object       ' TextBox containing delta angles for all polyhedra
Dim AllDeltaBoxCont As String   ' TextBox contents with delta angles for all polyhedra

  ' Standard shapes
  StdAng(1, 1) = 53.1
  StdAng(1, 2) = 53.1
  StdAng(1, 3) = 53.1
  Poly(1) = "D3h-TRIGONAL BIPYRAMID"
  StdAng(2, 1) = 0
  StdAng(2, 2) = 75.7
  StdAng(2, 3) = 75.7
  Poly(2) = "C4v-TETRAGONAL PYRAMID"

  IsTrBypFam = TRUE
  nE = 0
  For i = 1 To nvert
    count(i) = 0
    For j = 1 To nvert
      If (iedge(i, j) > 0) Then
        count(i) = count(i) + 1
      End If
    Next j
    If (count(i) = 4) Then
      If (nE < 3) Then
        nE = nE + 1
        E(nE) = i
      Else
        IsTrBypFam = FALSE
      End If
    End If
    If (count(i) > 4) Or (count(i) < 3) Then IsTrBypFam = FALSE
  Next i
  If (nE < 3) Then IsTrBypFam = FALSE

  ' For TBP family
  If (IsTrBypFam) Then

    b(1, 1) = E(1)
    b(2, 1) = E(2)
    b(1, 2) = E(1)
    b(2, 2) = E(3)
    b(1, 3) = E(2)
    b(2, 3) = E(3)

    dt(1) = dihedr(iedge(E(1), E(2)))
    dt(2) = dihedr(iedge(E(1), E(3)))
    dt(3) = dihedr(iedge(E(2), E(3)))

    ' Rearrange dt by increase
    For i = 1 To 2
      For j = i + 1 To 3
        If (dt(i) > dt(j)) Then
          Call swapd(dt(i), dt(j))
          Call swap(b(1, i), b(1, j))
          Call swap(b(2, i), b(2, j))
        End If
      Next j
    Next i

    DeltaBoxCont = "Dihedral angles delta belonging to edges of the type e:"
    For i = 1 To 3
      DeltaBoxCont = DeltaBoxCont & CR & "delta" & Str(i) & " (belonging to the edge " & Str(atnum(b(1, i))) & "-" & Str(atnum(b(2, i))) & ") = " & FStr(dt(i), 2, 1) & " deg."
    Next i
    DeltaBox = rspage.TextBoxes.AddEmpty
    DeltaBox.SetContent (DeltaBoxCont)
    DeltaBox.SetBound(1200,700,1100,200)

    AllDeltaBoxCont = "Delta angles in the ideal polyhedra:" & CR
    For i = 1 To 2
      AllDeltaBoxCont = AllDeltaBoxCont & Poly(i) & ":" & CR
      For j = 1 To 3
        AllDeltaBoxCont = AllDeltaBoxCont & FStr(StdAng(i, j), 2, 1)
        If (j < 3) Then AllDeltaBoxCont = AllDeltaBoxCont & ", " Else AllDeltaBoxCont = AllDeltaBoxCont & " deg." & CR
      Next j
    Next i
    AllDeltaBox = rspage.TextBoxes.AddEmpty
    AllDeltaBox.SetContent (AllDeltaBoxCont)
    AllDeltaBox.SetBound(1200,1400,1000,350)

    For i = 1 To 2
      Dev(i) = 0
      For j = 1 To 3
        Dev(i) = Dev(i) + Square(dt(j) - StdAng(i, j)) / (3 - 1)
      Next j
      Print Dev(i)
    Next i
    PolyhBoxCont = "The proximate polyhedron is" & CR
    If (Dev(1) <= Dev(2)) Then PolyhBoxCont = PolyhBoxCont & Poly(1) Else PolyhBoxCont = PolyhBoxCont & Poly(2)


 Else   'IsTrBypFam


    PolyhBoxCont = "This polyhedron does not belong" & CR & "to the trigonal bipyramid family"

 End If 'IsTrBypFam



 Set PolyhBox = rspage.TextBoxes.AddEmpty
 PolyhBox.SetContent (PolyhBoxCont)
 PolyhBox.SetBound(1600,500,700,200)

End Sub 'Proceed5Coord



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Proceed7Coord(iedge() As Integer, dihedr() As Double, atnum() As Integer, x() As Double, y() As Double, z() As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Treat 7-coordination shapes                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i, j As Integer
Dim E(3) As Integer             ' 5-edg vertices
Dim nE As Integer               ' number of 5-edg vertices
Dim Dev(3) As Double            ' squared deviation from ideal polyhedra
Dim dt(3) As Double             ' delta angles
Dim b(2, 3) As Integer          ' 3 pairs of atoms formative edges u
Dim MinDt As Double
Dim IdPol As Integer
Dim StdAng(3, 4) As Double      ' standard delta1-2 for trig. prism and octahedron and number of
Dim IsCpOctFam As Boolean       ' does polyhedron belongs to the capped octahedron family
Dim IsInReacPth As Boolean      ' does polyhedron shape lie on the reaction path C3v-C2v-D5h
Dim count(nvert) As Integer     ' number of eges at each vertex
Dim Poly(3) As String           ' possible polyhedron
Dim DeltaBox As Object          ' TextBox containing delta angles
Dim DeltaBoxCont As String      ' TextBox contents with delta angles
Dim PolyhBox As Object          ' TextBox containing polyhedron information
Dim PolyhBoxCont As String      ' TextBox contents with polyhedron information
Dim AllDeltaBox As Object       ' TextBox containing delta angles for all polyhedra
Dim AllDeltaBoxCont As String   ' TextBox contents with delta angles for all polyhedra

  ' Standard shapes
  StdAng(1, 1) = 24.2
  StdAng(1, 2) = 24.2
  StdAng(1, 3) = 24.2
  StdAng(1, 4) = 3
  Poly(1) = "C3v-CAPPED OCTAHEDRON"
  StdAng(2, 1) = 0
  StdAng(2, 2) = 0
  StdAng(2, 3) = 41.5
  StdAng(2, 4) = 3
  Poly(2) = "C2v-MONOCAPPED TRIGONAL PRISM"
  StdAng(3, 1) = 54.4
  StdAng(3, 2) = 54.4
  ' StdAng(3, 3) = 0 ' does not exist !!!
  StdAng(3, 4) = 2
  Poly(3) = "D5h-PENTAGONAL BIPYRAMID"

  IsInReacPth = TRUE
  IsCpOctFam = TRUE
  nE = 0
  For i = 1 To nvert
    count(i) = 0
    For j = 1 To nvert
      If (iedge(i, j) > 0) Then
        count(i) = count(i) + 1
      End If
    Next j
    If (count(i) = 5) Then
      If (nE < 3) Then
        nE = nE + 1
        E(nE) = i
      Else
        IsCpOctFam = False
      End If
    End If
    If (count(i) > 5) Or (count(i) < 3) Then IsCpOctFam = False
  Next i
  If (nE < 2) Then IsCpOctFam = FALSE 'maybe needs additional control if belongs to the family


  ' For capped octahedron family
  If (IsCpOctFam) Then
    If (nE = 3) Then ' C3v or C2v
      b(1, 1) = E(1)
      b(2, 1) = E(2)
      b(1, 2) = E(1)
      b(2, 2) = E(3)
      b(1, 3) = E(2)
      b(2, 3) = E(3)
      dt(1) = dihedr(iedge(E(1), E(2)))
      dt(2) = dihedr(iedge(E(1), E(3)))
      dt(3) = dihedr(iedge(E(2), E(3)))
    Else ' D5h
      ' find dt1
      MinDt = 2 * Pi * 180 / Pi
      For i = 1 To 2
        For j = 1 To nvert
          If (iedge(E(i), j) > 0) Then
            If (dihedr(iedge(E(i), j)) < MinDt) Then
              MinDt = dihedr(iedge(E(i), j))
              b(2, 1) = j
              b(1, 1) = E(i)
            End If
          End If
        Next j
      Next i
      dt(1) = MinDt
      ' Find dt2
      b(1, 2) = b(1, 1)
      MinDt = 2 * Pi * 180 / Pi
      For j = 1 To nvert
        If (iedge(E(1), j) > 0) Then
          If (dihedr(iedge(b(1, 2), j)) < MinDt) And (iedge(b(2, 1), j) = 0) And (b(2, 1) <> j) Then
            MinDt = dihedr(iedge(b(1, 2), j))
            b(2, 2) = j
          End If
        End If
      Next j
      dt(2) = MinDt
    End If 'nE

    ' Rearrange dt by increase
    For i = 1 To nE - 1
      For j = i + 1 To nE
        If (dt(i) > dt(j)) Then
          Call swapd(dt(i), dt(j))
          Call swap(b(1, i), b(1, j))
          Call swap(b(2, i), b(2, j))
        End If
      Next j
    Next i

    DeltaBoxCont = "Dihedral angles delta belonging to edges of the type u:"
    For i = 1 To nE
      DeltaBoxCont = DeltaBoxCont & CR & "delta" & Str(i) & " (belonging to the edge " & Str(atnum(b(1, i))) & "-" & Str(atnum(b(2, i))) & ") = " & FStr(dt(i), 2, 1) & " deg."
    Next i
    DeltaBox = rspage.TextBoxes.AddEmpty
    AllDeltaBoxCont = "Delta angles in the ideal polyhedra:" & CR
    For i = 1 To 3
      AllDeltaBoxCont = AllDeltaBoxCont & Poly(i) & ":" & CR
      For j = 1 To Fix(StdAng(i, 4))
        AllDeltaBoxCont = AllDeltaBoxCont & FStr(StdAng(i, j), 2, 1)
        If (j < Fix(StdAng(i, 4))) Then AllDeltaBoxCont = AllDeltaBoxCont & ", " Else AllDeltaBoxCont = AllDeltaBoxCont & " deg." & CR
      Next j
    Next i
    AllDeltaBox = rspage.TextBoxes.AddEmpty
    AllDeltaBox.SetContent (AllDeltaBoxCont)
    AllDeltaBox.SetBound(1200,1400,1000,350)
    For i = 1 To 3
      Dev(i) = 0
      For j = 1 To minii(nE, Fix(StdAng(i, 4)))
        Dev(i) = Dev(i) + Square(dt(j) - StdAng(i, j)) / (minii(nE, Fix(StdAng(i, 4))) - 1)
      Next j
      Print Dev(i)
    Next i
    PolyhBoxCont = "The proximate polyhedron is" & CR
    If (Dev(1) <= Dev(2)) And (Dev(1) <= Dev(3)) Then IdPol = 1
    If (Dev(2) <= Dev(1)) And (Dev(2) <= Dev(3)) Then IdPol = 2
    If (Dev(3) <= Dev(1)) And (Dev(3) <= Dev(2)) Then IdPol = 3
    PolyhBoxCont = PolyhBoxCont & Poly(IdPol)
    If (IdPol < 3) And (nE < 3) Then IsInReacPth = FALSE
    If (IsInReacPth = False) Then
      DeltaBoxCont = DeltaBoxCont & CR & CR & "Warning! The shape of the polyhedron" & CR & " deviates from the reaction pathway " & CR & Poly(1) & " <-> " & CR & Poly(2) & " <-> " & CR & Poly(3)
      DeltaBoxCont = DeltaBoxCont & CR & "Check your polyhedron for the" & CR & Poly(2)
    End If
    DeltaBox.SetContent (DeltaBoxCont)
    DeltaBox.SetBound(1200,700,1100,550)

  Else

    PolyhBoxCont = "This polyhedron does not belong" & CR & "to the capped octahedron family"

  End If 'IsCpOctFam


  PolyhBox = rspage.TextBoxes.AddEmpty
  PolyhBox.SetContent (PolyhBoxCont)
  PolyhBox.SetBound(1600,500,700,100)
End Sub 'Proceed7Coord



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function GetPhi(A1 As Integer, A2 As Integer, B1 As Integer, B2 As Integer, x() As Double, y() As Double, z() As Double) As Double
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim x3, y3, z3, Ap1, Bp1, Cp1, Ap2, Bp2, Cp2, D As Double
 x3 = (x(B1) + x(B2)) / 2
 y3 = (y(B1) + y(B2)) / 2
 z3 = (z(B1) + z(B2)) / 2
 Call GetPlane(Ap1, Bp1, Cp1, D, x(A1), x(A2), x3, y(A1), y(A2), y3, z(A1), z(A2), z3)
 x3 = (x(A1) + x(A2)) / 2
 y3 = (y(A1) + y(A2)) / 2
 z3 = (z(A1) + z(A2)) / 2
 Call GetPlane(Ap2, Bp2, Cp2, D, x(B1), x(B2), x3, y(B1), y(B2), y3, z(B1), z(B2), z3)
 GetPhi = Pi - GetAngle(Ap1, Bp1, Cp1, Ap2, Bp2, Cp2)
End Function 'GetPhi



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub GetPlane(A As Double, b As Double, C As Double, D As Double, x1 As Double, x2 As Double, x3 As Double, y1 As Double, y2 As Double, y3 As Double, z1 As Double, z2 As Double, z3 As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 A = (y2 - y1) * (z3 - z1) - (z2 - z1) * (y3 - y1)
 b = (x3 - x1) * (z2 - z1) - (z3 - z1) * (x2 - x1)
 C = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)
 D = x1 * (-(y2 - y1) * (z3 - z1) + (z2 - z1) * (y3 - y1)) + y1 * (-(x3 - x1) * (z2 - z1) + (z3 - z1) * (x2 - x1)) + z1 * (-(x2 - x1) * (y3 - y1) + (y2 - y1) * (x3 - x1))
End Sub 'GetPlane



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function GetAngle(A1 As Double, B1 As Double, C1 As Double, A2 As Double, B2 As Double, C2 As Double) As Double
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim cosine As Double
Dim IsGot As Boolean
 IsGot = False
 cosine = (A1 * A2 + B1 * B2 + C1 * C2) / Sqrt((A1 * A1 + B1 * B1 + C1 * C1) * (A2 * A2 + B2 * B2 + C2 * C2))
 If (cosine >= one) Then
  IsGot = True
  GetAngle = Pi - zero
 End If
 If (cosine <= -one) Then
  IsGot = True
  GetAngle = Pi - Pi
 End If
  If (IsGot = False) Then GetAngle = Pi - (Atan(-cosine / Sqrt(-cosine * cosine + 1)) + 2 * Atan(1)) ' arccos
End Function 'GetAngle



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub GetIntersPoint(xi As Double, yi As Double, zi As Double, A As Double, b As Double, C As Double, D As Double, x() As Double, y() As Double, z() As Double, p1 As Integer, p2 As Integer)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim l,m,n,p,q As Double
  l = x(p2) - x(p1)
  m = y(p2) - y(p1)
  n = z(p2) - z(p1)
  q = A * l + b * m + C * n
  ' If q equals to zero (or near to) then intersection point
  ' is in infinity (signum does not matter)
  If (Abs(q) < DOUBLE_ZERO) Then p = 1 / DOUBLE_ZERO Else p = (A * x(p1) + b * y(p1) + C * z(p1) + D) / q
  xi = x(p1) - l * p
  yi = y(p1) - m * p
  zi = z(p1) - n * p
End Sub 'GetIntersPoint




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function mini(val1 As Double, val2 As Double) As Double
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  If (val1 > val2) Then mini = val2 Else mini = val1
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function minii(ByVal val1 As Integer, ByVal val2 As Integer) As Integer
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  If (val1 > val2) Then minii = val2 Else minii = val1
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function maxi(val1 As Double, val2 As Double) As Double
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  If (val1 < val2) Then maxi = val2 Else maxi = val1
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Square(ByVal val1 As Double) As Double
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  Square = val1 * val1
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub swap(val1 As Integer, val2 As Integer)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim tmp As Integer
  tmp = val1
  val1 = val2
  val2 = tmp
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub swapd(val1 As Double, val2 As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim tmp As Double
  tmp = val1
  val1 = val2
  val2 = tmp
End Sub



'***LIBRARY PROCEDURES BEGIN



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub RefreshDiagram(diag As Object,strmol As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' LIB0.BAS PROCEDURE                                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Redraws the diagram with a molecule or structure object             '
'                                                                     '
' ENTER                                                               '
'     diag            object of type CS_DIAGRAM                       '
'     strmol          object of type CB_MOLECULE or CB_STRUCTURE      '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim l,t,w,h,w1,h1 As Integer
  diag.GetBound(l,t,w,h)
  diag.Depict(strmol)
  diag.GetBound(w,h,w1,h1)
  diag.SetBound(l,t,w1,h1)
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function LabelDiagramWithNumbers(diag As Object) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' LIB0.BAS PROCEDURE                                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Redraws the diagram showing order numbers at atoms                  '
' EXIT                                                                '
'     returns TRUE at success otherwise FALSE                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim asm,struc,at As Object, i,nat As Integer
  LabelDiagramWithNumbers=FALSE
  asm=Assemblies.AddFromCS(diag)
  If asm=NULL Then Exit Function
  struc=asm.Structures.Item(1)
  If struc=NULL Then Exit Function
  ' Supply atomic labels
  With asm
    nat=.Count
    For i=1 To nat
      at=.Item(i)
      at.SetName(Str(i))
    Next i
  End With
  ' Show labelled diagram
  RefreshDiagram(diag,struc)
  LabelDiagramWithNumbers=TRUE
End Function



'***LIBRARY PROCEDURES END

'@@@@@@