home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Educational / RasMol / Source / molecule.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-10-12  |  48.2 KB  |  2,039 lines
  1. /* molecule.c
  2.  * RasMol2 Molecular Graphics
  3.  * Roger Sayle, August 1995
  4.  * Version 2.6
  5.  */
  6. #include "rasmol.h"
  7.  
  8. #ifdef IBMPC
  9. #include <windows.h>
  10. #include <malloc.h>
  11. #endif
  12. #ifdef APPLEMAC
  13. #include <Types.h>
  14. #endif
  15. #ifndef sun386
  16. #include <stdlib.h>
  17. #endif
  18.  
  19. #include <string.h>
  20. #include <ctype.h>
  21. #include <stdio.h>
  22. #include <math.h>
  23.  
  24. #define MOLECULE
  25. #include "molecule.h"
  26. #include "command.h"
  27. #include "abstree.h"
  28. #include "transfor.h"
  29. #include "render.h"
  30.  
  31. #define HBondPool   32
  32. #define BondPool    32
  33. #define AtomPool    32
  34.  
  35. #define NoLadder     0x00
  36. #define ParaLadder   0x01
  37. #define AntiLadder   0x02
  38.  
  39. #define Cos70Deg     0.34202014332567
  40.  
  41. #define MaxHBondDist   ((Long)300*300)
  42. #define MaxBondDist    ((Long)475*475)
  43. #define MinBondDist    ((Long)100*100)
  44. #define AbsMaxBondDist 600
  45.  
  46. #ifdef APPLEMAC
  47. #define AllocSize   256
  48. typedef struct _AllocRef {
  49.     struct _AllocRef *next;
  50.     void *data[AllocSize];
  51.     int count;
  52.     } AllocRef;
  53. static AllocRef *AllocList;  
  54. #endif
  55.  
  56.  
  57. typedef struct {
  58.           char name[4];
  59.           int code;
  60.       } SynonymTable;
  61.  
  62. #define RESSYNMAX 16
  63. static SynonymTable ResSynonym[RESSYNMAX] = {
  64.     { "ADE", 24 },  /*   A : Adenosine   */
  65.     { "CPR", 11 },  /* PRO : Cis-proline */
  66.     { "CSH", 17 },  /* CYS : Cystine     */
  67.     { "CSM", 17 },  /* CYS : Cystine     */
  68.     { "CYH", 17 },  /* CYS : Cystine     */
  69.     { "CYT", 25 },  /*   C : Cytosine    */
  70.     { "D2O", 47 },  /* DOD : Heavy Water */
  71.     { "GUA", 26 },  /*   G : Guanosine   */
  72.     { "H2O", 46 },  /* HOH : Solvent     */
  73.     { "SOL", 46 },  /* HOH : Solvent     */
  74.     { "SUL", 48 },  /* SO4 : Sulphate    */
  75.     { "THY", 27 },  /*   T : Thymidine   */
  76.     { "TIP", 46 },  /* HOH : Water       */
  77.     { "TRY", 20 },  /* TRP : Tryptophan  */
  78.     { "URI", 28 },  /*   U : Uridine     */
  79.     { "WAT", 46 }   /* HOH : Water       */
  80.         };
  81.  
  82.  
  83. static Molecule __far *FreeMolecule;
  84. static HBond __far *FreeHBond;
  85. static Chain __far *FreeChain;
  86. static Group __far *FreeGroup;
  87. static Atom __far *FreeAtom;
  88. static Bond __far *FreeBond;
  89.  
  90. static IntCoord __far *IntPrev;
  91. static HBond __far * __far *CurHBond;
  92. static int MemSize;
  93.  
  94.  
  95. /* Macros for commonly used loops */
  96. #define ForEachAtom  for(chain=Database->clist;chain;chain=chain->cnext) \
  97.              for(group=chain->glist;group;group=group->gnext)    \
  98.              for(aptr=group->alist;aptr;aptr=aptr->anext)
  99. #define ForEachBond  for(bptr=Database->blist;bptr;bptr=bptr->bnext)
  100.  
  101.  
  102. /* Forward Reference */
  103. void DestroyDatabase();
  104.  
  105.  
  106. #ifdef APPLEMAC
  107. /* External RasMac Function Declaration! */
  108. void SetFileInfo( char*, OSType, OSType, short );
  109. #endif
  110.  
  111.  
  112. static void FatalDataError(ptr)
  113.     char *ptr;
  114. {
  115.     char buffer[80];
  116.  
  117.     sprintf(buffer,"Database Error: %s!",ptr);
  118.     RasMolFatalExit(buffer);
  119. }
  120.  
  121.  
  122.  
  123. void DescribeMolecule()
  124. {
  125.     char buffer[40];
  126.  
  127.     if( CommandActive )
  128.     WriteChar('\n');
  129.     CommandActive=False;
  130.  
  131.     if( *InfoMoleculeName )
  132.     {   WriteString("Molecule name ....... ");
  133.     WriteString(InfoMoleculeName);
  134.     WriteChar('\n');
  135.     }
  136.  
  137.     if( *InfoClassification )
  138.     {   WriteString("Classification ...... ");
  139.     WriteString(InfoClassification);
  140.     WriteChar('\n');
  141.     }
  142.  
  143.     if( Database && (MainGroupCount>1) )
  144.     {   WriteString("Secondary Structure . ");
  145.     if( InfoStrucSrc==SourceNone )
  146.     {   WriteString("No Assignment\n");
  147.     } else if( InfoStrucSrc==SourcePDB )
  148.     {   WriteString("PDB Data Records\n");
  149.     } else WriteString("Calculated\n");
  150.     }
  151.  
  152.  
  153.     if( *InfoIdentCode )
  154.     {   WriteString("Brookhaven Code ..... ");
  155.     WriteString(InfoIdentCode);
  156.     WriteChar('\n');
  157.     }
  158.  
  159.     if( InfoChainCount>1 )
  160.     {   sprintf(buffer,"Number of Chains .... %d\n",InfoChainCount);
  161.     WriteString(buffer);
  162.     }
  163.  
  164.     sprintf(buffer,"Number of Groups .... %d",MainGroupCount);
  165.     WriteString(buffer);
  166.     if( HetaAtomCount )
  167.     {   sprintf(buffer," (%d)\n",HetaGroupCount);
  168.     WriteString(buffer);
  169.     } else WriteChar('\n');
  170.  
  171.     sprintf(buffer,"Number of Atoms ..... %ld",(long)MainAtomCount);
  172.     WriteString(buffer);
  173.     if( HetaAtomCount )
  174.     {   sprintf(buffer," (%d)\n",HetaAtomCount);
  175.     WriteString(buffer);
  176.     } else WriteChar('\n');
  177.  
  178.     if( InfoBondCount )
  179.     {   sprintf(buffer,"Number of Bonds ..... %ld\n",(long)InfoBondCount);
  180.     WriteString(buffer);
  181.     }
  182.  
  183.     if( InfoSSBondCount != -1 )
  184.     {   WriteString("Number of Bridges ... ");
  185.     sprintf(buffer,"%d\n\n",InfoSSBondCount);
  186.     WriteString(buffer);
  187.     }
  188.  
  189.     if( InfoHBondCount != -1 )
  190.     {   WriteString("Number of H-Bonds ... ");
  191.     sprintf(buffer,"%d\n",InfoHBondCount);
  192.     WriteString(buffer);
  193.     }
  194.  
  195.     if( InfoHelixCount != -1 )
  196.     {   WriteString("Number of Helices ... ");
  197.     sprintf(buffer,"%d\n",InfoHelixCount);
  198.     WriteString(buffer);
  199.  
  200.     WriteString("Number of Strands ... ");
  201.     sprintf(buffer,"%d\n",InfoLadderCount);
  202.     WriteString(buffer);
  203.  
  204.     WriteString("Number of Turns ..... ");
  205.     sprintf(buffer,"%d\n",InfoTurnCount);
  206.     WriteString(buffer);
  207.     }
  208. }
  209.  
  210.  
  211. #ifdef APPLEMAC
  212. /* Avoid System Memory Leaks! */
  213. static void RegisterAlloc( data )
  214.     void *data;
  215. {
  216.     register AllocRef *ptr;
  217.     
  218.     if( !AllocList || (AllocList->count==AllocSize) )
  219.     {   ptr = (AllocRef *)_fmalloc( sizeof(AllocRef) );
  220.     if( !ptr ) FatalDataError("Memory allocation failed");
  221.     
  222.     ptr->next = AllocList;
  223.     ptr->data[0] = data;
  224.     ptr->count = 1;
  225.     AllocList = ptr;
  226.     } else AllocList->data[AllocList->count++] = data;
  227. }
  228. #else
  229. #define RegisterAlloc(x)
  230. #endif
  231.  
  232.  
  233. /*==================================*/
  234. /* Group & Chain Handling Functions */
  235. /*==================================*/
  236.  
  237. void CreateChain( ident )
  238.     int ident;
  239. {
  240.     register Chain __far *prev;
  241.  
  242.     if( !CurMolecule )
  243.     {   if( !(CurMolecule = FreeMolecule) )
  244.     {   MemSize += sizeof(Molecule);
  245.         CurMolecule = (Molecule __far *)_fmalloc(sizeof(Molecule));
  246.         if( !CurMolecule ) FatalDataError("Memory allocation failed");
  247.         RegisterAlloc( CurMolecule );
  248.     } else FreeMolecule = (void __far*)0;
  249.  
  250.     CurChain = (void __far*)0;
  251.     CurMolecule->slist = (void __far*)0;
  252.     CurMolecule->hlist = (void __far*)0;
  253.     CurMolecule->blist = (void __far*)0;
  254.     CurMolecule->clist = (void __far*)0;
  255.     Database = CurMolecule;
  256.     }
  257.  
  258.     /* Handle chain breaks! */
  259.     if( !(prev=CurChain) )
  260.     if( (prev=CurMolecule->clist) )
  261.         while( prev->cnext )
  262.         prev = prev->cnext;
  263.  
  264.     if( !(CurChain = FreeChain) )
  265.     {   MemSize += sizeof(Chain);
  266.     CurChain = (Chain __far *)_fmalloc(sizeof(Chain));
  267.     if( !CurChain ) FatalDataError("Memory allocation failed");
  268.     RegisterAlloc( CurChain );
  269.     } else FreeChain = FreeChain->cnext;
  270.  
  271.     if( prev )
  272.     {   prev->cnext = CurChain;
  273.     } else CurMolecule->clist = CurChain;
  274.     CurChain->cnext = (void __far*)0;
  275.      
  276.     CurChain->ident = ident;
  277.     CurChain->model = NMRModel;
  278.     CurChain->glist = (void __far*)0;
  279.     CurChain->blist = (void __far*)0;
  280.     ConnectAtom = (void __far*)0;
  281.     CurGroup = (void __far*)0;
  282.     InfoChainCount++;
  283. }
  284.  
  285.  
  286. void CreateGroup( pool )
  287.     int pool;
  288. {
  289.     register Group __far *ptr;
  290.     register int i;
  291.  
  292.     if( !(ptr = FreeGroup) )
  293.     {   MemSize += pool*sizeof(Group);
  294.     ptr = (Group __far *)_fmalloc( pool*sizeof(Group) );
  295.     if( !ptr ) FatalDataError("Memory allocation failed");
  296.     RegisterAlloc( ptr );
  297.     for( i=1; i<pool; i++ )
  298.     {   ptr->gnext = FreeGroup;
  299.         FreeGroup = ptr++;
  300.     } 
  301.     } else FreeGroup = ptr->gnext;
  302.     
  303.     if( CurGroup )
  304.     {   ptr->gnext = CurGroup->gnext;
  305.         CurGroup->gnext = ptr;
  306.     } else 
  307.     {   ptr->gnext = CurChain->glist;
  308.         CurChain->glist = ptr;
  309.     }
  310.     CurGroup = ptr;
  311.  
  312.     CurAtom = (void __far*)0;
  313.     ptr->alist = (void __far*)0;
  314.     ptr->insert = ' ';
  315.     ptr->struc = 0;
  316.     ptr->flag = 0;
  317.     ptr->col1 = 0;
  318.     ptr->col2 = 0;
  319. }
  320.  
  321.  
  322. int FindResNo( ptr )
  323.     char *ptr;
  324. {
  325.     register int hi,lo;
  326.     register int refno;
  327.     register int flag;
  328.     register int mid;
  329.  
  330.     for( refno=0; refno<ResNo; refno++ )
  331.     if( !strncmp(Residue[refno],ptr,3) )
  332.         return( refno );
  333.  
  334.     lo = 0;
  335.     hi = RESSYNMAX;
  336.     while( lo < hi )
  337.     {   mid = (hi+lo)>>1;
  338.         flag = strncmp(ResSynonym[mid].name,ptr,3);
  339.         if( !flag ) return( ResSynonym[mid].code );
  340.  
  341.         /* Binary Search */
  342.         if( flag<0 )
  343.         {   lo = mid+1;
  344.         } else hi = mid;
  345.     }
  346.  
  347.     if( ResNo++ == MAXRES )
  348.     FatalDataError("Too many new residues");
  349.     Residue[refno][0] = *ptr++;
  350.     Residue[refno][1] = *ptr++;
  351.     Residue[refno][2] = *ptr;
  352.     return( refno );
  353. }
  354.  
  355.  
  356. void ProcessGroup( heta )
  357.     int heta;
  358. {
  359.     register int serno;
  360.  
  361.     serno = CurGroup->serno;
  362.     if( IsSolvent(CurGroup->refno) )
  363.         heta = True;
  364.  
  365.     if( heta )
  366.     {   HetaGroupCount++;
  367.         if( HMinMaxFlag )
  368.         {   if( serno > MaxHetaRes )
  369.             {   MaxHetaRes = serno;
  370.             } else if( serno < MinHetaRes )
  371.                 MinHetaRes = serno;
  372.         } else MinHetaRes = MaxHetaRes = serno;
  373.     } else
  374.     {   MainGroupCount++;
  375.         if( MMinMaxFlag )
  376.         {   if( serno > MaxMainRes )
  377.             {   MaxMainRes = serno;
  378.             } else if( serno < MinMainRes )
  379.                 MinMainRes = serno;
  380.         } else MinMainRes = MaxMainRes = serno;
  381.     }
  382. }
  383.  
  384.  
  385. void CreateMolGroup()
  386. {
  387.     strcpy(InfoFileName,DataFileName);
  388.  
  389.     CreateChain( ' ' );
  390.     CreateGroup( 1 );
  391.  
  392.     CurGroup->refno = FindResNo( "MOL" );
  393.     CurGroup->serno = 1;
  394.     
  395.     MinMainRes = MaxMainRes = 1;
  396.     MinHetaRes = MaxHetaRes = 0;
  397.     MainGroupCount = 1;
  398. }
  399.  
  400.  
  401. /*=========================*/
  402. /* Atom Handling Functions */
  403. /*=========================*/
  404.  
  405.  
  406. Atom __far *CreateAtom()
  407. {
  408.     register Atom __far *ptr;
  409.     register int i;
  410.  
  411.     if( !(ptr = FreeAtom) )
  412.     {   MemSize += AtomPool*sizeof(Atom);
  413.     ptr = (Atom __far *)_fmalloc( AtomPool*sizeof(Atom) );
  414.     if( !ptr ) FatalDataError("Memory allocation failed");
  415.     RegisterAlloc( ptr );
  416.     for( i=1; i<AtomPool; i++ )
  417.     {   ptr->anext = FreeAtom;
  418.         FreeAtom = ptr++;
  419.     } 
  420.     } else FreeAtom = ptr->anext;
  421.  
  422.     if( CurAtom )
  423.     {   ptr->anext = CurAtom->anext;
  424.         CurAtom->anext = ptr;
  425.     } else 
  426.     {   ptr->anext = CurGroup->alist;
  427.         CurGroup->alist = ptr;
  428.     }
  429.     CurAtom = ptr;
  430.  
  431.     SelectCount++;
  432.     ptr->flag = SelectFlag | NonBondFlag;
  433.     ptr->label = (void*)0;
  434.     ptr->radius = 375;
  435.     ptr->altl = ' ';
  436.     ptr->mbox = 0;
  437.     ptr->col = 0;
  438.  
  439.     return( ptr );
  440. }
  441.  
  442.  
  443. void ProcessAtom( ptr )
  444.     Atom __far *ptr;
  445. {
  446.     ptr->elemno = GetElemNumber(CurGroup,ptr);
  447.     if( ptr->elemno == 1 )
  448.     {   ptr->flag |= HydrogenFlag;
  449.     HasHydrogen = True;
  450.     }
  451.  
  452.     if( !IsSolvent(CurGroup->refno) )
  453.     {   if( !(ptr->flag&(HydrogenFlag|HeteroFlag)) )
  454.         ptr->flag |= NormAtomFlag;
  455.     } else ptr->flag |= HeteroFlag;
  456.  
  457. #ifdef INVERT
  458.     ptr->yorg = -ptr->yorg;
  459. #endif
  460.  
  461.     if( HMinMaxFlag || MMinMaxFlag )
  462.     {   if( ptr->xorg < MinX ) 
  463.     {   MinX = ptr->xorg;
  464.     } else if( ptr->xorg > MaxX ) 
  465.         MaxX = ptr->xorg;
  466.  
  467.     if( ptr->yorg < MinY ) 
  468.     {   MinY = ptr->yorg;
  469.     } else if( ptr->yorg > MaxY ) 
  470.         MaxY = ptr->yorg;
  471.  
  472.     if( ptr->zorg < MinZ ) 
  473.     {   MinZ = ptr->zorg;
  474.     } else if( ptr->zorg > MaxZ ) 
  475.         MaxZ = ptr->zorg;
  476.     } else 
  477.     {   MinX = MaxX = ptr->xorg;
  478.     MinY = MaxY = ptr->yorg;
  479.     MinZ = MaxZ = ptr->zorg;
  480.     }
  481.         
  482.     if( ptr->flag & HeteroFlag )
  483.     {   if( HMinMaxFlag )
  484.     {   if( ptr->temp < MinHetaTemp ) 
  485.         {   MinHetaTemp = ptr->temp;
  486.         } else if( ptr->temp > MaxHetaTemp ) 
  487.         MaxHetaTemp = ptr->temp;
  488.     } else MinHetaTemp = MaxHetaTemp = ptr->temp;
  489.     HMinMaxFlag = True;
  490.     HetaAtomCount++;
  491.     } else
  492.     {   if( MMinMaxFlag )
  493.     {   if( ptr->temp < MinMainTemp ) 
  494.         {   MinMainTemp = ptr->temp;
  495.         } else if( ptr->temp > MaxMainTemp ) 
  496.         MaxMainTemp = ptr->temp;
  497.     } else MinMainTemp = MaxMainTemp = ptr->temp;
  498.     MMinMaxFlag = True;
  499.     MainAtomCount++;
  500.     }
  501. }
  502.  
  503.  
  504. Atom __far *FindGroupAtom( group, n )
  505.     Group __far *group;  Byte n;
  506. {
  507.     register Atom __far *ptr;
  508.  
  509.     for( ptr=group->alist; ptr; ptr=ptr->anext )
  510.         if( ptr->refno == n ) return( ptr );
  511.     return( (Atom __far*)0 );
  512. }
  513.  
  514.  
  515. int NewAtomType( ptr )
  516.     char *ptr;
  517. {
  518.     register int refno;
  519.     register int i;
  520.  
  521.     for( refno=0; refno<ElemNo; refno++ )
  522.     if( !strncmp(ElemDesc[refno],ptr,4) )
  523.         return(refno);
  524.  
  525.     if( ElemNo++ == MAXELEM )
  526.         FatalDataError("Too many new atom types");
  527.  
  528.     for( i=0; i<4; i++ )
  529.         ElemDesc[refno][i] = ptr[i];
  530.     return( refno );
  531. }
  532.  
  533.  
  534. int SimpleAtomType( type )
  535.     char *type;
  536. {
  537.     char name[4];
  538.  
  539.     name[2] = name[3] = ' ';
  540.     if( type[1] && (type[1]!=' ') )
  541.     {   name[0] = ToUpper(type[0]);
  542.     name[1] = ToUpper(type[1]);
  543.     } else
  544.     {   name[1] = ToUpper(type[0]);
  545.     name[0] = ' ';
  546.     }
  547.     return( NewAtomType(name) );
  548. }
  549.  
  550.  
  551. int ComplexAtomType( ptr )
  552.     char *ptr;
  553. {
  554.     static char name[4];
  555.     register int i;
  556.  
  557.     if( isdigit(ptr[1]) )
  558.     {   name[1] = ToUpper(*ptr);
  559.     name[2] = name[3] = ' ';
  560.     name[0] = ' ';
  561.     } else for( i=0; i<4; i++ )
  562.     name[i] = ToUpper(ptr[i]);
  563.  
  564.  
  565.     /* Handle Unconventional Naming */
  566.     if( IsProtein(CurGroup->refno) )
  567.     {   if( name[0]=='H' )
  568.         {   name[0]=' ';
  569.             name[1]='H';
  570.         }
  571.     } else if( IsNucleo(CurGroup->refno) )
  572.     {   if( name[3]=='\'' )
  573.             name[3] = '*';
  574.  
  575.         if( (name[1]=='O') && (name[2]=='P') )
  576.         {   if( !strncmp(name," OP1",4) ||
  577.                 !strncmp(name,"1OP ",4) )
  578.                 return( 8 );
  579.             if( !strncmp(name," OP2",4) ||
  580.                 !strncmp(name,"2OP ",4) )
  581.                 return( 9 );
  582.         }
  583.     }
  584.     return( NewAtomType(name) );
  585. }
  586.  
  587.  
  588.  
  589.  
  590.  
  591. /*===============================*/
  592. /* Z-Matrix Conversion Functions */
  593. /*===============================*/
  594.  
  595. #ifdef FUNCPROTO
  596. static IntCoord __far* GetInternalCoord( int );
  597. #endif
  598.  
  599. void InitInternalCoords()
  600. {
  601.     IntList = (IntCoord __far*)0;
  602.     IntPrev = (IntCoord __far*)0;
  603. }
  604.  
  605.  
  606. IntCoord __far* AllocInternalCoord()
  607. {
  608.     register IntCoord __far *ptr;
  609.  
  610.     ptr = (IntCoord __far*)_fmalloc(sizeof(IntCoord));
  611.     if( !ptr ) FatalDataError("Memory allocation failed");
  612.     ptr->inext = (IntCoord __far*)0;
  613.  
  614.     if( IntPrev )
  615.     {   IntPrev->inext = ptr;
  616.     } else IntList = ptr;
  617.     IntPrev = ptr;
  618.     return( ptr );
  619. }
  620.  
  621.  
  622. static IntCoord __far* GetInternalCoord( index )
  623.     int index;
  624. {
  625.     register IntCoord __far *ptr;
  626.  
  627.     ptr = IntList;
  628.     while( (index>1) && ptr->inext )
  629.     {   ptr = ptr->inext;
  630.         index--;
  631.     }
  632.     return( ptr );
  633. }
  634.  
  635.  
  636. void FreeInternalCoords()
  637. {
  638.     register IntCoord __far *ptr;
  639.  
  640.     while( (ptr = IntList) )
  641.     {    IntList = ptr->inext;
  642.          _ffree( ptr );
  643.     }
  644. }
  645.  
  646.  
  647. int ConvertInternal2Cartesian()
  648. {
  649.     register IntCoord __far *ptr;
  650.     register IntCoord __far *na;
  651.     register IntCoord __far *nb;
  652.     register IntCoord __far *nc;
  653.     register double cosph,sinph,costh,sinth,coskh,sinkh;
  654.     register double cosa,sina,cosd,sind;
  655.     register double dist,angle,dihed;
  656.  
  657.     register double xpd,ypd,zpd,xqd,yqd,zqd;
  658.     register double xa,ya,za,xb,yb,zb;
  659.     register double rbc,xyb,yza,temp;
  660.     register double xpa,ypa,zqa;
  661.     register double xd,yd,zd;
  662.     register int flag;
  663.  
  664.  
  665.     /* Atom #1 */
  666.     ptr = IntList;
  667.     ptr->dist  = 0.0;
  668.     ptr->angle = 0.0;
  669.     ptr->dihed = 0.0;
  670.  
  671.     if( !(ptr=ptr->inext) )
  672.         return( True );
  673.  
  674.     /* Atom #2 */
  675.     ptr->angle = 0.0;
  676.     ptr->dihed = 0.0;
  677.  
  678.     if( !(ptr=ptr->inext) )
  679.         return( True );
  680.  
  681.     /* Atom #3 */
  682.     dist = ptr->dist;
  683.     angle = Deg2Rad*ptr->angle;
  684.     cosa = cos(angle);
  685.     sina = sin(angle);
  686.     if( ptr->na == 1 )
  687.     {   na = IntList;
  688.         ptr->dist = na->dist + cosa*dist;
  689.     } else /* ptr->na == 2 */
  690.     {   na = IntList->inext;
  691.         ptr->dist = na->dist - cosa*dist;
  692.     }
  693.     ptr->angle = sina*dist;
  694.     ptr->dihed = 0.0;
  695.  
  696.     while( ptr=ptr->inext )
  697.     {   dist = ptr->dist;
  698.         angle = Deg2Rad*ptr->angle;
  699.         dihed = Deg2Rad*ptr->dihed;
  700.  
  701.         /* Optimise this access?? */
  702.         na = GetInternalCoord(ptr->na);
  703.         nb = GetInternalCoord(ptr->nb);
  704.         nc = GetInternalCoord(ptr->nc);
  705.  
  706.         xb = nb->dist  - na->dist;
  707.         yb = nb->angle - na->angle;
  708.         zb = nb->dihed - na->dihed;
  709.  
  710.         rbc = xb*xb + yb*yb + zb*zb;
  711.         if( rbc < 0.0001 )
  712.             return( False );
  713.         rbc= 1.0/sqrt(rbc);
  714.  
  715.         cosa = cos(angle);
  716.         sina = sin(angle);
  717.  
  718.  
  719.         if( fabs(cosa) >= 0.999999 )
  720.         {   /* Colinear */
  721.             temp = dist*rbc*cosa;
  722.             ptr->dist  = na->dist  + temp*xb;
  723.             ptr->angle = na->angle + temp*yb;
  724.             ptr->dihed = na->dihed + temp*zb;
  725.         } else
  726.         {   xa = nc->dist  - na->dist;
  727.             ya = nc->angle - na->angle;
  728.             za = nc->dihed - na->dihed;
  729.  
  730.             sind = -sin(dihed);
  731.             cosd = cos(dihed);
  732.  
  733.             xd = dist*cosa;
  734.             yd = dist*sina*cosd;
  735.             zd = dist*sina*sind;
  736.  
  737.             xyb = sqrt(xb*xb + yb*yb);
  738.             if( xyb < 0.1 )
  739.             {   /* Rotate about y-axis! */
  740.                 temp = za; za = -xa; xa = temp;
  741.                 temp = zb; zb = -xb; xb = temp;
  742.                 xyb = sqrt(xb*xb + yb*yb);
  743.                 flag = True;
  744.             } else flag = False;
  745.  
  746.             costh = xb/xyb;
  747.             sinth = yb/xyb;
  748.             xpa = costh*xa + sinth*ya;
  749.             ypa = costh*ya - sinth*xa;
  750.  
  751.             sinph = zb*rbc;
  752.             cosph = sqrt(1.0 - sinph*sinph);
  753.             zqa = cosph*za  - sinph*xpa;
  754.  
  755.             yza = sqrt(ypa*ypa + zqa*zqa);
  756.  
  757.             if( yza > 1.0E-10 )
  758.             {   coskh = ypa/yza;
  759.                 sinkh = zqa/yza;
  760.  
  761.                 ypd = coskh*yd - sinkh*zd;
  762.                 zpd = coskh*zd + sinkh*yd;
  763.             } else
  764.             {   /* coskh = 1.0; */
  765.                 /* sinkh = 0.0; */
  766.                 ypd = yd;
  767.                 zpd = zd;
  768.             }
  769.  
  770.             xpd = cosph*xd  - sinph*zpd;
  771.             zqd = cosph*zpd + sinph*xd;
  772.             xqd = costh*xpd - sinth*ypd;
  773.             yqd = costh*ypd + sinth*xpd;
  774.  
  775.             if( flag )
  776.             {   /* Rotate about y-axis! */
  777.                 ptr->dist  = na->dist  - zqd;
  778.                 ptr->angle = na->angle + yqd;
  779.                 ptr->dihed = na->dihed + xqd;
  780.             } else
  781.             {   ptr->dist  = na->dist  + xqd;
  782.                 ptr->angle = na->angle + yqd;
  783.                 ptr->dihed = na->dihed + zqd;
  784.             }
  785.         }
  786.     }
  787.     return( True );
  788. }
  789.  
  790.  
  791.  
  792. /*=========================*/
  793. /* Bond Handling Functions */
  794. /*=========================*/
  795.  
  796.  
  797. #ifdef FUNCPROTO
  798. Bond __far *ProcessBond( Atom __far*, Atom __far*, int );
  799. static void CreateHydrogenBond( Atom __far*, Atom __far*,
  800.                 Atom __far*, Atom __far*, int, int );
  801. #endif
  802.  
  803. Bond __far *ProcessBond( src, dst, flag )
  804.     Atom __far *src, __far *dst;
  805.     int flag;
  806. {
  807.     register Bond __far *ptr;
  808.     register int i;
  809.  
  810.     if( flag & (DoubBondFlag|TripBondFlag) )
  811.         DrawDoubleBonds = True;
  812.  
  813.     if( !(ptr = FreeBond) )
  814.     {   MemSize += BondPool*sizeof(Bond);
  815.     ptr = (Bond __far *)_fmalloc( BondPool*sizeof(Bond) );
  816.     if( !ptr ) FatalDataError("Memory allocation failed");
  817.     RegisterAlloc( ptr );
  818.     for( i=1; i<BondPool; i++ )
  819.     {   ptr->bnext = FreeBond;
  820.         FreeBond = ptr++;
  821.     } 
  822.     } else FreeBond = ptr->bnext;
  823.  
  824.     ptr->flag = flag | SelectFlag;
  825.     ptr->srcatom = src;
  826.     ptr->dstatom = dst;
  827.     ptr->radius = 0;
  828.     ptr->col = 0;
  829.  
  830.     return( ptr );
  831. }
  832.  
  833. static void CreateHydrogenBond( srcCA, dstCA, src, dst, energy, offset )
  834.     Atom __far *srcCA, __far *dstCA;
  835.     Atom __far *src, __far *dst;
  836.     int energy, offset;
  837. {
  838.     register HBond __far *ptr;
  839.     register int i,flag;
  840.  
  841.     if( !(ptr = FreeHBond) )
  842.     {   MemSize += HBondPool*sizeof(HBond);
  843.     ptr = (HBond __far *)_fmalloc( HBondPool*sizeof(HBond) );
  844.     if( !ptr ) FatalDataError("Memory allocation failed");
  845.     RegisterAlloc( ptr );
  846.     for( i=1; i<HBondPool; i++ )
  847.     {   ptr->hnext = FreeHBond;
  848.         FreeHBond = ptr++;
  849.     } 
  850.     } else FreeHBond = ptr->hnext;
  851.  
  852.     if( (offset>=-128) && (offset<127) )
  853.     {   ptr->offset = (Char)offset;
  854.     } else ptr->offset = 0;
  855.  
  856.     flag = ZoneBoth? src->flag&dst->flag : src->flag|dst->flag;
  857.     ptr->flag = flag & SelectFlag;
  858.  
  859.     ptr->src = src;
  860.     ptr->dst = dst;
  861.     ptr->srcCA = srcCA;
  862.     ptr->dstCA = dstCA;
  863.     ptr->energy = energy;
  864.     ptr->col = 0;
  865.  
  866.     *CurHBond = ptr;
  867.     ptr->hnext = (void __far*)0;
  868.     CurHBond = &ptr->hnext;
  869.     InfoHBondCount++;
  870. }
  871.  
  872.  
  873. void CreateBond( src, dst, flag )
  874.     int src, dst, flag;
  875. {
  876.     register Chain __far *chain;
  877.     register Group __far *group;
  878.     register Atom __far *aptr;
  879.     register Atom __far *sptr;
  880.     register Atom __far *dptr;
  881.     register Bond __far *bptr;
  882.     register int done;
  883.  
  884.  
  885.     if( src == dst )
  886.     return;
  887.  
  888.     sptr = (void __far*)0;
  889.     dptr = (void __far*)0;
  890.  
  891.     done = False;
  892.     for( chain=Database->clist; chain && !done; chain=chain->cnext )
  893.     for( group=chain->glist; group && !done; group=group->gnext )
  894.         for( aptr=group->alist; aptr; aptr=aptr->anext )
  895.         {   if( aptr->serno == src )
  896.         {   sptr = aptr;
  897.             if( dptr )
  898.             {   done = True;
  899.             break;
  900.             }
  901.         } else if( aptr->serno == dst )
  902.         {   dptr = aptr;
  903.             if( sptr )
  904.             {   done = True;
  905.             break;
  906.             }
  907.         }
  908.         }
  909.  
  910.  
  911.     /* Both found! */
  912.     if( done ) 
  913.     {   if( flag )
  914.         {   /* Reset Non-bonded flags! */
  915.         sptr->flag &= ~NonBondFlag;
  916.         dptr->flag &= ~NonBondFlag;
  917.  
  918.         bptr = ProcessBond( sptr, dptr, flag );
  919.         bptr->bnext = CurMolecule->blist;
  920.         CurMolecule->blist = bptr;
  921.         InfoBondCount++;
  922.  
  923.         } else /* Hydrogen Bond! */
  924.         {   if( InfoHBondCount<0 ) 
  925.             {   CurHBond = &CurMolecule->hlist;
  926.                 InfoHBondCount = 0;
  927.             }
  928.             CreateHydrogenBond( NULL, NULL, sptr, dptr, 0, 0 );
  929.         }
  930.     }
  931. }
  932.  
  933.  
  934. void CreateBondOrder( src, dst )
  935.     int src, dst;
  936. {
  937.     register Bond __far *bptr;
  938.     register int bs,bd;
  939.  
  940.     ForEachBond
  941.     {   bs = bptr->srcatom->serno;
  942.     bd = bptr->dstatom->serno;
  943.  
  944.     if( ((bs==src)&&(bd==dst)) || ((bs==dst)&&(bd==src)) )
  945.         {   DrawDoubleBonds = True;
  946.         if( bptr->flag & NormBondFlag )
  947.         {  /* Convert Single to Double */
  948.            bptr->flag &= ~(NormBondFlag);
  949.            bptr->flag |= DoubBondFlag;
  950.         } else if( bptr->flag & DoubBondFlag )
  951.         {  /* Convert Double to Triple */
  952.            bptr->flag &= ~(DoubBondFlag);
  953.            bptr->flag |= TripBondFlag;
  954.         }
  955.         return;
  956.     }
  957.     }
  958.     CreateBond( src, dst, NormBondFlag );
  959. }
  960.  
  961.  
  962. static void TestBonded( sptr, dptr, flag )
  963.     Atom __far *sptr, __far *dptr; 
  964.     int flag;
  965. {
  966.     register Bond __far *bptr;
  967.     register Long dx, dy, dz;
  968.     register Long max, dist;
  969.  
  970.     if( flag )
  971.     {    /* Sum of covalent radii with 0.56A tolerance */
  972.          dist = Element[sptr->elemno].covalrad + 
  973.                 Element[dptr->elemno].covalrad + 140;
  974.          max = dist*dist;  
  975.     } else 
  976.     {    /* Fast Bio-Macromolecule Bonding Calculation */
  977.          if( (sptr->flag|dptr->flag) & HydrogenFlag )
  978.      {      max = MaxHBondDist;
  979.          } else max = MaxBondDist;
  980.     }
  981.  
  982.     dx = sptr->xorg-dptr->xorg;   if( (dist=dx*dx)>max ) return;
  983.     dy = sptr->yorg-dptr->yorg;   if( (dist+=dy*dy)>max ) return;
  984.     dz = sptr->zorg-dptr->zorg;   if( (dist+=dz*dz)>max ) return;
  985.  
  986.     if( dist > MinBondDist )
  987.     {   /* Reset Non-bonded flags! */
  988.     sptr->flag &= ~NonBondFlag;
  989.     dptr->flag &= ~NonBondFlag;
  990.  
  991.     bptr = ProcessBond(sptr,dptr,NormBondFlag);
  992.     bptr->bnext = CurMolecule->blist;
  993.     CurMolecule->blist = bptr;
  994.     InfoBondCount++;
  995.     }
  996. }
  997.  
  998.  
  999. static void ReclaimHBonds( ptr )
  1000.     HBond __far *ptr;
  1001. {
  1002.     register HBond __far *temp;
  1003.  
  1004.     if( (temp = ptr) )
  1005.     {   while( temp->hnext )
  1006.         temp=temp->hnext;
  1007.     temp->hnext = FreeHBond;
  1008.     FreeHBond = ptr;
  1009.     }
  1010. }
  1011.  
  1012.  
  1013. static void ReclaimBonds( ptr )
  1014.     Bond __far *ptr;
  1015. {
  1016.     register Bond __far *temp;
  1017.  
  1018.     if( (temp = ptr) )
  1019.     {   while( temp->bnext )
  1020.         temp=temp->bnext;
  1021.     temp->bnext = FreeBond;
  1022.     FreeBond = ptr;
  1023.     }
  1024. }
  1025.  
  1026.  
  1027. void CreateMoleculeBonds( info, flag )
  1028.     int info, flag;
  1029. {
  1030.     register int i, x, y, z;
  1031.     register Long tx, ty, tz;
  1032.     register Long mx, my, mz; 
  1033.     register Long dx, dy, dz;
  1034.     register int lx, ly, lz, hx, hy, hz;
  1035.     register Atom __far *aptr, __far *dptr;
  1036.     register Chain __far *chain;
  1037.     register Group __far *group;
  1038.     char buffer[40];
  1039.  
  1040.  
  1041.     if( !Database ) 
  1042.     return;
  1043.  
  1044.     dx = (MaxX-MinX)+1;
  1045.     dy = (MaxY-MinY)+1;
  1046.     dz = (MaxZ-MinZ)+1;
  1047.  
  1048.     InfoBondCount = 0;
  1049.     ReclaimBonds( CurMolecule->blist );
  1050.     CurMolecule->blist = (void __far*)0;
  1051.     ResetVoxelData();
  1052.  
  1053.     for( chain=Database->clist; chain; chain=chain->cnext )
  1054.     {   ResetVoxelData();
  1055.     for( group=chain->glist; group; group=group->gnext )
  1056.         for( aptr=group->alist; aptr; aptr=aptr->anext )
  1057.         {   /* Initially non-bonded! */
  1058.         aptr->flag |= NonBondFlag;
  1059.  
  1060.         mx = aptr->xorg-MinX;
  1061.         my = aptr->yorg-MinY;
  1062.         mz = aptr->zorg-MinZ;
  1063.  
  1064.         tx = mx-AbsMaxBondDist;  
  1065.         ty = my-AbsMaxBondDist;  
  1066.         tz = mz-AbsMaxBondDist;  
  1067.  
  1068.         lx = (tx>0)? (int)((VOXORDER*tx)/dx) : 0;
  1069.         ly = (ty>0)? (int)((VOXORDER*ty)/dy) : 0;
  1070.         lz = (tz>0)? (int)((VOXORDER*tz)/dz) : 0;
  1071.  
  1072.         tx = mx+AbsMaxBondDist;  
  1073.         ty = my+AbsMaxBondDist;  
  1074.         tz = mz+AbsMaxBondDist;
  1075.  
  1076.         hx = (tx<dx)? (int)((VOXORDER*tx)/dx) : VOXORDER-1;
  1077.         hy = (ty<dy)? (int)((VOXORDER*ty)/dy) : VOXORDER-1;
  1078.         hz = (tz<dz)? (int)((VOXORDER*tz)/dz) : VOXORDER-1;
  1079.     
  1080.         for( x=lx; x<=hx; x++ )
  1081.         {   i = VOXORDER2*x + VOXORDER*ly;
  1082.             for( y=ly; y<=hy; y++ )
  1083.             {   for( z=lz; z<=hz; z++ )
  1084.                 if( (dptr = (Atom __far*)HashTable[i+z]) )
  1085.                 do { TestBonded(aptr,dptr,flag);
  1086.                 } while( (dptr = dptr->next) );
  1087.             i += VOXORDER;
  1088.             }
  1089.         }
  1090.         
  1091.         x = (int)((VOXORDER*mx)/dx);
  1092.         y = (int)((VOXORDER*my)/dy);
  1093.         z = (int)((VOXORDER*mz)/dz);
  1094.  
  1095.         i = VOXORDER2*x + VOXORDER*y + z;
  1096.         aptr->next = (Atom __far*)HashTable[i];
  1097.         HashTable[i] = (void __far*)aptr;
  1098.         }
  1099.     VoxelsClean = False;
  1100.     }
  1101.  
  1102.     if( info )
  1103.     {   if( CommandActive )
  1104.         WriteChar('\n');
  1105.     CommandActive=False;
  1106.     sprintf(buffer,"Number of Bonds ..... %ld\n\n",(long)InfoBondCount);
  1107.     WriteString(buffer);
  1108.     }
  1109. }
  1110.  
  1111.  
  1112. /*===============================*/
  1113. /* Disulphide bridging functions */
  1114. /*===============================*/
  1115.  
  1116. #ifdef FUNCPROTO
  1117. static Atom __far *FindCysSulphur(  Group __far* );
  1118. #endif
  1119.  
  1120. static Atom __far *FindCysSulphur( group )
  1121.     Group __far *group;
  1122. {
  1123.     register Atom __far *ptr;
  1124.     register char *elem;
  1125.  
  1126.     for( ptr=group->alist; ptr; ptr=ptr->anext )
  1127.     {   elem = ElemDesc[ptr->refno];
  1128.         if( (elem[1]=='S') && (elem[0]==' ')  )
  1129.             return( ptr );
  1130.     }
  1131.     return( (Atom __far*)0 );
  1132. }
  1133.  
  1134.  
  1135. static void TestDisulphideBridge( group1, group2, cys1 )
  1136.     Group __far *group1, __far *group2;  
  1137.     Atom __far *cys1;
  1138. {
  1139.     register HBond __far *ptr;
  1140.     register Atom __far *cys2;
  1141.     register int dx, dy, dz;
  1142.     register Long max,dist;
  1143.  
  1144.     if( !(cys2=FindCysSulphur(group2)) )
  1145.     return;
  1146.  
  1147.     max = (Long)750*750;
  1148.     dx = (int)(cys1->xorg-cys2->xorg);   if( (dist=(Long)dx*dx)>max ) return;
  1149.     dy = (int)(cys1->yorg-cys2->yorg);   if( (dist+=(Long)dy*dy)>max ) return;
  1150.     dz = (int)(cys1->zorg-cys2->zorg);   if( (dist+=(Long)dz*dz)>max ) return;
  1151.  
  1152.     if( !(ptr = FreeHBond) )
  1153.     {   MemSize += sizeof(HBond);
  1154.     ptr = (HBond __far*)_fmalloc(sizeof(HBond));
  1155.     if( !ptr ) FatalDataError("Memory allocation failed");
  1156.     RegisterAlloc( ptr );
  1157.     } else FreeHBond = ptr->hnext;
  1158.  
  1159.     ptr->dst = cys1;
  1160.     if( !(ptr->dstCA=FindGroupAtom(group1,1)) )
  1161.     ptr->dstCA = cys1;
  1162.  
  1163.     ptr->src = cys2;
  1164.     if( !(ptr->srcCA=FindGroupAtom(group2,1)) )
  1165.     ptr->srcCA = cys2;
  1166.  
  1167.     ptr->offset = 0;
  1168.     ptr->energy = 0;
  1169.     ptr->flag = 0;
  1170.     ptr->col = 0;
  1171.  
  1172.     ptr->hnext = CurMolecule->slist;
  1173.     CurMolecule->slist = ptr;
  1174.  
  1175.     group1->flag |= CystineFlag;
  1176.     group2->flag |= CystineFlag;
  1177.     InfoSSBondCount++;
  1178. }
  1179.  
  1180.  
  1181. void FindDisulphideBridges()
  1182. {
  1183.     register Chain __far *chn1;
  1184.     register Chain __far *chn2;
  1185.     register Group __far *group1;
  1186.     register Group __far *group2;
  1187.     register Atom __far *cys;
  1188.     char buffer[40];
  1189.  
  1190.     if( !Database ) return;
  1191.     ReclaimHBonds( CurMolecule->slist );
  1192.     InfoSSBondCount = 0;
  1193.  
  1194.     for(chn1=Database->clist;chn1;chn1=chn1->cnext)
  1195.     for(group1=chn1->glist;group1;group1=group1->gnext)
  1196.         if( IsCysteine(group1->refno) && (cys=FindCysSulphur(group1)) )
  1197.         {   for(group2=group1->gnext;group2;group2=group2->gnext)
  1198.             if( IsCysteine(group2->refno) )
  1199.             TestDisulphideBridge(group1,group2,cys);
  1200.  
  1201.         for(chn2=chn1->cnext;chn2;chn2=chn2->cnext)
  1202.             for(group2=chn2->glist;group2;group2=group2->gnext)
  1203.             if( IsCysteine(group2->refno) )
  1204.                 TestDisulphideBridge(group1,group2,cys);
  1205.         }
  1206.  
  1207.     if( FileDepth == -1 )
  1208.     {   if( CommandActive )
  1209.             WriteChar('\n');
  1210.         CommandActive=False;
  1211.     
  1212.         sprintf(buffer,"Number of Bridges ... %d\n\n",InfoSSBondCount);
  1213.         WriteString(buffer);
  1214.     }
  1215. }
  1216.  
  1217.  
  1218. /*=========================================*/
  1219. /* Kabsch & Sander Structure Determination */
  1220. /*=========================================*/
  1221.  
  1222. #ifdef FUNCPROTO
  1223. static int CalculateBondEnergy( Group __far* );
  1224. static void CalcProteinHBonds( Chain __far* );
  1225. static void CalcNucleicHBonds( Chain __far* );
  1226. static int IsHBonded( Atom __far*, Atom __far*, HBond __far* );
  1227. static void TestLadder( Chain __far* );
  1228. #endif
  1229.  
  1230.  
  1231. /* Coupling constant for Electrostatic Energy   */
  1232. /* QConst = -332 * 0.42 * 0.2 * 1000.0 [*250.0] */
  1233. #define QConst (-6972000.0)
  1234. #define MaxHDist ((Long)2250*2250)
  1235. #define MinHDist ((Long)125*125)
  1236.  
  1237.  
  1238. /* Protein Donor Atom Coordinates */
  1239. static int hxorg,hyorg,hzorg;
  1240. static int nxorg,nyorg,nzorg;
  1241. static Atom __far *best1CA;
  1242. static Atom __far *best2CA;
  1243. static Atom __far *best1;
  1244. static Atom __far *best2;
  1245. static Atom __far *optr;
  1246. static int res1,res2;
  1247. static int off1,off2;
  1248.  
  1249.  
  1250. static int CalculateBondEnergy( group )
  1251.     Group __far *group;
  1252. {
  1253.     register double dho,dhc;
  1254.     register double dnc,dno;
  1255.  
  1256.     register Atom __far *cptr;
  1257.     register Long dx,dy,dz;
  1258.     register Long dist;
  1259.     register int result;
  1260.  
  1261.     if( !(cptr=FindGroupAtom(group,2)) )  return(0);
  1262.     if( !(optr=FindGroupAtom(group,3)) )  return(0);
  1263.  
  1264.     dx = hxorg-optr->xorg;  
  1265.     dy = hyorg-optr->yorg;  
  1266.     dz = hzorg-optr->zorg;
  1267.     dist = dx*dx+dy*dy+dz*dz;
  1268.     if( dist < MinHDist ) 
  1269.     return( -9900 );
  1270.     dho = sqrt((double)dist);
  1271.  
  1272.     dx = hxorg-cptr->xorg;  
  1273.     dy = hyorg-cptr->yorg;  
  1274.     dz = hzorg-cptr->zorg;
  1275.     dist = dx*dx+dy*dy+dz*dz;
  1276.     if( dist < MinHDist ) 
  1277.     return( -9900 );
  1278.     dhc = sqrt((double)dist);
  1279.  
  1280.     dx = nxorg-cptr->xorg;  
  1281.     dy = nyorg-cptr->yorg;  
  1282.     dz = nzorg-cptr->zorg;
  1283.     dist = dx*dx+dy*dy+dz*dz;
  1284.     if( dist < MinHDist ) 
  1285.     return( -9900 );
  1286.     dnc = sqrt((double)dist);
  1287.  
  1288.     dx = nxorg-optr->xorg;  
  1289.     dy = nyorg-optr->yorg;  
  1290.     dz = nzorg-optr->zorg;
  1291.     dist = dx*dx+dy*dy+dz*dz;
  1292.     if( dist < MinHDist ) 
  1293.     return( -9900 );
  1294.     dno = sqrt((double)dist);
  1295.  
  1296.     result = (int)(QConst/dho - QConst/dhc + QConst/dnc - QConst/dno);
  1297.  
  1298.     if( result<-9900 ) 
  1299.     {   return( -9900 );
  1300.     } else if( result>-500 ) 
  1301.     return( 0 );
  1302.  
  1303.     return( result );
  1304. }
  1305.  
  1306.  
  1307. static void CalcProteinHBonds( chn1 )
  1308.     Chain __far *chn1;
  1309. {
  1310.     register int energy, offset;
  1311.     register Chain __far *chn2;
  1312.     register Group __far *group1;
  1313.     register Group __far *group2;
  1314.     register Atom __far *ca1;
  1315.     register Atom __far *ca2;
  1316.     register Atom __far *pc1;
  1317.     register Atom __far *po1;
  1318.     register Atom __far *n1;
  1319.     register int pos1,pos2;
  1320.     register int dx,dy,dz;
  1321.     register double dco;
  1322.     register Long dist;
  1323.  
  1324.     pos1 = 0;
  1325.     pc1 = po1 = (void __far*)0;
  1326.     for(group1=chn1->glist;group1;group1=group1->gnext)
  1327.     {   pos1++;
  1328.         if( pc1 && po1 )
  1329.         {   dx = (int)(pc1->xorg - po1->xorg);
  1330.             dy = (int)(pc1->yorg - po1->yorg);
  1331.             dz = (int)(pc1->zorg - po1->zorg);
  1332.         } else
  1333.         {   pc1 = FindGroupAtom(group1,2);
  1334.             po1 = FindGroupAtom(group1,3);
  1335.         continue;
  1336.     }
  1337.  
  1338.         pc1 = FindGroupAtom(group1,2);
  1339.         po1 = FindGroupAtom(group1,3);
  1340.  
  1341.         if( !IsAmino(group1->refno) || IsProline(group1->refno) )
  1342.             continue;
  1343.  
  1344.         if( !(ca1=FindGroupAtom(group1,1)) ) continue;
  1345.         if( !(n1=FindGroupAtom(group1,0)) )  continue;
  1346.  
  1347.     dist = (Long)dx*dx + (Long)dy*dy + (Long)dz*dz;
  1348.     dco = sqrt( (double)dist )/250.0;
  1349.  
  1350.     nxorg = (int)n1->xorg;   hxorg = nxorg + (int)(dx/dco);
  1351.     nyorg = (int)n1->yorg;   hyorg = nyorg + (int)(dy/dco);
  1352.     nzorg = (int)n1->zorg;   hzorg = nzorg + (int)(dz/dco);
  1353.     res1 = res2 = 0;
  1354.  
  1355.     /* Only Hydrogen Bond within a single chain!       */
  1356.     /* for(chn2=Database->clist;chn2;chn2=chn2->cnext) */
  1357.  
  1358.     chn2 = chn1;
  1359.     {   /* Only consider non-empty peptide chains! */
  1360.         /* if( !chn2->glist || !IsProtein(chn2->glist->refno) ) */
  1361.         /*     continue;                                        */
  1362.  
  1363.         pos2 = 0;
  1364.         for(group2=chn2->glist;group2;group2=group2->gnext)
  1365.         {   pos2++;
  1366.         if( (group2==group1) || (group2->gnext==group1) )
  1367.             continue;
  1368.  
  1369.         if( !IsAmino(group2->refno) ) 
  1370.             continue;
  1371.         if( !(ca2=FindGroupAtom(group2,1)) ) 
  1372.             continue;
  1373.  
  1374.         dx = (int)(ca1->xorg-ca2->xorg);
  1375.         if( (dist=(Long)dx*dx) > MaxHDist )
  1376.             continue;
  1377.  
  1378.         dy = (int)(ca1->yorg-ca2->yorg);
  1379.         if( (dist+=(Long)dy*dy) > MaxHDist )
  1380.             continue;
  1381.  
  1382.         dz = (int)(ca1->zorg-ca2->zorg);
  1383.         if( (dist+=(Long)dz*dz) > MaxHDist )
  1384.             continue;
  1385.  
  1386.         if( (energy = CalculateBondEnergy(group2)) )
  1387.         {   if( chn1 == chn2 )
  1388.             {   offset = pos1 - pos2;
  1389.             } else offset = 0;
  1390.  
  1391.             if( energy<res1 )
  1392.             {   best2CA = best1CA;  best1CA = ca2;
  1393.             best2 = best1;      best1 = optr;
  1394.             res2 = res1;        res1 = energy;
  1395.             off2 = off1;        off1 = offset;
  1396.             } else if( energy<res2 )
  1397.             {   best2CA = ca2;
  1398.             best2 = optr;
  1399.             res2 = energy;
  1400.             off2 = offset;
  1401.             }
  1402.         }
  1403.         }  /* group2 */
  1404.     }      /* chn2 */
  1405.  
  1406.     if( res1 ) 
  1407.     {   if( res2 ) 
  1408.         CreateHydrogenBond(ca1,best2CA,n1,best2,res2,off2);
  1409.         CreateHydrogenBond(ca1,best1CA,n1,best1,res1,off1);
  1410.     }
  1411.     }
  1412. }
  1413.  
  1414.  
  1415. static void CalcNucleicHBonds( chn1 )
  1416.     Chain __far *chn1;
  1417. {
  1418.     register Chain __far *chn2;
  1419.     register Group __far *group1;
  1420.     register Group __far *group2;
  1421.     register Group __far *best;
  1422.     register Atom __far *ca1;
  1423.     register Atom __far *ca2;
  1424.     register Atom __far *n1;
  1425.     register Long max,dist;
  1426.     register int dx,dy,dz;
  1427.     register int refno;
  1428.  
  1429.  
  1430.     for(group1=chn1->glist;group1;group1=group1->gnext)
  1431.     {   if( !IsPurine(group1->refno) ) continue;
  1432.     /* Find N1 of Purine Group */
  1433.     if( !(n1=FindGroupAtom(group1,21)) )
  1434.         continue;
  1435.  
  1436.     /* Maximum N1-N3 distance 5A */
  1437.     refno = NucleicCompl(group1->refno);
  1438.     max = (Long)1250*1250;
  1439.     best = (void __far*)0;
  1440.  
  1441.     for(chn2=Database->clist;chn2;chn2=chn2->cnext)
  1442.     {   /* Only consider non-empty nucleic acid chains! */
  1443.         if( (chn1==chn2) || !chn2->glist || 
  1444.         !IsDNA(chn2->glist->refno) )
  1445.         continue;
  1446.  
  1447.         for(group2=chn2->glist;group2;group2=group2->gnext)
  1448.         if( group2->refno == (Byte)refno )
  1449.         {   /* Find N3 of Pyramidine Group */
  1450.             if( !(ca1=FindGroupAtom(group2,23)) )
  1451.             continue;
  1452.  
  1453.             dx = (int)(ca1->xorg - n1->xorg);
  1454.             if( (dist=(Long)dx*dx) >= max ) 
  1455.             continue;
  1456.  
  1457.             dy = (int)(ca1->yorg - n1->yorg);
  1458.             if( (dist+=(Long)dy*dy) >= max ) 
  1459.             continue;
  1460.  
  1461.             dz = (int)(ca1->zorg - n1->zorg);
  1462.             if( (dist+=(Long)dz*dz) >= max )
  1463.             continue;
  1464.  
  1465.             best1 = ca1;
  1466.             best = group2;
  1467.             max = dist;
  1468.         }
  1469.     }
  1470.  
  1471.     if( best )
  1472.     {   /* Find the sugar phosphorous atoms */
  1473.         ca1 = FindGroupAtom( group1, 7 );
  1474.         ca2 = FindGroupAtom( best, 7 );
  1475.  
  1476.         CreateHydrogenBond( ca1, ca2, n1, best1, 0, 0 );
  1477.         if( IsGuanine(group1->refno) )
  1478.         {   /* Guanine-Cytosine */
  1479.         if( (ca1=FindGroupAtom(group1,22)) &&  /* G.N2 */
  1480.             (ca2=FindGroupAtom(best,26)) )     /* C.O2 */
  1481.             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1482.                     ca1, ca2, 0, 0 );
  1483.  
  1484.         if( (ca1=FindGroupAtom(group1,28)) &&  /* G.O6 */
  1485.             (ca2=FindGroupAtom(best,24)) )     /* C.N4 */
  1486.             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1487.                     ca1, ca2, 0, 0 );
  1488.  
  1489.         } else /* Adenine-Thymine */
  1490.         if( (ca1=FindGroupAtom(group1,25)) &&  /* A.N6 */
  1491.             (ca2=FindGroupAtom(best,27)) )     /* T.O4 */
  1492.             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1493.                     ca1, ca2, 0, 0 );
  1494.     }
  1495.     }
  1496. }
  1497.  
  1498.  
  1499. void CalcHydrogenBonds()
  1500. {
  1501.     register Chain __far *chn1;
  1502.     char buffer[40];
  1503.  
  1504.     if( !Database ) return;
  1505.     ReclaimHBonds( CurMolecule->hlist );
  1506.     CurMolecule->hlist = (void __far*)0;
  1507.     CurHBond = &CurMolecule->hlist;
  1508.     InfoHBondCount = 0;
  1509.  
  1510.     if( MainAtomCount > 10000 )
  1511.     {   if( CommandActive )
  1512.         WriteChar('\n');
  1513.     WriteString("Please wait... ");
  1514.     CommandActive=True;
  1515.     }
  1516.  
  1517.     for(chn1=Database->clist; chn1; chn1=chn1->cnext)
  1518.         if( chn1->glist )
  1519.         {   if( IsProtein(chn1->glist->refno) )
  1520.         {   CalcProteinHBonds(chn1);
  1521.         } else if( IsDNA(chn1->glist->refno) )
  1522.                 CalcNucleicHBonds(chn1);
  1523.         }
  1524.  
  1525.     if( FileDepth == -1 )
  1526.     {   if( CommandActive )
  1527.         WriteChar('\n');
  1528.         CommandActive=False;
  1529.     
  1530.         sprintf(buffer,"Number of H-Bonds ... %d\n",InfoHBondCount);
  1531.         WriteString(buffer);
  1532.     }
  1533. }
  1534.  
  1535.  
  1536. static int IsHBonded( src, dst, ptr )
  1537.     Atom __far *src, __far *dst;
  1538.     HBond __far *ptr;
  1539. {
  1540.     while( ptr && (ptr->srcCA==src) )
  1541.     if( ptr->dstCA == dst )
  1542.     {   return( True );
  1543.     } else ptr=ptr->hnext;
  1544.     return( False );
  1545. }
  1546.  
  1547.  
  1548. static void FindAlphaHelix( pitch, flag )
  1549.     int pitch, flag;
  1550. {
  1551.     register HBond __far *hbond;
  1552.     register Chain __far *chain;
  1553.     register Group __far *group;
  1554.     register Group __far *first;
  1555.     register Group __far *ptr;
  1556.     register Atom __far *srcCA;
  1557.     register Atom __far *dstCA;
  1558.     register int res,dist,prev;
  1559.  
  1560.     /* Protein chains only! */
  1561.     hbond = Database->hlist;
  1562.     for( chain=Database->clist; chain; chain=chain->cnext )
  1563.     if( (first=chain->glist) && IsProtein(first->refno) )
  1564.     {   prev = False; dist = 0;
  1565.     for( group=chain->glist; hbond && group; group=group->gnext )
  1566.     {   if( IsAmino(group->refno) && (srcCA=FindGroupAtom(group,1)) )
  1567.         {   if( dist==pitch )
  1568.         {   res = False;
  1569.             dstCA=FindGroupAtom(first,1);
  1570.  
  1571.             while( hbond && hbond->srcCA == srcCA )
  1572.             {   if( hbond->dstCA==dstCA ) res=True;
  1573.             hbond = hbond->hnext;
  1574.             }
  1575.  
  1576.             if( res )
  1577.             {   if( prev )
  1578.             {   if( !(first->struc&HelixFlag) ) 
  1579.                 InfoHelixCount++;
  1580.  
  1581.                 ptr = first;
  1582.                 do {
  1583.                 ptr->struc |= flag;
  1584.                 ptr = ptr->gnext;
  1585.                 } while( ptr != group );
  1586.             } else prev = True;
  1587.             } else prev = False;
  1588.         } else while( hbond && hbond->srcCA==srcCA )
  1589.             hbond = hbond->hnext;
  1590.         } else prev = False;
  1591.  
  1592.         if( group->struc&HelixFlag )
  1593.         {   first = group; prev = False; dist = 1;
  1594.         } else if( dist==pitch )
  1595.         {   first = first->gnext;
  1596.         } else dist++;
  1597.     }
  1598.     } else if( first && IsNucleo(first->refno) )
  1599.     while( hbond && !IsAminoBackbone(hbond->src->refno) )
  1600.         hbond = hbond->hnext;
  1601. }
  1602.  
  1603.  
  1604. static Atom __far *cprevi, __far *ccurri, __far *cnexti;
  1605. static HBond __far *hcurri, __far *hnexti;
  1606. static Group __far *curri, __far *nexti;
  1607.  
  1608.  
  1609.  
  1610. static void TestLadder( chain )
  1611.     Chain __far *chain;
  1612. {
  1613.     register Atom __far *cprevj, __far *ccurrj, __far *cnextj;
  1614.     register HBond __far *hcurrj, __far *hnextj;
  1615.     register Group __far *currj, __far *nextj;
  1616.     register int count, result, found;
  1617.  
  1618.     /* Already part of atleast one ladder */
  1619.     found = curri->flag & SheetFlag;
  1620.     nextj = nexti->gnext;
  1621.  
  1622.     hnextj = hnexti;
  1623.     while( hnextj && hnextj->srcCA==cnexti )
  1624.     hnextj = hnextj->hnext;
  1625.  
  1626.     while( True )
  1627.     {   if( nextj )
  1628.         if( IsProtein(chain->glist->refno) )
  1629.         {   count = 1;
  1630.         do {
  1631.             cnextj = FindGroupAtom(nextj,1);
  1632.             if( count == 3 )
  1633.             {   if( IsHBonded(cnexti,ccurrj,hnexti) &&
  1634.                 IsHBonded(ccurrj,cprevi,hcurrj) )
  1635.             {   result = ParaLadder;
  1636.             } else if( IsHBonded(cnextj,ccurri,hnextj) &&
  1637.                    IsHBonded(ccurri,cprevj,hcurri) )
  1638.             {   result = ParaLadder;
  1639.             } else if( IsHBonded(cnexti,cprevj,hnexti) &&
  1640.                    IsHBonded(cnextj,cprevi,hnextj) )
  1641.             {   result = AntiLadder;
  1642.             } else if( IsHBonded(ccurrj,ccurri,hcurrj) &&
  1643.                    IsHBonded(ccurri,ccurrj,hcurri) )
  1644.             {   result = AntiLadder;
  1645.             } else result = NoLadder;
  1646.  
  1647.             if( result )
  1648.             {   curri->struc |= SheetFlag;
  1649.                 currj->struc |= SheetFlag;
  1650.                 if( found ) return;
  1651.                 found = True;
  1652.             }
  1653.             } else count++;
  1654.  
  1655.             cprevj = ccurrj; ccurrj = cnextj; 
  1656.             currj = nextj;   hcurrj = hnextj;
  1657.  
  1658.             while( hnextj && hnextj->srcCA==cnextj )
  1659.             hnextj = hnextj->hnext;
  1660.         } while( (nextj = nextj->gnext) );
  1661.  
  1662.         } else if( IsNucleo(chain->glist->refno) )
  1663.         while( hnextj && !IsAminoBackbone(hnextj->src->refno) )
  1664.             hnextj = hnextj->hnext;
  1665.  
  1666.     if( (chain = chain->cnext) ) 
  1667.     {   nextj = chain->glist;
  1668.     } else return;
  1669.     }
  1670. }
  1671.  
  1672.  
  1673. static void FindBetaSheets()
  1674. {
  1675.     register Chain __far *chain;
  1676.     register int ladder;
  1677.     register int count;
  1678.  
  1679.     hnexti = Database->hlist;
  1680.     for( chain=Database->clist; chain; chain=chain->cnext )
  1681.     if( (nexti = chain->glist) )
  1682.         if( IsProtein(nexti->refno) )
  1683.         {   count = 1;
  1684.         ladder = False;
  1685.         do {
  1686.             cnexti = FindGroupAtom(nexti,1);
  1687.  
  1688.             if( count == 3 )
  1689.             {   TestLadder( chain );
  1690.             if( curri->struc & SheetFlag )
  1691.             {   if( !ladder )
  1692.                 {   InfoLadderCount++;
  1693.                 ladder = True;
  1694.                 }
  1695.             } else ladder = False;
  1696.             } else count++;
  1697.  
  1698.             cprevi = ccurri; ccurri = cnexti; 
  1699.             curri = nexti;   hcurri = hnexti;
  1700.             while( hnexti && hnexti->srcCA==cnexti )
  1701.             hnexti = hnexti->hnext;
  1702.         } while( (nexti = nexti->gnext) );
  1703.  
  1704.         } else if( IsNucleo(nexti->refno) )
  1705.         while( hnexti && !IsAminoBackbone(hnexti->src->refno) )
  1706.             hnexti = hnexti->hnext;
  1707. }
  1708.  
  1709.  
  1710. static void FindTurnStructure()
  1711. {
  1712.     static Atom __far *aptr[5];
  1713.     register Chain __far *chain;
  1714.     register Group __far *group;
  1715.     register Group __far *prev;
  1716.     register Atom __far *ptr;
  1717.     register Long ux,uy,uz,mu;
  1718.     register Long vx,vy,vz,mv;
  1719.     register int i,found,len;
  1720.     register Real CosKappa;
  1721.  
  1722.     for( chain=Database->clist; chain; chain=chain->cnext )
  1723.     if( chain->glist && IsProtein(chain->glist->refno) )
  1724.     {   len = 0;  found = False;
  1725.         for( group=chain->glist; group; group=group->gnext )
  1726.         {    ptr = FindGroupAtom(group,1);
  1727.          if( ptr && (ptr->flag&BreakFlag) )
  1728.          {   found = False;
  1729.              len = 0;
  1730.          } else if( len==5 )
  1731.          {   for( i=0; i<4; i++ )
  1732.              aptr[i] = aptr[i+1];
  1733.              len = 4;
  1734.          } else if( len==2 )
  1735.              prev = group;
  1736.  
  1737.          aptr[len++] = ptr;
  1738.          if( len==5 ) 
  1739.          {   if( !(prev->struc&(HelixFlag|SheetFlag)) &&
  1740.              aptr[0] && aptr[2] && aptr[4] )
  1741.              {   ux = aptr[2]->xorg - aptr[0]->xorg;
  1742.              uy = aptr[2]->yorg - aptr[0]->yorg;
  1743.              uz = aptr[2]->zorg - aptr[0]->zorg;
  1744.  
  1745.              vx = aptr[4]->xorg - aptr[2]->xorg;
  1746.              vy = aptr[4]->yorg - aptr[2]->yorg;
  1747.              vz = aptr[4]->zorg - aptr[2]->zorg;
  1748.  
  1749.              mu = ux*ux + uy*uy + uz*uz;
  1750.              mv = vx*vx + vz*vz + vy*vy;
  1751.              if( mu && mv )
  1752.              {   CosKappa = (Real)(ux*vx + uy*vy + uz*vz);
  1753.                  CosKappa /= sqrt( (Real)mu*mv );
  1754.                  if( CosKappa<Cos70Deg )
  1755.                  {   if( !found )
  1756.                      InfoTurnCount++;
  1757.                  prev->struc |= TurnFlag;
  1758.                  }
  1759.              }
  1760.              }
  1761.              found = prev->struc&TurnFlag;
  1762.              prev = prev->gnext;
  1763.          } /* len==5 */
  1764.         }
  1765.     }
  1766. }
  1767.  
  1768.  
  1769. static void FindBetaTurns()
  1770. {
  1771.     static Atom __far *aptr[4];
  1772.     register Chain __far *chain;
  1773.     register Group __far *group;
  1774.     register Group __far *prev;
  1775.     register Group __far *next;
  1776.     register Atom __far *ptr;
  1777.     register Long dx,dy,dz;
  1778.     register int found,len;
  1779.     register int flag;
  1780.  
  1781.     for( chain=Database->clist; chain; chain=chain->cnext )
  1782.         if( chain->glist && IsProtein(chain->glist->refno) )
  1783.         {   prev = chain->glist;  
  1784.             len = 0;  found = False;
  1785.             for( next=chain->glist; next; next=next->gnext )
  1786.             {   ptr = FindGroupAtom(next,1);
  1787.                 if( ptr && (ptr->flag&BreakFlag) )
  1788.                 {   found = False;
  1789.                     prev = next;
  1790.                     len = 0;
  1791.                 } else if( len==4 )
  1792.                 {   aptr[0] = aptr[1];
  1793.                     aptr[1] = aptr[2];
  1794.                     aptr[2] = aptr[3];
  1795.                     aptr[3] = ptr;
  1796.  
  1797.                 } else aptr[len++] = ptr;
  1798.                 if( len==4 ) 
  1799.                 {   flag = False;
  1800.                     if( aptr[0] && aptr[3] )
  1801.                     {   dx = aptr[3]->xorg - aptr[0]->xorg;
  1802.                         dy = aptr[3]->yorg - aptr[0]->yorg;
  1803.                         dz = aptr[3]->zorg - aptr[0]->zorg;
  1804.                         if( dx*dx + dy*dy + dz*dz < (Long)1750*1750 )
  1805.                         {   group = prev;
  1806.                             while( group!=next->gnext )
  1807.                             {   if( !(group->struc&(HelixFlag|SheetFlag)) )
  1808.                                 {   group->struc |= TurnFlag;
  1809.                                     flag = True;
  1810.                                 }
  1811.                                 group = group->gnext;
  1812.                             }
  1813.                             if( !found && flag ) 
  1814.                                 InfoTurnCount++;
  1815.                         }
  1816.                     }
  1817.                     prev = prev->gnext;   
  1818.                     found = flag;
  1819.                 } /* len==4 */
  1820.             }
  1821.         }
  1822. }
  1823.  
  1824.  
  1825. void DetermineStructure( flag )
  1826.     int flag;
  1827. {
  1828.     register Chain __far *chain;
  1829.     register Group __far *group;
  1830.     char buffer[40];
  1831.  
  1832.     if( !Database )
  1833.     return;
  1834.  
  1835.     if( InfoHBondCount<0 )
  1836.     CalcHydrogenBonds();
  1837.  
  1838.     if( InfoHelixCount>=0 )
  1839.     for( chain=Database->clist; chain; chain=chain->cnext )
  1840.         for( group=chain->glist; group; group=group->gnext )
  1841.         group->struc = 0;
  1842.  
  1843.     InfoStrucSrc = SourceCalc;
  1844.     InfoLadderCount = 0;
  1845.     InfoHelixCount = 0;
  1846.     InfoTurnCount = 0;
  1847.  
  1848.     if( InfoHBondCount )
  1849.     {   FindAlphaHelix(4,Helix4Flag);
  1850.     FindBetaSheets();
  1851.     FindAlphaHelix(3,Helix3Flag);
  1852.     FindAlphaHelix(5,Helix5Flag);
  1853.  
  1854.         if( !flag )
  1855.     {   FindTurnStructure();
  1856.         } else FindBetaTurns();
  1857.     }
  1858.  
  1859.     if( FileDepth == -1 )
  1860.     {   if( CommandActive )
  1861.         WriteChar('\n');
  1862.         CommandActive=False;
  1863.  
  1864.         sprintf(buffer,"Number of Helices ... %d\n",InfoHelixCount);
  1865.         WriteString(buffer);
  1866.         sprintf(buffer,"Number of Strands ... %d\n",InfoLadderCount);
  1867.         WriteString(buffer);
  1868.         sprintf(buffer,"Number of Turns ..... %d\n",InfoTurnCount);
  1869.         WriteString(buffer);
  1870.     }
  1871. }
  1872.  
  1873.  
  1874. void RenumberMolecule( start )
  1875.     int start;
  1876. {
  1877.     register Chain __far *chain;
  1878.     register Group __far *group;
  1879.     register int hinit, minit;
  1880.     register int resno;
  1881.  
  1882.     if( !Database )
  1883.     return;
  1884.  
  1885.     hinit = minit = False;
  1886.     for( chain=Database->clist; chain; chain=chain->cnext )
  1887.     {   resno = start;
  1888.     for( group=chain->glist; group; group=group->gnext )
  1889.     {   if( group->alist->flag & HeteroFlag )
  1890.         {   if( hinit )
  1891.         {   if( resno > MaxHetaRes )
  1892.             {   MaxHetaRes = resno;
  1893.             } else if( resno < MinHetaRes )
  1894.             MinHetaRes = resno;
  1895.         } else MinHetaRes = MaxHetaRes = resno;
  1896.         hinit = True;
  1897.         } else
  1898.         {   if( minit )
  1899.         {   if( resno > MaxMainRes )
  1900.             {   MaxMainRes = resno;
  1901.             } else if( resno < MinMainRes )
  1902.             MinMainRes = resno;
  1903.         } else MinMainRes = MaxMainRes = resno;
  1904.         minit = True;
  1905.         }
  1906.         group->serno = resno++;
  1907.     }
  1908.     }
  1909. }
  1910.  
  1911.  
  1912. /*===============================*/
  1913. /* Molecule Database Maintenance */
  1914. /*===============================*/
  1915.  
  1916. static void ReclaimAtoms( ptr )
  1917.     Atom __far *ptr;
  1918. {
  1919.     register Atom __far *temp;
  1920.  
  1921.     if( (temp = ptr) )
  1922.     {   while( temp->anext )
  1923.         temp=temp->anext;
  1924.     temp->anext = FreeAtom;
  1925.     FreeAtom = ptr;
  1926.     }
  1927. }
  1928.  
  1929. static void ResetDatabase()
  1930. {
  1931.     Database = CurMolecule = (void __far*)0;
  1932.     MainGroupCount = HetaGroupCount = 0;
  1933.     InfoChainCount = HetaAtomCount = 0;
  1934.     MainAtomCount = InfoBondCount = 0;  
  1935.     SelectCount = 0;
  1936.  
  1937.     InfoStrucSrc = SourceNone;
  1938.     InfoSSBondCount = InfoHBondCount = -1;
  1939.     InfoHelixCount = InfoLadderCount = -1;
  1940.     InfoTurnCount = -1;
  1941.  
  1942.     CurGroup = (void __far*)0;
  1943.     CurChain = (void __far*)0;
  1944.     CurAtom = (void __far*)0;
  1945.  
  1946.     MinX = MinY = MinZ = 0;
  1947.     MaxX = MaxY = MaxZ = 0;
  1948.  
  1949.     MinMainTemp = MaxMainTemp = 0;
  1950.     MinHetaTemp = MaxHetaTemp = 0;
  1951.     MinMainRes = MaxMainRes = 0;
  1952.     MinHetaRes = MaxHetaRes = 0;
  1953.  
  1954.     *InfoMoleculeName = 0;
  1955.     *InfoClassification = 0;
  1956.     *InfoIdentCode = 0;
  1957.     *InfoSpaceGroup = 0;
  1958.     *InfoFileName = 0;
  1959.  
  1960.     VoxelsClean = False;
  1961.     HMinMaxFlag = False;
  1962.     MMinMaxFlag = False;
  1963.     HasHydrogen = False;
  1964.     ElemNo = MINELEM;
  1965.     ResNo = MINRES;
  1966.     MaskCount = 0;
  1967.     NMRModel = 0;
  1968. }
  1969.  
  1970.  
  1971. void DestroyDatabase()
  1972. {
  1973.     register void __far *temp;
  1974.     register Group __far *gptr;
  1975.  
  1976.     if( Database )
  1977.     {   ReclaimHBonds( Database->slist );
  1978.     ReclaimHBonds( Database->hlist );
  1979.     ReclaimBonds( Database->blist );
  1980.  
  1981.     while( Database->clist )
  1982.     {   ReclaimBonds(Database->clist->blist);
  1983.         if( (gptr = Database->clist->glist) )
  1984.         {   ReclaimAtoms(gptr->alist);
  1985.         while( gptr->gnext )
  1986.         {   gptr = gptr->gnext;
  1987.             ReclaimAtoms(gptr->alist);
  1988.         }
  1989.         gptr->gnext = FreeGroup;
  1990.         FreeGroup = Database->clist->glist;
  1991.         }
  1992.         temp = Database->clist->cnext;
  1993.         Database->clist->cnext = FreeChain;
  1994.         FreeChain = Database->clist;
  1995.         Database->clist = temp;
  1996.     }
  1997.  
  1998.     FreeMolecule = Database;
  1999.     Database = (void __far*)0;
  2000.     }
  2001.     ResetDatabase();
  2002. }
  2003.  
  2004.  
  2005. void PurgeDatabase()
  2006. {
  2007. #ifdef APPLEMAC
  2008.     register AllocRef *ptr;
  2009.     register AllocRef *tmp;
  2010.     register int i;
  2011.     
  2012.     /* Avoid Memory Leaks */
  2013.     for( ptr=AllocList; ptr; ptr=tmp )
  2014.     {   for( i=0; i<ptr->count; i++ )
  2015.         _ffree( ptr->data[i] );
  2016.     tmp = ptr->next;
  2017.     _ffree( ptr );
  2018.     }
  2019. #endif
  2020. }
  2021.  
  2022.  
  2023. void InitialiseDatabase()
  2024. {
  2025.     FreeMolecule = (void __far*)0;
  2026.     FreeHBond = (void __far*)0;
  2027.     FreeChain = (void __far*)0;
  2028.     FreeGroup = (void __far*)0;
  2029.     FreeAtom = (void __far*)0;
  2030.     FreeBond = (void __far*)0;
  2031.  
  2032. #ifdef APPLEMAC
  2033.     AllocList = (void*)0;
  2034. #endif
  2035.  
  2036.     ResetDatabase();
  2037. }
  2038.  
  2039.