'
' Peptide Shortcut Expander - 3D Builder
'  
' Revision:     v1.2
' Author:       Robin Martin
' Date:         February 21, 2006
' Description:  Added information for Ornithine        

' Revision:     v1.1
' Author:       Robin Martin
' Date:         March 12, 2004 
' Description:  Added new data forms to ask user to select using text on page or fill in a form.





'Global consts and vars

CONST MaxAA=21                                              'No of reference AA
Dim AAName(MaxAA),AAName1(MaxAA),AAName3(MaxAA) As String   'AA name, 1-letter and 3-letter code
Dim AAFName(MaxAA) As String, AAStructure(MaxAA) As Object  'AA file name and ref. structure
Dim AACalpha(MaxAA),AACend(MaxAA),AANend(MaxAA) As Integer  'AA particular atom numbers
Dim AAOx(MaxAA),AAOhxl(MaxAA),AAH1N(MaxAA),AAH2N(MaxAA) As Integer

' Misc

CONST TITLE="Peptide Builder"                               'ProgName
CONST MaxResidues=1000                                      'MaxSize of peptide
CONST RAD_TO_DEG = 57.29577951                              'Evident?
CONST RIGHT_ALPHA_HELIX=0
CONST BETA_EXTENDED=1
CONST BETA_TURN_ANY=2
Dim NResToShowInfo As Integer

CONST iomd1="Use existing textbox on page"                  'RAM modified
CONST iomd2="Enter code with a form"
CONST fnm1="pepchoose.frm"
CONST fnm2="singleinput.frm"
CONST fnm3="threeinput.frm"
CONST codcho1="SINGLE letter codes"
CONST codcho3="THREE letter codes"

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Main As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim Peptide,PhiAts(4,MaxResidues),PsiAts(4,MaxResidues),diag As Object
Dim IPeptide(MaxResidues),MotifCode(MaxResidues), NResidues, ans,n,i As Integer
Dim OK As Boolean, SPeptide,S1Peptide,mess As String


  MAIN="Failed or nothing to do!"

  'Initialize
  Call Init
  For i=1 To MaxResidues
    MotifCode(i)=RIGHT_ALPHA_HELIX
  Next i

  'Get the shorthand formula with secondary structure labels if any
  S1Peptide=Get1LetterCode(MotifCode)
  If S1Peptide="" Then Exit Function


  'Decrypt shorthand formula
  NResidues=Convert1LetterCodeToNumber(S1Peptide,IPeptide)
  If NResidues<=0 Then Exit Function


  'Get residue structures
  OK=GetResidueStructures(S1Peptide)
  If Not OK Then Exit Function


  'Assign residue atoms
  OK=AssignResidueAtoms
  If Not OK Then Exit Function


  'Prepare to build (ask user for screen updates)
  ans=MessageBox("Ready to build the peptide of"+Chr(13)+Str(NResidues)+" residues",TITLE, MBB_OKCANCEL+MBI_INFORMATION)
  If ans=MBR_CANCEL Then Exit Function
  NResToShowInfo=32000
  If Nresidues>30 Then
    mess="Building may require some time."+Chr(13)+"Report on the progress upon constructing 3D-model?"
    ans=MessageBox(mess,TITLE, MBB_YESNO + MBI_QUESTION)
    If ans=MBR_YES Then
      If NResidues>30 Then NResToShowInfo=5
      If NResidues>50 Then NResToShowInfo=10
      If NResidues>100 Then NResToShowInfo=20
      If NResidues>300 Then NResToShowInfo=50
    End If
  End If


  diag=ActiveDocument.ActivePage.Diagrams.Addempty


  'Build the peptide
  Peptide=BuildThePeptide(IPeptide,NResidues,PhiAts,PsiAts,diag)


  'Shape the peptide
  ans=MessageBox("Building completed successfully."+Chr(13)+"Ready to set a secondary structure",TITLE, MBI_INFORMATION)
  Call SetOverallSecondaryStructure(Peptide,NResidues,PhiAts,PsiAts,diag,MotifCode)


  '3D-Opt?
  ans=MessageBox("Perform 3D-Optimization?",TITLE, MBB_YESNO + MBI_QUESTION)
  If ans=MBR_YES Then
    If Peptide.Assembly.Count<200 Then
      Peptide=Peptide.Do3DOptimize(0.1)
    Else
      MessageBox("The structure is too large to be able to perform optimization! ",TITLE, MBI_INFORMATION)
    End If
  End If


  'Show the peptide
  ShowStruct(diag,Peptide)


  'Export announcement
  If NResidues>30 Then
    mess="The structure is very large so you may wish export it to a MOL file"+Chr(13)+"for use with third party viewing software."
    ans=MessageBox(mess,TITLE,MBI_INFORMATION)
  End If

  'Say good bye
  mess="To improve the viewing of the structure, try removing the explicit hydrogens,"+Chr(13)+"set zoom level to Fit All, and/or copy to 3D-Viewer!"
  MessageBox(mess,TITLE, MBI_INFORMATION)

  MAIN="Completed."

End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Init
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Initialize reference data structures                                        '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  AAName3(1)= "Ala" : AAName1(1)= "A" : AAName(1)= "Alanine"        : AAFName(1)= "AASTRUCT/Ala.mol"
  AAName3(2)= "Glu" : AAName1(2)= "E" : AAName(2)= "Glutamic acid"  : AAFName(2)= "AASTRUCT/Glu.mol"
  AAName3(3)= "Gln" : AAName1(3)= "Q" : AAName(3)= "Glutamine"      : AAFName(3)= "AASTRUCT/Gln.mol"
  AAName3(4)= "Asp" : AAName1(4)= "D" : AAName(4)= "Aspartic acid" : AAFName(4)= "AASTRUCT/Asp.mol"
  AAName3(5)= "Asn" : AAName1(5)= "N" : AAName(5)= "Asparagine"     : AAFName(5)= "AASTRUCT/Asn.mol"
  AAName3(6)= "Leu" : AAName1(6)= "L" : AAName(6)= "Leucine"        : AAFName(6)= "AASTRUCT/Leu.mol"
  AAName3(7)= "Gly" : AAName1(7)= "G" : AAName(7)= "Glycine"        : AAFName(7)= "AASTRUCT/Gly.mol"
  AAName3(8)= "Lys" : AAName1(8)= "K" : AAName(8)= "Lysine"         : AAFName(8)= "AASTRUCT/Lys.mol"
  AAName3(9)= "Ser" : AAName1(9)= "S" : AAName(9)= "Serine"         : AAFName(9)= "AASTRUCT/Ser.mol"
  AAName3(10)="Val" : AAName1(10)="V" : AAName(10)="Valine"         : AAFName(10)="AASTRUCT/Val.mol"
  AAName3(11)="Arg" : AAName1(11)="R" : AAName(11)="Arginine"       : AAFName(11)="AASTRUCT/Arg.mol"
  AAName3(12)="Thr" : AAName1(12)="T" : AAName(12)="Threonine"      : AAFName(12)="AASTRUCT/Thr.mol"
  AAName3(13)="Pro" : AAName1(13)="P" : AAName(13)="Proline"        : AAFName(13)="AASTRUCT/Pro.mol"
  AAName3(14)="Ile" : AAName1(14)="I" : AAName(14)="Isoleucine"     : AAFName(14)="AASTRUCT/Ile.mol"
  AAName3(15)="Met" : AAName1(15)="M" : AAName(15)="Methionine"     : AAFName(15)="AASTRUCT/Met.mol"
  AAName3(16)="Phe" : AAName1(16)="F" : AAName(16)="Phenylalanine"  : AAFName(16)="AASTRUCT/Phe.mol"
  AAName3(17)="Tyr" : AAName1(17)="Y" : AAName(17)="Tyrosine"       : AAFName(17)="AASTRUCT/Tyr.mol"
  AAName3(18)="Cys" : AAName1(18)="C" : AAName(18)="Cysteine"       : AAFName(18)="AASTRUCT/Cys.mol"
  AAName3(19)="Trp" : AAName1(19)="W" : AAName(19)="Tryptophan"    : AAFName(19)="AASTRUCT/Trp.mol"
  AAName3(20)="His" : AAName1(20)="H" : AAName(20)="Histidine"      : AAFName(20)="AASTRUCT/His.mol"
  AAName3(21)="Orn" : AAName1(21)="O" : AAName(21)="Ornithine"      : AAFName(21)="AASTRUCT/Orn.mol"
End Sub



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Get1LetterCode(MotifCode() As Integer) As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Get shorthand formula (and convert it to) 1-letter code                     '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim TBoxes As Object, Text,s,mess As String
Dim cdch,imod as string
Dim Frm1,Frm2 as object



  Get1LetterCode=""
'
'The following section added by Robin Martin on March 12, 2004
'
Frm1=ReadForm(fnm1)
if Frm1.ExecForm then
	cdch=Frm1.GetStrValue("CodeChoice")
	imod=Frm1.GetStrValue("InputMode")
	if imod=iomd2 then                     'Use form input
		if cdch=codcho1 then           'Use Single Letter Code form
			Frm2=ReadForm(fnm2)
			if Frm2.ExecForm then
				Text=Frm2.GetStrValue("SingleCode")
			else
				MessageBox("You pressed Cancel on the form. The program will exit!","Cancel",MBI_EXCLAMATION)
				Exit Function
			end if
			
		else                            'Use Three Letter code form
			Frm2=ReadForm(fnm3)
			if Frm2.ExecForm then
				Text=Frm2.GetStrValue("TripleCode")
			else
				MessageBox("You pressed Cancel on the form. The program will exit!","Cancel",MBI_EXCLAMATION)
				Exit Function
			end if
		end if
		Goto carryon
	end if
else
	MessageBox("You pressed Cancel on the form. The program will exit!","Cancel",MBI_EXCLAMATION)
	Exit Function	
end if ' Frm1
'
' End of modification
'
  TBoxes=ActiveDocument.ActivePage.TextBoxes

  If TBoxes.Count<>1 Then
    mess="A page should include a single text box with peptide shorthand"+Chr(13)
    mess=mess+"(e.g., Thr-Phe-Lys-Thr-Val or TFKTV)."+Chr(13)
    mess=mess+"Please correct and run the program again."
    MessageBox(mess,TITLE, MBB_OK + MBI_STOP)
    Exit Function
  End If

  Text=TBoxes.Item(1).GetContent
carryon:
  If Text="" Then Exit Function

  If (InStr(1,Text,"-")>0) Or (InStr(1,Text,"~")>0) Or (InStr(1,Text,"^")>0)Then
    'Have (probably) 3-letter code; convert it to 1-letter code
    s=Convert3To1LetterCode(Text,MotifCode)
  Else
    'Have (probably) 1-letter code
    s=Text
  End If

  Get1LetterCode=DelNonPrintable(s)

End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function DelNonPrintable(s0 As String) As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Remove CR, etc.                                                             '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim s,c As String, i,l As Integer
  l=Len(s0) : s=""
  For i=1 To l
    c=Mid(s0,i,1) : If c<>Chr(13) And c<>Chr(10)Then s=s+c
  Next i
  DelNonPrintable=s
End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Convert3To1LetterCode(s0 As String,MotifCode() As Integer) As String
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Converts 3- to 1-letter notation                                            '
' Returns 1-letter form or empty string @ fail                                '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim s,s3,sslabel As String, i,j,l,n As Integer

  Convert3To1LetterCode=""
  l=Len(s0)
  s=""
  n=0

  For i=1 To l-2 Step 4
    s3=Mid(s0,i,3)

    For j=1 To MaxAA

      If AAName3(j)=s3 Then

        s=s+AAName1(j)

        'check secndry struct label
        If i<l-2 Then
          n=n+1
          sslabel=Mid(s0,i+3,1)
          If sslabel="~" Then MotifCode(n)=BETA_EXTENDED
          If sslabel="^" Then MotifCode(n)=BETA_TURN_ANY
        End If

        GoTo nexti

      End If

    Next j

    MessageBox("Could not recognize the code: "+Chr(13)+Chr(13)+s3,TITLE, MBB_OK + MBI_STOP)
    Exit Function
  nexti:
  Next i


  Convert3To1LetterCode=s
End Function



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Convert1LetterCodeToNumber(s0 As String,IPeptide() As Integer) As Integer
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Converts 1-letter notation to internal AA type numbers                      '
' Returns N of residues @success or 0  @fail                                  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim s1 As String, i,j,l As Integer

  Convert1LetterCodeToNumber=0

  l=Len(s0)

  For i=1 To l
    s1=Mid(s0,i,1)

    For j=1 To MaxAA
      If AAName1(j)=s1 Then
        IPeptide(i)=j
        GoTo nexti
      End If
    Next j

    MessageBox("Could not recognize the code: "+Chr(13)+Chr(13)+s1,TITLE, MBB_OK + MBI_STOP)
    Exit Function
  nexti:
  Next i

  Convert1LetterCodeToNumber=Len(s0)
End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function GetResidueStructures(S1Peptide As String) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Imports (necessary) structure blocks from file data                         '
' Returns True @success or False  @fail                                       '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim s As String, i,j As Integer, Asm As Object

  GetResidueStructures=False

  For i=1 To MaxAA
    AAStructure(i)=NULL

    If (InStr(1,S1Peptide,AAName1(i))>0) Then
      Asm=Assemblies.AddFromFile(AAFName(i),1)
      If Asm<>NULL Then  AAStructure(i)=Asm.Structures.Item(1)
      If (AAStructure(i)=NULL) Then
        MessageBox("Could not find or treat AA data file <"+AAFName(i)+"> for "+AAName(i),TITLE, MBB_OK + MBI_STOP)
        Exit Function
      End If
    End If

  Next i

  GetResidueStructures=True
End Function





'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function AssignResidueAtoms As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Assigns Calphas, etc. for necessary AA struct. units                        '
' Returns True @success or False  @fail                                       '
'                                                                             '
' TRICK: we save atoms as their integer indices, as we'll further CopyTree the'
'        AA reference structures with (FALSE) argument; this means that atom  '
'        handles will change. However (and the trick is that) CopyTree(FALSE) '
'        will preserve atomic numbering                                       '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i As Integer, struc,asm,mol,at,ca,nend,cend As Object

  AssignResidueAtoms=False

  For i=1 To MaxAA

    struc=AAStructure(i)
    If struc<>NULL Then
      mol=struc.Molecule
      asm=struc.Assembly
      ca=AssignCAlpha(mol,nend,cend)
      If ca=NULL Then Exit Function
      If cend=NULL Then Exit Function
      If nend=NULL Then Exit Function
      AACalpha(i)=asm.Index(ca) : AACend(i)=asm.Index(cend) : AANend(i)=asm.Index(nend)
      Call OxAt(mol,cend,at)    : If at=NULL Then Exit Function
      AAOx(i)=asm.Index(at)
      Call OhxlAt(mol,cend,at)  : If at=NULL Then Exit Function
      AAOhxl(i)=asm.Index(at)
      Call HAt(mol,nend,at)     : If at=NULL Then Exit Function
      AAH1N(i)=asm.Index(at)
    End If

  Next i

  AssignResidueAtoms=True
End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub RecallResidueAtoms(struc  As Object,ByVal k As Integer,ca As Object,cend As Object,nend As Object,ox As Object,ohxl As Object,h1n As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Recall particular atoms for ref. AA #k (k=1-20) and get handles from struc  '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  With struc.Assembly
    ca=.Item(AACalpha(k))
    cend=.Item(AACend(k))
    nend=.Item(AANend(k))
    ox=.Item(AAOx(k))
    ohxl=.Item(AAOhxl(k))
    h1n=.Item(AAH1N(k))
  End With
  'If ca=NULL Or cend=NULL Or nend=NULL Or ox=NULL Or ohxl=NULL Or h1n=NULL Then Exit Function
End Sub






'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function BuildThePeptide(IPeptide() As Integer,ByRef NResidues As Integer,PhiAts() As Object,PsiAts() As Object,diag As Object) As Object
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Building the peptide 3D structure going from N- to C-terminus               '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,ans,nn,AANum As Integer, Peptide, NextAA As Object
Dim ca,nend,cend,ox,ohxl,h1n,ca1,nend1,cend1,ox1,ohxl1,h1n1 As Object


  BuildThePeptide=NULL

  AANum=IPeptide(1)
  Peptide=AAStructure(AANum).Assembly.CopyTree(FALSE).Structures.Item(1)
  Call RecallResidueAtoms(Peptide,AANum,ca,cend,nend,ox,ohxl,h1n)

  'Save atoms forming (Psi,Phi) torsions                                        '
  PhiAts(1,1)=h1n : PhiAts(2,1)=nend : PhiAts(3,1)=ca   : PhiAts(4,1)=cend
  PsiAts(1,1)=nend   : PsiAts(2,1)=ca   : PsiAts(3,1)=cend : PsiAts(4,1)=ox
  'TRICK: our PHI AND PSI defined not as in IUPAC-IUB,                          '
  'but the torsions come to 0 at same values (IUP-IUB consider trans as 0 deg :)'
  '                                                                             '
  ' Suppose we have                                                             '
  '                                                                             '
  ' Cend(i-1) -   N(i)       -   Calpha(i) -  Cend(i)  -  Ocarboxyl(i)          '
  '                                                                             '
  'Our agreement is as follows:                                                 '
  '    phi =   Cend(i-1) -   N(i)       -   Calpha(i) -  Cend(i)                '
  '    psi =   Nend(i)   -   Calpha(i)  -   Cend(i)   -  Ocarboxyl(i)           '
  '                                                                             '
  'Another trick: pass H bound to N(i) at place of Cend(i-1) for 1st residue    '



  'Grow the chain
  nn=1
  For i=2 To Nresidues

    AANum=IPeptide(i)
    NextAA=AAStructure(AANum).Assembly.CopyTree(FALSE).Structures.Item(1)

    Call RecallResidueAtoms(NextAA,AANum,ca1,cend1,nend1,ox1,ohxl1,h1n1)
    PhiAts(1,i)=cend   : PhiAts(2,i)=nend1 : PhiAts(3,i)=ca1   : PhiAts(4,i)=cend1
    PsiAts(1,i)=nend1  : PsiAts(2,i)=ca1   : PsiAts(3,i)=cend1 : PsiAts(4,i)=ox1

    'Extend the chain at C-terminus...
    Call AddAA(Peptide,NextAA,ca,cend,nend,ox,ohxl,h1n,ca1,cend1,nend1,ox1,ohxl1,h1n1)
    'and update a growing end of the chain
    ca=ca1 : nend=nend1 : cend=cend1 : ox=ox1 : ohxl=ohxl1 : h1n=h1n1

    'Refresh a screen if necessary
    nn=nn+1
    if nn>=NResToShowInfo Then
      nn=0
      'ShowStruct(diag,Peptide)
      ans=MessageBox("Building in progress: "+Str(i)+" of "+Str(NResidues)+Chr(13)+"Continue?",TITLE, MBB_YESNO + MBI_QUESTION)
      If ans=MBR_NO Then
        BuildThePeptide=Peptide
        NResidues=i
        Exit Function
      End If
    End If

  Next i


  BuildThePeptide=Peptide

End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub AddAA(ByRef Peptide As Object,AA As Object,ca As Object,cend As Object,nend As Object,ox As Object,ohxl As Object,h1n As Object,ca1 As Object,cend1 As Object,nend1 As Object,ox1 As Object,ohxl1 As Object,h1n1 As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Add amino acid residue to the chain                                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim H_nend,O_cend1,H,Ox1 As Object


  'Add the new residue "mechanically"
  Call Peptide.Assembly.Merge(AA)


  'Make a bond between the chain and new residue
  ReplaceGroup(Peptide, AA, cend, ohxl, nend1, h1n1)


  If (ca<>NULL And ca1<>NULL And nend<>NULL And cend1<>NULL) Then

    'Set Pauling-Corey geometry for peptide bond
    Peptide.SetBLen(cend,nend1,1.32)
    Peptide.SetVAngle(cend,nend1,ca1,126.0/RAD_TO_DEG)
    Peptide.SetVAngle(ca,cend,nend1,118.0/RAD_TO_DEG)
    Peptide.SetVAngle(ox,cend,nend1,123.0/RAD_TO_DEG)

    ' Ensure that peptide bond is trans
    SmartSetTAngle(ca1,nend1,cend,ca,Peptide,180.0/RAD_TO_DEG)
      'Ough... note that the order of atoms is significant if we treat proline
      '(see comments in SmartSetTAngle body


    ' Recall that in the reference structure we have amiNe and H was
    ' at pyramidal nitrogen; now ensure that we have planar amiDe nitrogen and
    ' reposition H
    Call HAt(Peptide, nend1, H) 'could not simply use h1n1 - we've just removed it
    If H<> NULL Then
      Peptide.SetTAngle(H,nend1,cend,ca,0.0)                ' H-N-C"-Ca" --->cis
      Peptide.SetVAngle(H,nend1,ca1,113.0/RAD_TO_DEG)       ' Set Pauling-Corey geometry
    Else
      'a rather special case of proline
    End If

  End If


End Sub




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function ReplaceGroup(Entity As Object, Subst As Object, atE As Object, repE As Object, atS As Object, repS As Object) As Boolean
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Appends Subst to Entity by replacing the group                              '
' at Entity's atom repE with Subst                                            '
'                                                                             '
' Atoms atE and atS will become bonded in (thus grown up) Entity              '
' Atoms repE and repS will be deleted with all their terminal atoms           '
'                                                                             '
' Bending angles at atoms atE and atS keep the source values                  '
' Torsional angle around atE-atS will have arbitrary value                    '
' If bondlen>0 then bond atE-atS length will be equal to bondlen              '
'                                                                             '
'NB: Entity and Subst should be both either molecules or structures           '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim EType As Integer, OK As Boolean, EMol, Bnd As Object
Dim xre,yre,zre,xrs,yrs,zrs,xas,yas,zas,dx,dy,dz,D As Double

  ReplaceGroup = False


  With Entity

    'Check types
    EType = .GetType
    If Subst.GetType <> EType Then Exit Function
    Select Case EType
      Case CB_MOLECULE
        EMol = Entity
      Case CB_STRUCTURE
        EMol = .Molecule
      Case Else
        Exit Function
    End Select

    Call .Assembly.Merge(Subst)           'append
    If EType = CB_STRUCTURE Then          'if we deal with structure do some movements
      Bnd = EMol.AddBond(atE, repS, 1)    'make a bond
      Entity.SetBLen(atE, repS, 0.0)

      'If we deal with structure adjust a valence dash
      Entity.GetAtomXyz(atS,xas,yas,zas)
      Entity.GetAtomXyz(repS,xrs,yrs,zrs)
      Entity.GetAtomXyz(repE,xre,yre,zre)
      dx=(yas-yrs)*(zre-zrs)-(yre-yrs)*(zas-zrs)
      dy=(zas-zrs)*(xre-xrs)-(zre-zrs)*(xas-xrs)
      dz=(xas-xrs)*(yre-yrs)-(xre-xrs)*(yas-yrs)
      D=Sqrt(dx*dx+dy*dy+dz*dz)
      If D>1e-10 Then 'if D<1e-10 then valence position is adjusted
        dx=dx/D: dy=dy/D: dz=dz/D
        Entity.SetAtomXyz(repS,xrs+dx,yrs+dy,zrs+dz)
        Entity.SetTAngle(atS, repS, atE, repE, 0.0)
      End If
      Kill(Bnd)
    End If


    'Delete leaving groups
    Call DeleteTerminalGroup(Entity, repE)
    Call DeleteTerminalGroup(Entity, repS)
    Bnd = EMol.AddBond(atE, atS, 1)   'Make a bond
    ReplaceGroup = TRUE

  End With


End Function 'ReplaceGroup




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub DeleteTerminalGroup(Entity As Object, Terminal As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Deletes Terminal atom and all associated single atoms from Entity           '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim Asm As Object, Environment As Object, Atom As Object

  With Entity

    Asm=.Assembly
    Environment =.Molecule.AssocAtoms(Terminal)

    'Delete associated single atoms
    For Each Atom In Environment
      If .AssocAtoms(Atom).Count = 1 Then Call Asm.AtRemove(Asm.Index(Atom))
    Next Atom

    'Delete central atom itself
    Call Asm.AtRemove(Asm.Index(Terminal))

  End With


End Sub




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function AssignCAlpha(AAMol As Object,ByRef nend As Object,ByRef cend As Object) As Object
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Returns CAlpha atom of AAMol or NULL @fail                                  '
' Also finds C-end and N-end atoms                                            '
' CAlpha is a carbon connected to both NHR and C(O)OR                         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' I do not like this procedure; it should be much more short/nice.            '
' but it works...                                                             '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim asm,at,assoc,neighbor,bonds,b,b1,b2,ab1,ab2 As Object
Dim n, n1, bo1,bo2,bo3 As Integer, Found,FoundCOO,FoundNH,cond1,cond2 As Boolean

  AssignCalpha=NULL
  nend=NULL
  cend=NULL
  asm=AAMol.Assembly


  'Check each atom for probable being Calpha
  For Each at In asm
    FoundCOO=False
    FoundNH=False

    'Check for being carbon of proper coord. number
    If at.GetElNumber<>6 Then Goto nextat
    assoc=AAMol.AssocAtoms(at)
    If (assoc.Count<>4) Then Goto nextat

    'Likely CAlpha, check the neighbours
    For Each neighbor In assoc

      'Check for carbon of C(O)O
      If Not FoundCOO Then
        'Check for carbon of proper coord. number
        If (neighbor.GetElNumber=6) Then
          bonds=AAMol.AssocBonds(neighbor)
          If (bonds.Count<>3) Then Goto nextnbr
          'Remove bond to at (possible CAlpha) from collection
          For Each b In bonds
            If (GetOtherAtom(b,neighbor)=at) Then
              bonds.AtRemove(bonds.Index(b)) : Exit For
            End If
          Next b
          'Check two other bonds
          b1=bonds.Item(1)
          b2=bonds.Item(2)
          bo1=b1.GetBondOrder
          bo2=b2.GetBondOrder
          ab1=GetOtherAtom(b1,neighbor)
          ab2=GetOtherAtom(b2,neighbor)
          cond1 = (bo1=BO_DOUBLE) And (ab1.GetElNumber=8)
          cond1 = cond1 And (bo2=BO_SINGLE) And (ab2.GetElNumber=8)
          cond2 = (bo2=BO_DOUBLE) And (ab2.GetElNumber=8)
          cond2 = cond2 And (bo1=BO_SINGLE) And (ab1.GetElNumber=8)
          FoundCOO = cond1 Or cond2

          'OK?
          If FoundCOO Then
            cend=neighbor
            GoTo nextnbr
          End If
        End If
      End If

      'Check for NHR
      If Not FoundNH Then
        'Check for nitrogen of proper coord. number
        If (neighbor.GetElNumber=7) Then
          bonds=AAMol.AssocBonds(neighbor)
          'Check for all single bonds to neighbours
          FoundNH=True
          For Each b In bonds
            FoundNH=FoundNH And (b.GetBondOrder=BO_SINGLE)
          Next b
          'OK?
          If FoundNH Then
            nend=neighbor
            GoTo nextnbr
          End If
        End If
      End If

      nextnbr:
      If (FoundCOO And FoundNH) Then Exit For
    Next neighbor

    nextat:
    If (FoundCOO And FoundNH) Then
      AssignCAlpha=at
      Exit Function
    End If

  Next at


End Function 'AssignCalpha




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function Calpha(AANum As Integer,ByRef nend As Object,ByRef cend As Object) As Object
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Returns Calpha atom of AA reference unit no. AANum or NULL @fail            '
' Also finds C-end and N-end atoms                                            '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  With AAStructure(AANum).Assembly
    Calpha=.Item(AACalpha(AANum))
    cend=.Item(AACend(AANum))
    nend=.Item(AANend(AANum))
  End With
End Function 'CAlpha




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub HAt(Entity As Object, CentralAtom As Object, Hydrogen As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Returns a single hydrogen atached to CentralAtom or Null if no H present    '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim Atom As Object, Environment As Object

  Environment=Entity.AssocAtoms(CentralAtom)

  For Each Atom In Environment

    If Atom.GetElNumber = 1 Then
      Hydrogen = Atom
      Exit Sub
    End If

  Next Atom

  Hydrogen = NULL  'found nothing
End Sub




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub OhxlAt(Entity As Object, CentralAtom As Object, Oxygen As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Returns oxygen atom of hydroxyl attached to Atom in Entity or Null          '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim Atom As Object, AtomAtO As Object, Environment As Object, OEnvironment As Object

  Environment = Entity.AssocAtoms(CentralAtom)

  For Each Atom In Environment

    If Atom.GetElNumber = 8 Then              'Is oxygen there?
      OEnvironment=Entity.AssocAtoms(Atom)

      For Each AtomAtO In OEnvironment        'Is there H's at it?
        If (AtomAtO <> Atom) And (AtomAtO.GetElNumber = 1) Then
          Oxygen = Atom : Exit Sub
        End If
      Next AtomAtO

    End If

  Next Atom

  Oxygen = NULL  'found nothing
End Sub 'OhxlAt





'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub OxAt(Entity As Object, CentralAtom As Object, Oxygen As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Returns oxygen atom of carbony centred at Atom in Entity or Null            '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim Bond, BEnvironment, At As Object

  BEnvironment = Entity.AssocBonds(CentralAtom)

  For Each Bond In BEnvironment

    If Bond.GetBondOrder=BO_DOUBLE Then

      At=GetOtherAtom(Bond,CentralAtom)
      If At.GetElNumber = 8 Then
        Oxygen=At
        Exit Sub
      End If

    End If

  Next Bond

  Oxygen = NULL  'found nothing
End Sub 'OhxlAt




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function GetOtherAtom(Bond As Object, Atom1 As Object) As Object
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  GetOtherAtom=NULL
  If Bond.Atom1=Atom1 Then
      GetOtherAtom=Bond.Atom2
  Else
    If Bond.Atom2=Atom1 Then GetOtherAtom=Bond.Atom1
  End If
End Function




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub SetOverallSecondaryStructure(Peptide As Object,NResidues As Integer,PhiAts()  As Object,PsiAts() As Object,diag As Object,MotifCode() As Integer)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim i,ans,nn As Integer,phi0,psi0 As Double


  'Go trough a chain
  nn=0
  For i=1 To Nresidues
    Select Case MotifCode(i)
      Case RIGHT_ALPHA_HELIX
                              phi0=-60.0
                              psi0=-60.0
      Case BETA_EXTENDED
                              phi0=-90.0
                              psi0=120.0
      Case BETA_TURN_ANY
                              phi0=-90.0
                              psi0=-30.0
    End Select

    'correct for our vs. IUPAC-IUB definitions
    phi0=phi0
    psi0=-psi0
                'Anybody more experiencied -- check and correct if necessary?


    phi0=phi0/RAD_TO_DEG
    psi0=psi0/RAD_TO_DEG
    SmartSetTAngle(PsiAts(1,i),PsiAts(2,i),PsiAts(3,i),PsiAts(4,i),Peptide,psi0)
    SmartSetTAngle(PhiAts(1,i),PhiAts(2,i),PhiAts(3,i),PhiAts(4,i),Peptide,phi0)

    nn=nn+1
    if nn>=NResToShowInfo Then
      nn=0
      ans=MessageBox("Shaping in progress: "+Str(i)+" of "+Str(NResidues)+Chr(13)+"Continue?",TITLE, MBB_YESNO + MBI_QUESTION)
      If ans=MBR_NO Then Exit Sub
    End If

  Next i

End Sub




'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub ShowStruct(diag As Object,Peptide As Object)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Dim l,t,w,h As Integer
  diag.Depict(Peptide)
  diag.GetBound(l,t,w,h)
  diag.SetBound(10,10,w,h)
End Sub



'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub SmartSetTAngle(a1 As Object,a2 As Object,a3 As Object,a4 As Object,struc As Object,ByVal v As Double)
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'                                                                             '
' This procedure sets TAngle in somewhat more usual way than SetTangle.       '
'                                                                             '
' To clarify: consider a1-a2-a3-a4-R                                          '
'                               \                                             '
'                                a5-R"                                        '
' First, SetTAngle(a1,a2,a3,a4) will rotate a1 leaving a2-a3-a4 in place;     '
'   So,  use SetTAngle(a4,a3,a2,a1) if you wish rotating a4                   '
'   Fortunately, it is not the matter at all in most cases                    '
' Second. At calling SetTAngle(a4,a3,a2,a1), a4 and its chain R               '
'   is rotated at a2-a3; however, a5 and its chain R" will not!               '
'   So, one MUST call SetTAngle(a1,a2,a3,a5) explicitly -- otherwise          '
'   the bonds directionality at a3 will be corrupted.                         '
'                                                                             '
' SmartSetTangle rotates all chains; it also presumes inverse                 '
' (==common) order of passed atoms, a1-a2-a3-a4 for the scheme above.         '
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

Dim a5 As Object,dv As Double
  dv=v-Struc.GetTAngle(a4,a3,a2,a1)
  Struc.SetTAngle(a4,a3,a2,a1,v)
  For Each a5 In struc.Molecule.AssocAtoms(a3)
    If a5<>a4 Then Struc.SetTAngle(a5,a3,a2,a1,Struc.GetTAngle(a5,a3,a2,a1)+dv)
  Next a5
End Sub 'SmartSetTAngle

'***LIBRARY PROCEDURES BEGIN
'@@@@@@