'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' ACD CHEMBASIC DEMO PROGRAM                                          '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'                                                                     '
' Molecular 3D Editor//PUCKER.BAS                                     '
'                                                                     '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'                                                                     '
' The utility measures Cremer-Pople puckering parameters for the ring '
'                                                                     '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''


CONST TITLE="ChemBasic Molecular Editor // Ring puckering"
CONST PI=3.1415926
CONST TWOPI=6.2831852
CONST LTB=200           ' textboxes layout
CONST TTB=300
CONST HTB=120

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Main As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' PUCKER.BAS                                                          '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim page,diag,asm,struc As Object, natring,iatring() As Integer
Dim squery As String, OK As Boolean

  MAIN="Failed or nothing to do!"

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

  ' Do the job
  OK=LabelDiagramWithNumbers(diag)
  If Not OK Then Exit Function
  squery=UCase(UserIOBox("Please supply ring atoms" ,TITLE ,  ""))
  OK=ParseQueryRing(struc,squery,natring,iatring)
  If Not OK Then Exit Function
  OK=CremerPoplePuckering(diag,struc,natring,iatring)
  Call ClearDiagramLabels(diag)
  Main="Completed."


End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function ParseQueryRing(struc As object,ByVal s As String,nr As Integer, iatring() As Integer) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Extract (and check for validity) ring atom numbers from a string    '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim nat,k As Integer, ss(25) As String, OK As Boolean, asm As Object

  ParseQueryRing=False

  If s="" Then Exit Function

  ' Extract the numbers
  nr=SubStrings(Trim(s)," :-;,",ss)
  ReDim iatring(nr) As Integer

  If (nr>=5) And (nr<=25) Then
    For k=1 to nr
      iatring(k)=Fix(Val(ss(k)))
    Next k
  Else
    If nr>25 Then
      MessageBox("Could not analyze such a large ring", TITLE, MBB_OK + MBI_EXCLAMATION)
    Else
      MessageBox("Could not extract [more than 4] atoms", TITLE, MBB_OK + MBI_EXCLAMATION)
    End If
    Exit Function
  End If

  ' Check them for validity
  asm=struc.Assembly
  nat=asm.Count
  For k=1 to nr
    If (iatring(k)<0) Or (iatring(k)>nat) Then
      MessageBox("Bad atom numbers sequence ("+s+")", TITLE, MBB_OK + MBI_EXCLAMATION)
      Exit Function
    End If
    If Not struc.IsRing(asm.Item(iatring(k))) Then
      MessageBox("Atom "+Str(iatring(k))+" is outside ring!", TITLE, MBB_OK + MBI_EXCLAMATION)
      Exit Function
    End If

  Next k
  ' here have not checked atoms for equal numbers, sorry....

  ParseQueryRing=True

End Function




''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function CremerPoplePuckering(diag As Object,struc As Object,natring As Integer, iatring() As Integer) As Boolean
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Follows the definition given by Cremer & Pople in JACS,v.97,p.1354(1975)'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,j,m,pivot,npar,k,FF,l,t,w,h As Integer
Dim vect(3),x,y,z,x1,y1,z1,x2,y2,z2,x3,y3,z3,nvect(3),R1(3),R2(3),sm,cm,argp As Double
Dim dcos,dsin,delta,qcos(),qsin(),q(),psi(),qlast,qtot As Double
Dim squery As String
Dim diag1,ring,atring(1),pout,tbox As Object

  CremerPoplePuckering=FALSE

  ' Collect the atoms
  ring=struc.Assembly.CopyTree(True).Structures.Item(1)
  ReDim atring(natring) As Object
  For i=1 to natring
    j=iatring(i)
    atring(i)=ring.Assembly.Item(j)
  Next i

  'Move origin(zero-point) to the geom. center
  AdjustXYZToCenter(ring)

  ' Define R1, R2 and N
  R1(1)=0 :R1(2)=0 : R1(3)=0
  R2(1)=0 :R2(2)=0 : R2(3)=0
  For i=1 To natring
    ring.GetAtomXYZ(atring(i),x,y,z)
    sm=sin(TWOPI*Dbl(i-1)/Dbl(natring))
    cm=cos(TWOPI*Dbl(i-1)/Dbl(natring))
    R1(1)=R1(1)+x*sm
    R1(2)=R1(2)+y*sm
    R1(3)=R1(3)+z*sm
    R2(1)=R2(1)+x*cm
    R2(2)=R2(2)+y*cm
    R2(3)=R2(3)+z*cm
  Next i
  nvect(1)=R1(2)*R2(3)-R2(2)*R1(3)
  nvect(2)=-(R1(1)*R2(3)-R2(1)*R1(3))
  nvect(3)=R1(1)*R2(2)-R1(2)*R2(1)


  'Project Z along nvect
  x1=nvect(1)
  y1=nvect(2)
  z1=Abs(nvect(3))
  x2=0
  y2=0
  z2=SQRT(nvect(1)*nvect(1)+nvect(2)*nvect(2)+nvect(3)*nvect(3))
  vect(1)=(x1+x2)/2.0
  vect(2)=(y1+y2)/2.0
  vect(3)=(z1+z2)/2.0
  For i=1 to natring
    RotateAtomPi(ring,atring(i),vect)
  Next i

  'Project Y along projection of user defined atom on main plane
  pivot=Fix(Val(UserIOBox("Please supply a pivot ring atoms",TITLE,Str(iatring(1)))))
  ring.GetAtomXYZ(ring.Assembly.Item(pivot),x1,y1,z1)
  z1=0
  x2=0
  y2=SQRT(x1*x1+y1*y1)
  z2=0
  vect(1)=(x1+x2)/2.0
  vect(2)=(y1+y2)/2.0
  vect(3)=(z1+z2)/2.0
  For i=1 To natring
    RotateAtomPi(ring,atring(i),vect)
  Next i
  vect(1)=0.0
  vect(2)=1.0
  vect(3)=0.0
  For i=1 To natring
    RotateAtomPi(ring,atring(i),vect)
  Next i

  'Define a positive direction of Z axis
  If DefPositivZ(ring,atring,natring)>0 Then
    vect(1)=1.0
    vect(2)=0.0
    vect(3)=0.0
    For i=1 to natring
      RotateAtomPi(ring,atring(i),vect)
    Next i
  End If

  pout=ActiveDocument.AddEmpty

  ' Draw a copy structure
  diag1=diag.LoadOnTo(pout)
  diag1.GetBound(l,t,w,h)
  diag1.SetBound(1350,TTB,w,h)

  'Output final cartesian coordinates
  tbox=pout.TextBoxes.AddEmpty
  tbox.SetContent("CREMER AND POPLE PUCKERING CORDINATES")
  tbox.SetBound(LTB,TTB,1000,50)
  tbox=pout.TextBoxes.AddEmpty
  tbox.SetContent("FINAL CARTESIAN COORDINATES FOR THE RING ARE:")
  tbox.SetBound(LTB,TTB+HTB,1000,50)
  k=TTB+HTB*2
  For i=1 to natring
    ring.GetAtomXYZ(atring(i),x,y,z)
    tbox=pout.TextBoxes.AddEmpty
    tbox.SetContent(Str(i)+"  "+atring(i).ElSymbol+"  "+FStr(x,12,5)+"   "+FStr(y,12,5)+"   "+FStr(z,12,5))
    tbox.SetBound(LTB,k,1500,50)
    k=k+HTB\2
  Next i

  'Calculate all the puckering parameters
  k=k+100
  npar=0
  If (natring Mod 2)=0 Then
    npar=Int(natring/2-1)
  Else
    npar=Int((natring-1)/2)
  End If

  ReDim qcos(npar)
  ReDim qsin(npar)
  ReDim psi(npar)
  ReDim q(npar)

  For m=2 to npar
    ' Compute
    qcos(m)=0
    qsin(m)=0
    For j=1 to natring
      argp=TWOPI*Dbl(m)*Dbl((j-1))/Dbl(natring)
      ring.GetAtomXYZ(atring(j),x1,y1,z1)
      dcos= SQRT(2/natring)*z1*cos(argp)
      dsin=-SQRT(2/natring)*z1*sin(argp)
      qcos(m)=qcos(m)+dcos
      qsin(m)=qsin(m)+dsin
    Next j
    psi(m)=atan(qsin(m)/qcos(m))
    q(m)=qcos(m)/cos(psi(m))
    ' Print
    tbox=pout.TextBoxes.AddEmpty
    tbox.SetContent("psi("+Str(m)+")="+FStr(psi(m),12,5))
    tbox.SetBound(LTB,k,1500,50)
    tbox=pout.TextBoxes.AddEmpty
    tbox.SetContent("q("+Str(m)+")="+FStr(q(m),12,5))
    tbox.SetBound(LTB,k+HTB\2,1500,50)
    k=k+HTB
  Next m

  '
  if (natring Mod 2)=0 then
    qlast=0
    For j=1 to natring
      ring.GetAtomXYZ(atring(j),x1,y1,z1)
      delta=(1/SQRT(natring))*z1*cos((j-1)*PI)
      qlast=qlast+delta
    Next j
    x1=natring/2
    tbox=pout.TextBoxes.AddEmpty
    tbox.SetContent("q("+Str(x1)+")="+FStr(qlast,12,5))
    tbox.SetBound(LTB,k+HTB\2,1500,50)
  End If

  'Calculate the total puckering amplitude
  qtot=0
  For j=1 To natring
    ring.GetAtomXYZ(atring(j),x1,y1,z1)
    qtot=qtot+z1*z1
  Next j
  qtot=sqrt(qtot)
  tbox=pout.TextBoxes.AddEmpty
  tbox.SetContent("Q(total puckering amplitude)="+FStr(qtot,12,5))
  tbox.SetBound(LTB,k+HTB,1500,50)

  CremerPoplePuckering=TRUE

End Function



''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub RotateAtomPi(struc As Object,atm as Object,vector() As Double)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D,Dist,t,calpha,cbeta,cgamma as Double
  struc.GetAtomXYZ(atm,x1,y1,z1)

  a=vector(1)
  b=vector(2)
  c=vector(3)
  d=-(a*x1+b*y1+c*z1)
  t=-d/(a*a+b*b+c*c)
  x2=a*t
  y2=b*t
  z2=c*t
  dist=Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1))
  If dist<=0.0 Then Exit Sub
  calpha=(x2-x1)/Dist
  cbeta=(y2-y1)/Dist
  cgamma=(z2-z1)/Dist
  x3=x1+Dist*(1-cos(3.1415))*calpha
  y3=y1+Dist*(1-cos(3.1415))*cbeta
  z3=z1+Dist*(1-cos(3.1415))*cgamma

  struc.SetAtomXYZ(Atm,x3,y3,z3)
End Sub



''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function DefPositivZ(ring As Object,Arr() As Object,n As Integer) As Integer
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim a1,a2,a3 As Object, a,b,c,d,dd  As Double, i,j,k,Dolli As Integer

  Dolli=0
  For i=1 to n
    j=i+1 : k=i+2
    If i=n-1 Then
     j=n : k=1
    End If
    If i=n Then
     j=1  : k=2
    End If
    a1=Arr(i) : a2=Arr(j) : a3=Arr(k)
    DefPlane(ring,A1,A2,A3,A,B,C,D)
    DD=A*0+B*0+C*100+D
    'Print DD
    If DD>0.0 Then Dolli=Dolli+1
    if DD<0.0 Then Dolli=Dolli-1
  Next i
  'Print Dolli

  DefPositivZ=Dolli
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''




''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub DefPlane(struct as object,Atom1 as object,Atom2 as object,Atom3 as object,A as Double,B as Double,C as Double,D as Double)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim x1,y1,z1,x2,y2,z2,x3,y3,z3,A1,B1,C1,A2,B2,C2 As Double
  'Define coordinates
  struct.GetAtomXYZ(Atom1,x1,y1,z1):z1=1
  struct.GetAtomXYZ(Atom2,x2,y2,z2):z2=1
  struct.GetAtomXYZ(Atom3,x3,y3,z3):z3=1
  'Define plane
  A1=x2-x1:B1=y2-y1:C1=z2-z1
  A2=x3-x1:B2=y3-y1:C2=z3-z1
  A=B1*C2-C1*B2
  B=-(A1*C2-C1*A2)
  C=A1*B2-B1*A2
  D=-(x1*A+y1*B+z1*C)
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



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function ClearDiagramLabels(diag As Object) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' LIB0.BAS PROCEDURE                                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Clears the atomic labels and re-draws diagram                       '
' EXIT                                                                '
'     returns TRUE at success otherwise FALSE                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim asm,struc,at As Object, i,nat As Integer
  ClearDiagramLabels=FALSE
  asm=Assemblies.AddFromCS(diag)
  If asm=NULL Then Exit Function
  struc=asm.Structures.Item(1)
  If struc=NULL Then Exit Function
  ' Clear atomic labels
  For Each at In asm
    at.SetName("")
  Next at
  ' Show delabelled diagram
  RefreshDiagram(diag,struc)
  ClearDiagramLabels=TRUE
End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub AdjustXYZToCenter(strconf As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' LIB0.BAS PROCEDURE                                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Go to center-of-molecule coordinate system                          '
'                                                                     '
' ENTER                                                               '
'     strconf         object of type CB_CONForMATION or CB_STRUCTURE  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim x(1),y(1),z(1),x0,y0,z0,xx,yy,zz As Double, i,natoms As Integer, asm As Object

  If strconf.GetType<>CB_STRUCTURE And strconf.GetType<>CB_CONForMATION Then Exit Sub

  asm=strconf.Assembly
  With asm

    natoms=.Count
    ReDim x(natoms) : ReDim y(natoms) : ReDim z(natoms)
    x0=0.0 : y0=0.0 : z0=0.0
    For i=1 to natoms
      strconf.GetAtomXyz(.Item(i),x(i),y(i),z(i))
      x0=x0+x(i) : y0=y0+y(i) : z0=z0+z(i)
    Next i
    x0=x0/Dbl(natoms) : y0=y0/Dbl(natoms) : z0=z0/Dbl(natoms)
    For i=1 to natoms
      strconf.SetAtomXyz(.Item(i),x(i)-x0,y(i)-y0,z(i)-z0)
    Next i

  End With

End Sub



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function SubStrings(ByVal s As String, ByVal sc As String, ss() As String) As Integer
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' LIB0.BAS PROCEDURE                                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Splits the string onto substrings separated with character          '
' and stores substrings in a string aray                              '
'                                                                     '
' ENTER                                                               '
'     s            source string                                      '
'     sc           separators string (e.g., ".,- " means that         '
'                  '.' ',' '-' and ' ' are possible separators)       '
'                  (CR is always a separator)                         '
' EXIT                                                                '
'     returns number of substrings                                    '
'     ss() is properly re-dimensioned array of sub-strings            '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,le,ns As Integer,  c As string, copy As Boolean
  SubStrings=0
  le=Len(s)
  if le<1 Then Exit Function

  copy=False
  ns=0
  For i=1 To le
    c=Mid(s,i,1)
    'If c=sc Then
    If Instr(1,sc,c)>0 Then
    'separator occurred, toggle copying off or simply skip a character
      If copy Then copy=False
    Else
    'treat normal char
      If Not copy Or i=1 Then
        ns=ns+1
        copy=True
        ss(ns)=""
      End If
      ss(ns)=ss(ns)+c
    End If
  Next i
  SubStrings=ns
End Function


'***LIBRARY PROCEDURES END

'@@@@@@