home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 March / macformat-022.iso / Shareware City / Science / RasMol2 / molecule.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-28  |  63.8 KB  |  2,759 lines  |  [TEXT/KAHL]
  1. /* molecule.c
  2.  * RasMol2 Molecular Graphics
  3.  * Roger Sayle, October 1994
  4.  * Version 2.5
  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.  
  32. #define GroupPool    8
  33. #define HBondPool   32
  34. #define BondPool    32
  35. #define AtomPool    32
  36.  
  37. #define NoLadder     0x00
  38. #define ParaLadder   0x01
  39. #define AntiLadder   0x02
  40.  
  41. #define Cos70Deg     0.34202014332567
  42.  
  43. #define FeatHelix    1
  44. #define FeatSheet    2
  45. #define FeatTurn     3
  46.  
  47. #define SourceNone   0
  48. #define SourcePDB    1
  49. #define SourceCalc   2
  50.  
  51. #define MaxHBondDist   ((Long)300*300)
  52. #define MaxBondDist    ((Long)475*475)
  53. #define MinBondDist    ((Long)100*100)
  54. #define AbsMaxBondDist 500
  55.  
  56. typedef struct {
  57.     int init, term;
  58.     char chain;
  59.     char type;
  60.     } FeatEntry;
  61.  
  62. #ifdef APPLEMAC
  63. #define AllocSize   256
  64. typedef struct _AllocRef {
  65.     struct _AllocRef *next;
  66.     void *data[AllocSize];
  67.     int count;
  68.     } AllocRef;
  69. static AllocRef *AllocList;  
  70. #endif
  71.  
  72.  
  73. #define FeatSize    32
  74. typedef struct _Feature {
  75.         struct _Feature __far *fnext;
  76.     FeatEntry data[FeatSize];
  77.     int count;
  78.     } Feature;
  79.  
  80. typedef struct {
  81.       char src[4];
  82.       char dst[4];
  83.       } ConvTable;
  84.  
  85. #define MAXALCATOM   5
  86. static ConvTable AlcAtomTable[MAXALCATOM] = {
  87.     { { 'S', 'O', '2', ' ' }, { ' ', 'S', '2', ' ' } },  /*  1 */
  88.     { { 'C', 'A', 'R', ' ' }, { ' ', 'C', ' ', ' ' } },  /*  2 */
  89.     { { 'N', 'A', 'R', ' ' }, { ' ', 'N', ' ', ' ' } },  /*  3 */
  90.     { { 'N', 'A', 'M', ' ' }, { ' ', 'N', ' ', ' ' } },  /*  4 */
  91.     { { 'N', 'P', 'L', '3' }, { ' ', 'N', '3', ' ' } },  /*  5 */
  92.                  };
  93.  
  94. static Molecule __far *FreeMolecule;
  95. static HBond __far *FreeHBond;
  96. static Chain __far *FreeChain;
  97. static Group __far *FreeGroup;
  98. static Atom __far *FreeAtom;
  99. static Bond __far *FreeBond;
  100.  
  101. static char PDBInsert;
  102. static Feature __far *FeatList;
  103. static HBond __far * __far *CurHBond;
  104. static Molecule __far *CurMolecule;
  105. static Chain __far *CurChain;
  106. static Group __far *CurGroup;
  107. static Atom __far *CurAtom;
  108. static Atom __far *ConnectAtom;
  109. static int InfoStrucSource;
  110. static int HMinMaxFlag;
  111. static int MMinMaxFlag;
  112. static int MemSize;
  113.  
  114. static char Record[82];
  115. static FILE *DataFile;
  116.  
  117. /* Macros for commonly used loops */
  118. #define ForEachAtom  for(chain=Database->clist;chain;chain=chain->cnext) \
  119.              for(group=chain->glist;group;group=group->gnext)    \
  120.              for(aptr=group->alist;aptr;aptr=aptr->anext)
  121. #define ForEachBond  for(bptr=Database->blist;bptr;bptr=bptr->bnext)
  122.  
  123.  
  124. #ifdef APPLEMAC
  125. /* External RasMac Function Declaration! */
  126. void SetFileInfo( char*, OSType, OSType, short );
  127. #endif
  128.  
  129.  
  130. static void FatalDataError(ptr)
  131.     char *ptr;
  132. {
  133.     char buffer[80];
  134.  
  135.     sprintf(buffer,"Database Error: %s!",ptr);
  136.     RasMolFatalExit(buffer);
  137. }
  138.  
  139.  
  140. static void FetchRecord()
  141. {
  142.     register char *ptr;
  143.     register int ch;
  144.  
  145.     ptr = Record + 1;
  146.     if( !feof(DataFile) )
  147.     {   do {
  148.         ch = getc(DataFile);
  149.         if( (ch=='\n') || (ch==EOF) )
  150.         {   *ptr = 0;
  151.         return;
  152.             } else if( (ch=='\r') )
  153.             {   ch = getc(DataFile);
  154.                 if( ch != '\n' )
  155.                     ungetc(ch,DataFile);
  156.                 *ptr = 0;
  157.                 return;
  158.         } else *ptr++ = ch;
  159.     } while( ptr < Record+80 );
  160.  
  161.     /* skip to the end of the line! */
  162.     do { ch = getc(DataFile);
  163.     } while( (ch!='\n') && (ch!='\r') && (ch!=EOF) );
  164.  
  165.         if( ch=='\r' )
  166.         {   ch = getc(DataFile);
  167.             if( ch != '\n' )
  168.                 ungetc(ch,DataFile);
  169.         }
  170.     }
  171.     *ptr = 0;
  172. }
  173.  
  174.  
  175. static void ExtractString( len, src, dst )
  176.     int len;  char *src, *dst;
  177. {
  178.     register char *ptr;
  179.     register char ch;
  180.     register int i;
  181.  
  182.     ptr = dst;
  183.     for( i=0; i<len; i++ )
  184.     {   if( (ch = *src++) )
  185.     {   if( (*dst++ = ch) != ' ' ) 
  186.         ptr = dst;
  187.     } else break;
  188.     }
  189.     *ptr = 0;
  190. }
  191.  
  192.  
  193. static Long ReadValue( pos, len )
  194.     int pos, len;
  195. {
  196.     register Long result;
  197.     register char *ptr;
  198.     register char ch;
  199.     register int neg;
  200.  
  201.     result = 0;
  202.     neg = False;
  203.     ptr = Record+pos;
  204.     while( len-- )
  205.     {   ch = *ptr++;
  206.     if( (ch>='0') && (ch<='9') )
  207.     {   result = (10*result)+(ch-'0');
  208.     } else if( ch=='-' )
  209.         neg = True;
  210.     }
  211.     return( neg? -result : result );
  212. }
  213.  
  214.  
  215. void DescribeMolecule()
  216. {
  217.     char buffer[40];
  218.  
  219.     if( CommandActive )
  220.     WriteChar('\n');
  221.     CommandActive=False;
  222.  
  223.     if( *InfoMoleculeName )
  224.     {   WriteString("Molecule name ....... ");
  225.     WriteString(InfoMoleculeName);
  226.     WriteChar('\n');
  227.     }
  228.  
  229.     if( *InfoClassification )
  230.     {   WriteString("Classification ...... ");
  231.     WriteString(InfoClassification);
  232.     WriteChar('\n');
  233.     }
  234.  
  235.     if( Database && (MainGroupCount>1) )
  236.     {   WriteString("Secondary Structure . ");
  237.     if( InfoStrucSource==SourceNone )
  238.     {   WriteString("No Assignment\n");
  239.     } else if( InfoStrucSource==SourcePDB )
  240.     {   WriteString("PDB Data Records\n");
  241.     } else WriteString("Calculated\n");
  242.     }
  243.  
  244.  
  245.     if( *InfoIdentCode )
  246.     {   WriteString("Brookhaven Code ..... ");
  247.     WriteString(InfoIdentCode);
  248.     WriteChar('\n');
  249.     }
  250.  
  251.     if( InfoChainCount>1 )
  252.     {   sprintf(buffer,"Number of Chains .... %d\n",InfoChainCount);
  253.     WriteString(buffer);
  254.     }
  255.  
  256.     sprintf(buffer,"Number of Groups .... %d",MainGroupCount);
  257.     WriteString(buffer);
  258.     if( HetaAtomCount )
  259.     {   sprintf(buffer," (%d)\n",HetaGroupCount);
  260.     WriteString(buffer);
  261.     } else WriteChar('\n');
  262.  
  263.     sprintf(buffer,"Number of Atoms ..... %ld",MainAtomCount);
  264.     WriteString(buffer);
  265.     if( HetaAtomCount )
  266.     {   sprintf(buffer," (%d)\n",HetaAtomCount);
  267.     WriteString(buffer);
  268.     } else WriteChar('\n');
  269.  
  270.     if( InfoBondCount )
  271.     {   sprintf(buffer,"Number of Bonds ..... %ld\n",InfoBondCount);
  272.     WriteString(buffer);
  273.     }
  274.  
  275.     if( InfoSSBondCount != -1 )
  276.     {   WriteString("Number of Bridges ... ");
  277.     sprintf(buffer,"%d\n\n",InfoSSBondCount);
  278.     WriteString(buffer);
  279.     }
  280.  
  281.     if( InfoHBondCount != -1 )
  282.     {   WriteString("Number of H-Bonds ... ");
  283.     sprintf(buffer,"%d\n",InfoHBondCount);
  284.     WriteString(buffer);
  285.     }
  286.  
  287.     if( InfoHelixCount != -1 )
  288.     {   WriteString("Number of Helices ... ");
  289.     sprintf(buffer,"%d\n",InfoHelixCount);
  290.     WriteString(buffer);
  291.  
  292.     WriteString("Number of Strands ... ");
  293.     sprintf(buffer,"%d\n",InfoLadderCount);
  294.     WriteString(buffer);
  295.  
  296.     WriteString("Number of Turns ..... ");
  297.     sprintf(buffer,"%d\n",InfoTurnCount);
  298.     WriteString(buffer);
  299.     }
  300. }
  301.  
  302.  
  303. #ifdef APPLEMAC
  304. /* Avoid System Memory Leaks! */
  305. static void RegisterAlloc( data )
  306.     void *data;
  307. {
  308.     register AllocRef *ptr;
  309.     
  310.     if( !AllocList || (AllocList->count==AllocSize) )
  311.     {   ptr = (AllocRef *)_fmalloc( sizeof(AllocRef) );
  312.     if( !ptr ) FatalDataError("Memory allocation failed");
  313.     
  314.     ptr->next = AllocList;
  315.     ptr->data[0] = data;
  316.     ptr->count = 1;
  317.     AllocList = ptr;
  318.     } else AllocList->data[AllocList->count++] = data;
  319. }
  320. #else
  321. #define RegisterAlloc(x)
  322. #endif
  323.  
  324.  
  325. static void CreateChain( ident )
  326.     char ident;
  327. {
  328.     register Chain __far *prev;
  329.  
  330.     if( !CurMolecule )
  331.     {   if( !(CurMolecule = FreeMolecule) )
  332.     {   MemSize += sizeof(Molecule);
  333.         CurMolecule = (Molecule __far *)_fmalloc(sizeof(Molecule));
  334.         if( !CurMolecule ) FatalDataError("Memory allocation failed");
  335.         RegisterAlloc( CurMolecule );
  336.     } else FreeMolecule = (void __far*)0;
  337.  
  338.     CurChain = (void __far*)0;
  339.     CurMolecule->slist = (void __far*)0;
  340.     CurMolecule->hlist = (void __far*)0;
  341.     CurMolecule->blist = (void __far*)0;
  342.     CurMolecule->clist = (void __far*)0;
  343.     Database = CurMolecule;
  344.     }
  345.  
  346.     /* Handle chain breaks! */
  347.     if( !(prev=CurChain) )
  348.     if( (prev=CurMolecule->clist) )
  349.         while( prev->cnext )
  350.         prev = prev->cnext;
  351.  
  352.     if( !(CurChain = FreeChain) )
  353.     {   MemSize += sizeof(Chain);
  354.     CurChain = (Chain __far *)_fmalloc(sizeof(Chain));
  355.     if( !CurChain ) FatalDataError("Memory allocation failed");
  356.     RegisterAlloc( CurChain );
  357.     } else FreeChain = FreeChain->cnext;
  358.  
  359.     if( prev )
  360.     {   prev->cnext = CurChain;
  361.     } else CurMolecule->clist = CurChain;
  362.     CurChain->cnext = (void __far*)0;
  363.      
  364.     CurChain->ident = ident;
  365.     CurChain->glist = (void __far*)0;
  366.     CurChain->blist = (void __far*)0;
  367.     ConnectAtom = (void __far*)0;
  368.     CurGroup = (void __far*)0;
  369.     InfoChainCount++;
  370. }
  371.  
  372.  
  373. static void CreateGroup( pool )
  374.     int pool;
  375. {
  376.     register Group __far *ptr;
  377.     register int i;
  378.  
  379.     if( !(ptr = FreeGroup) )
  380.     {   MemSize += pool*sizeof(Group);
  381.     ptr = (Group __far *)_fmalloc( pool*sizeof(Group) );
  382.     if( !ptr ) FatalDataError("Memory allocation failed");
  383.     RegisterAlloc( ptr );
  384.     for( i=1; i<pool; i++ )
  385.     {   ptr->gnext = FreeGroup;
  386.         FreeGroup = ptr++;
  387.     } 
  388.     } else FreeGroup = ptr->gnext;
  389.     
  390.     if( CurGroup )
  391.     {   CurGroup->gnext = ptr;
  392.     } else CurChain->glist = ptr;
  393.     CurGroup = ptr;
  394.  
  395.     CurAtom = (void __far*)0;
  396.     ptr->gnext = (void __far*)0;
  397.     ptr->alist = (void __far*)0;
  398.     ptr->struc = 0;
  399.     ptr->flag = 0;
  400.     ptr->col1 = 0;
  401.     ptr->col2 = 0;
  402. }
  403.  
  404.  
  405. static Atom __far *CreateAtom()
  406. {
  407.     register Atom __far *ptr;
  408.     register int i;
  409.  
  410.     if( !(ptr = FreeAtom) )
  411.     {   MemSize += AtomPool*sizeof(Atom);
  412.     ptr = (Atom __far *)_fmalloc( AtomPool*sizeof(Atom) );
  413.     if( !ptr ) FatalDataError("Memory allocation failed");
  414.     RegisterAlloc( ptr );
  415.     for( i=1; i<AtomPool; i++ )
  416.     {   ptr->anext = FreeAtom;
  417.         FreeAtom = ptr++;
  418.     } 
  419.     } else FreeAtom = ptr->anext;
  420.  
  421.     if( CurAtom )
  422.     {   CurAtom->anext = ptr;
  423.     } else CurGroup->alist = ptr;
  424.     ptr->anext = (void __far*)0;
  425.     CurAtom = ptr;
  426.  
  427.     SelectCount++;
  428.     ptr->flag = SelectFlag | NonBondFlag;
  429.     ptr->label = (void*)0;
  430.     ptr->radius = 375;
  431.     ptr->altl = ' ';
  432.     ptr->mbox = 0;
  433.     ptr->col = 0;
  434.  
  435.     return( ptr );
  436. }
  437.  
  438.  
  439. static void ProcessAtom( ptr )
  440.     Atom __far *ptr;
  441. {
  442.     if( GetElemNumber(ptr) == 1 )
  443.     {   ptr->flag |= HydrogenFlag;
  444.     HasHydrogen = True;
  445.     }
  446.  
  447.     if( !(ptr->flag&(HydrogenFlag|HeteroFlag)) )
  448.     ptr->flag |= NormAtomFlag;
  449.  
  450. #ifdef INVERT
  451.     ptr->yorg = -ptr->yorg;
  452. #endif
  453.  
  454.     if( HMinMaxFlag || MMinMaxFlag )
  455.     {   if( ptr->xorg < MinX ) 
  456.     {   MinX = ptr->xorg;
  457.     } else if( ptr->xorg > MaxX ) 
  458.         MaxX = ptr->xorg;
  459.  
  460.     if( ptr->yorg < MinY ) 
  461.     {   MinY = ptr->yorg;
  462.     } else if( ptr->yorg > MaxY ) 
  463.         MaxY = ptr->yorg;
  464.  
  465.     if( ptr->zorg < MinZ ) 
  466.     {   MinZ = ptr->zorg;
  467.     } else if( ptr->zorg > MaxZ ) 
  468.         MaxZ = ptr->zorg;
  469.     } else 
  470.     {   MinX = MaxX = ptr->xorg;
  471.     MinY = MaxY = ptr->yorg;
  472.     MinZ = MaxZ = ptr->zorg;
  473.     }
  474.         
  475.     if( ptr->flag & HeteroFlag )
  476.     {   if( HMinMaxFlag )
  477.     {   if( ptr->temp < MinHetaTemp ) 
  478.         {   MinHetaTemp = ptr->temp;
  479.         } else if( ptr->temp > MaxHetaTemp ) 
  480.         MaxHetaTemp = ptr->temp;
  481.     } else MinHetaTemp = MaxHetaTemp = ptr->temp;
  482.     HMinMaxFlag = True;
  483.     HetaAtomCount++;
  484.     } else
  485.     {   if( MMinMaxFlag )
  486.     {   if( ptr->temp < MinMainTemp ) 
  487.         {   MinMainTemp = ptr->temp;
  488.         } else if( ptr->temp > MaxMainTemp ) 
  489.         MaxMainTemp = ptr->temp;
  490.     } else MinMainTemp = MaxMainTemp = ptr->temp;
  491.     MMinMaxFlag = True;
  492.     MainAtomCount++;
  493.     }
  494. }
  495.  
  496.  
  497. static int FindResNo( ptr )
  498.     char *ptr;
  499. {
  500.     register int refno;
  501.  
  502.     for( refno=0; refno<ResNo; refno++ )
  503.     if( !strncmp(Residue[refno],ptr,3) )
  504.         return( refno );
  505.  
  506.     if( !strncmp(ptr,"CSH",3) ||
  507.     !strncmp(ptr,"CYH",3) ||
  508.     !strncmp(ptr,"CSM",3) )
  509.     {   return(17);  /* cystine */
  510.     } else if( !strncmp(ptr,"WAT",3) ||
  511.            !strncmp(ptr,"H2O",3) ||
  512.            !strncmp(ptr,"SOL",3) ||
  513.            !strncmp(ptr,"TIP",3) )
  514.     /* check HHO, OHH, 0H2? */
  515.     {   return(31);  /* Water (HOH) */
  516.     } else if( !strncmp(ptr,"D2O",3) )
  517.     {   return(32);  /* Heavy water (DOD) */
  518.     } else if( !strncmp(ptr,"SUL",3) )
  519.     {   return(33);  /* Sulphate (SO4) */
  520.     } else if( !strncmp(ptr,"CPR",3) )
  521.     {   return(11);  /* cis-proline */
  522.     } else if( !strncmp(ptr,"TRY",3) )
  523.     return(15);  /* tryptophan */
  524.     
  525.  
  526.     if( ResNo++ == MAXRES )
  527.     FatalDataError("Too many new residues");
  528.     Residue[refno][0] = *ptr++;
  529.     Residue[refno][1] = *ptr++;
  530.     Residue[refno][2] = *ptr;
  531.     return( refno );
  532. }
  533.  
  534.  
  535. static Long ReadPDBCoord( offset )
  536.     int offset;
  537. {
  538.     register int len,neg;
  539.     register Long result;
  540.     register char *ptr;
  541.     register char ch;
  542.  
  543.     result = 0;
  544.     neg = False;
  545.     len = 8;
  546.  
  547.     ptr = Record+offset;
  548.     while( len-- )
  549.     {   ch = *ptr++;
  550.     if( (ch>='0') && (ch<='9') )
  551.     {   result = (10*result)+(ch-'0');
  552.     } else if( ch=='-' )
  553.         neg = True;
  554.     }
  555.  
  556.     /* Handle Chem3D PDB Files! */
  557.     if( Record[offset+3]=='.' )
  558.     result /= 10;
  559.     return( neg? -result : result );
  560. }
  561.  
  562.  
  563. static void ProcessPDBColourMask()
  564. {
  565.     register MaskDesc *ptr;
  566.     register char *mask;
  567.     register int i;
  568.  
  569.     if( MaskCount==MAXMASK )
  570.     FatalDataError("Too many COLOR records in file");
  571.     ptr = &UserMask[MaskCount];
  572.     mask = ptr->mask;
  573.  
  574.  
  575.     ptr->flags = 0;
  576.     for( i=7; i<12; i++ )
  577.     if( (*mask++ = Record[i]) != '#' )
  578.         ptr->flags |= SerNoFlag;
  579.  
  580.     for( i=13; i<21; i++ )
  581.     *mask++ = Record[i];
  582.     *mask++ = Record[22];
  583.  
  584.     for( i=23; i<27; i++ )
  585.     if( (*mask++ = Record[i]) != '#' )
  586.         ptr->flags |= ResNoFlag;
  587.     *mask++ = Record[27];
  588.  
  589.     ptr->r = (int)(ReadPDBCoord(31)>>2) + 5;
  590.     ptr->g = (int)(ReadPDBCoord(39)>>2) + 5;
  591.     ptr->b = (int)(ReadPDBCoord(47)>>2) + 5;
  592.     ptr->radius = (short)(5*ReadValue(55,6))>>1;
  593.     MaskCount++;
  594. }
  595.  
  596.  
  597. static Bond __far *ProcessBond( src, dst, flag )
  598.     Atom __far *src, __far *dst;
  599.     int flag;
  600. {
  601.     register Bond __far *ptr;
  602.     register int i;
  603.  
  604.     if( !(ptr = FreeBond) )
  605.     {   MemSize += BondPool*sizeof(Bond);
  606.     ptr = (Bond __far *)_fmalloc( BondPool*sizeof(Bond) );
  607.     if( !ptr ) FatalDataError("Memory allocation failed");
  608.     RegisterAlloc( ptr );
  609.     for( i=1; i<BondPool; i++ )
  610.     {   ptr->bnext = FreeBond;
  611.         FreeBond = ptr++;
  612.     } 
  613.     } else FreeBond = ptr->bnext;
  614.     
  615.     ptr->flag = flag | SelectFlag;
  616.     ptr->srcatom = src;
  617.     ptr->dstatom = dst;
  618.     ptr->radius = 0;
  619.     ptr->col = 0;
  620.  
  621.     return( ptr );
  622. }
  623.  
  624.  
  625. static void ConnectAtoms( src, dst, flag )
  626.     int src, dst, flag;
  627. {
  628.     register Chain __far *chain;
  629.     register Group __far *group;
  630.     register Atom __far *aptr;
  631.     register Atom __far *sptr;
  632.     register Atom __far *dptr;
  633.     register Bond __far *bptr;
  634.     register int done;
  635.  
  636.     if( src == dst )
  637.     return;
  638.  
  639.     sptr = (void __far*)0;
  640.     dptr = (void __far*)0;
  641.  
  642.     done = False;
  643.     for( chain=Database->clist; chain && !done; chain=chain->cnext )
  644.     for( group=chain->glist; group && !done; group=group->gnext )
  645.         for( aptr=group->alist; aptr; aptr=aptr->anext )
  646.         {   if( aptr->serno == src )
  647.         {   sptr = aptr;
  648.             if( dptr )
  649.             {   done = True;
  650.             break;
  651.             }
  652.         } else if( aptr->serno == dst )
  653.         {   dptr = aptr;
  654.             if( sptr )
  655.             {   done = True;
  656.             break;
  657.             }
  658.         }
  659.         }
  660.  
  661.  
  662.     /* Both found! */
  663.     if( done ) 
  664.     {   /* Reset Non-bonded flags! */
  665.     sptr->flag &= ~NonBondFlag;
  666.     dptr->flag &= ~NonBondFlag;
  667.  
  668.     bptr = ProcessBond( sptr, dptr, flag );
  669.     bptr->bnext = CurMolecule->blist;
  670.     CurMolecule->blist = bptr;
  671.     InfoBondCount++;
  672.     }
  673. }
  674.  
  675.  
  676. static void PDBConnectAtoms( src, dst )
  677.     int src, dst;
  678. {
  679.     register Bond __far *bptr;
  680.     register int bs,bd;
  681.  
  682.     ForEachBond
  683.     {   bs = bptr->srcatom->serno;
  684.     bd = bptr->dstatom->serno;
  685.  
  686.     if( ((bs==src)&&(bd==dst)) || ((bs==dst)&&(bd==src)) )
  687.     {   if( bptr->flag & NormBondFlag )
  688.         {  /* Convert Single to Double */
  689.            bptr->flag &= ~(NormBondFlag);
  690.            bptr->flag |= DoubBondFlag;
  691.         } else if( bptr->flag & DoubBondFlag )
  692.         {  /* Convert Double to Triple */
  693.            bptr->flag &= ~(DoubBondFlag);
  694.            bptr->flag |= TripBondFlag;
  695.         }
  696.         return;
  697.     }
  698.     }
  699.  
  700.     ConnectAtoms( src, dst, NormBondFlag );
  701. }
  702.  
  703.  
  704.  
  705. static void ProcessPDBGroup( heta, serno )
  706.     int heta, serno;
  707. {
  708.     PDBInsert = Record[27];
  709.     if( !CurChain || (CurChain->ident!=Record[22]) )
  710.     CreateChain( Record[22] );
  711.     CreateGroup( GroupPool );
  712.  
  713.     CurGroup->refno = FindResNo( &Record[18] );
  714.     CurGroup->serno = serno;
  715.  
  716.     /* Solvents should be hetero! */
  717.     if( IsSolvent(CurGroup->refno) )
  718.     heta = True;
  719.  
  720.     if( heta )
  721.     {   HetaGroupCount++;
  722.     if( HMinMaxFlag )
  723.     {   if( serno > MaxHetaRes ) 
  724.         {   MaxHetaRes = serno;
  725.         } else if( serno < MinHetaRes )
  726.         MinHetaRes = serno;
  727.     } else MinHetaRes = MaxHetaRes = serno;
  728.     } else 
  729.     {   MainGroupCount++;
  730.     if( MMinMaxFlag )
  731.     {   if( serno > MaxMainRes )
  732.         {   MaxMainRes = serno;
  733.         } else if( serno < MinMainRes )
  734.         MinMainRes = serno;
  735.     } else MinMainRes = MaxMainRes = serno; 
  736.     }
  737. }
  738.  
  739.  
  740. static void ProcessPDBAtom( heta )
  741.     int heta;
  742. {
  743.     auto char name[4];
  744.     register Bond __far *bptr;
  745.     register Atom __far *ptr;
  746.     register Long dx,dy,dz;
  747.     register int serno;
  748.     register int refno;
  749.     register int i;
  750.  
  751.     /* Ignore Pseudo Atoms!! */
  752.     if( (Record[13]==' ') && (Record[14]=='Q') )
  753.     return; 
  754.  
  755.     dx = ReadPDBCoord(31);
  756.     dy = ReadPDBCoord(39);
  757.     dz = ReadPDBCoord(47);
  758.  
  759.     /* Ignore XPLOR Pseudo Atoms!! */
  760.     if( (dx==9999000L) && (dy==9999000L) && (dz==9999000L) )
  761.     return;
  762.  
  763.     serno = (int)ReadValue(23,4);
  764.     if( !CurGroup || (CurGroup->serno!=serno) 
  765.     || (CurChain->ident!=Record[22]) 
  766.     || (PDBInsert!=Record[27]) )
  767.     ProcessPDBGroup( heta, serno );
  768.  
  769.     /* Solvents should be hetero! */
  770.     if( IsSolvent(CurGroup->refno) )
  771.     heta = True;
  772.  
  773.  
  774.     ptr = CreateAtom();
  775.     if( isdigit(Record[14]) )
  776.     {   name[1] = ToUpper(Record[13]);
  777.     name[2] = name[3] = ' ';
  778.     name[0] = ' ';
  779.  
  780.     } else for( i=0; i<4; i++ )
  781.     name[i] = ToUpper(Record[i+13]);
  782.  
  783.     for( refno=0; refno<ElemNo; refno++ )
  784.     if( !strncmp(ElemDesc[refno],name,4) )
  785.         break;
  786.  
  787.     if( refno==ElemNo )
  788.     {   if( ElemNo++ == MAXELEM )
  789.         FatalDataError("Too many new atom types");
  790.  
  791.     for( i=0; i<4; i++ )
  792.         ElemDesc[refno][i] = name[i];
  793.     }
  794.  
  795.     ptr->refno = refno;
  796.     ptr->altl = Record[17];
  797.     ptr->serno = (int)ReadValue(7,5);
  798.     ptr->temp = (int)ReadValue(61,6);
  799.  
  800.     ptr->xorg =  dx/4;
  801.     ptr->yorg =  dy/4;
  802.     ptr->zorg = -dz/4;
  803.  
  804.     if( heta ) 
  805.     ptr->flag |= HeteroFlag;
  806.  
  807.     ProcessAtom( ptr );
  808.     if( IsAlphaCarbon(refno) && IsProtein(CurGroup->refno) )
  809.     {   if( ConnectAtom )
  810.     {   dx = ConnectAtom->xorg - ptr->xorg;
  811.         dy = ConnectAtom->yorg - ptr->yorg;
  812.         dz = ConnectAtom->zorg - ptr->zorg;
  813.  
  814.         /* Break backbone if CA-CA > 7.00A */
  815.         if( dx*dx+dy*dy+dz*dz < (Long)1750*1750 )
  816.         {   bptr = ProcessBond(ptr,ConnectAtom,NormBondFlag);
  817.         bptr->bnext = CurChain->blist;
  818.         CurChain->blist = bptr;
  819.         } else ptr->flag |= BreakFlag;
  820.     }
  821.     ConnectAtom = ptr;
  822.     } else if( IsSugarPhosphate(refno) && IsNucleo(CurGroup->refno) )
  823.     {   if( ConnectAtom )
  824.     {   bptr = ProcessBond(ConnectAtom,ptr,NormBondFlag);
  825.         bptr->bnext = CurChain->blist;
  826.         CurChain->blist = bptr;
  827.     }
  828.     ConnectAtom = ptr;
  829.     }
  830. }
  831.  
  832.  
  833. static FeatEntry __far *AllocFeature()
  834. {
  835.     register Feature __far *ptr;
  836.  
  837.     if( !FeatList || (FeatList->count==FeatSize) )
  838.     {   ptr = (Feature __far*)_fmalloc(sizeof(Feature));
  839.     if( !ptr ) FatalDataError("Memory allocation failed");
  840.     /* Features are always deallocated! */
  841.     
  842.     ptr->fnext = FeatList;
  843.     ptr->count = 0;
  844.     FeatList = ptr;
  845.     } else ptr = FeatList;
  846.  
  847.     return( &(ptr->data[ptr->count++]) );
  848. }
  849.  
  850.  
  851. #ifdef FUNCPROTO
  852. static void UpdateFeature( FeatEntry __far*, int );
  853. #endif
  854.  
  855. static void UpdateFeature( ptr, mask )
  856.     FeatEntry __far *ptr;  int mask;
  857. {
  858.     register Chain __far *chain;
  859.     register Group __far *group;
  860.  
  861.     for( chain=Database->clist; chain; chain=chain->cnext )
  862.     if( chain->ident == ptr->chain )
  863.     {   group=chain->glist;
  864.         while( group && (group->serno<ptr->init) )
  865.         group = group->gnext;
  866.  
  867.         while( group && (group->serno<=ptr->term) )
  868.         {   group->struc |= mask;
  869.         group = group->gnext;
  870.         }
  871.         return;
  872.     }
  873. }
  874.  
  875.  
  876. static void ProcessPDBFeatures()
  877. {
  878.     register Feature __far *next;
  879.     register Feature __far *ptr;
  880.     register int i;
  881.  
  882.     InfoTurnCount = 0;
  883.     InfoHelixCount = 0;
  884.     InfoLadderCount = 0;
  885.     InfoStrucSource = SourcePDB;
  886.  
  887.     for( ptr=FeatList; ptr; ptr=next )
  888.     {    if( Database )
  889.          for( i=0; i<ptr->count; i++ )
  890.          if( ptr->data[i].type==FeatHelix )
  891.          {   UpdateFeature( &ptr->data[i], HelixFlag );
  892.              InfoHelixCount++;
  893.          } else if( ptr->data[i].type==FeatSheet )
  894.          {   UpdateFeature( &ptr->data[i], SheetFlag );
  895.              InfoLadderCount++;
  896.          } else /* FeatTurn */
  897.          {   UpdateFeature( &ptr->data[i], TurnFlag );
  898.              InfoTurnCount++;
  899.          }
  900.  
  901.      /* Deallocate Memory */
  902.      next = ptr->fnext;
  903.      _ffree( ptr );
  904.     }
  905. }
  906.  
  907.  
  908. int LoadPDBMolecule( fp )
  909.     FILE *fp;
  910. {
  911.     register FeatEntry __far *ptr;
  912.     register int srcatm, dstatm;
  913.     register char *src, *dst;
  914.     register int model;
  915.  
  916.     model = True;
  917.     FeatList = (void __far*)0;
  918.     DataFile = fp;
  919.  
  920.  
  921.     do {   
  922.     FetchRecord();
  923.  
  924.     if( !strncmp("ATOM",Record+1,4) )
  925.     {   if( model ) ProcessPDBAtom(False);
  926.     } else if( !strncmp("HETA",Record+1,4) )
  927.     {   if( model ) ProcessPDBAtom(True);
  928.     } else if( !strncmp("SHEE",Record+1,4) )
  929.     {   if( !model ) continue;
  930.     
  931.         /* Remaining SHEET record fields   */
  932.         /* 39-40 .... Strand Parallelism   */
  933.         /* 33 ....... Same Chain as 22?    */
  934.         ptr = AllocFeature();
  935.         ptr->type = FeatSheet;
  936.         ptr->chain = Record[22];
  937.         ptr->init = (int)ReadValue(23,4);
  938.         ptr->term = (int)ReadValue(34,4);
  939.        
  940.     } else if( !strncmp("HELI",Record+1,4) )
  941.     {   if( !model ) continue;
  942.     
  943.         /* Remaining HELIX record fields   */
  944.         /* 39-40 .... Helix Classification */
  945.         /* 32 ....... Same Chain as 20?    */
  946.         ptr = AllocFeature();
  947.         ptr->type = FeatHelix;
  948.         ptr->chain = Record[20];
  949.         ptr->init = (int)ReadValue(22,4);
  950.         ptr->term = (int)ReadValue(34,4);
  951.  
  952.     } else if( !strncmp("TURN",Record+1,4) )
  953.     {   if( !model ) continue;
  954.     
  955.         ptr = AllocFeature();
  956.         ptr->type = FeatTurn;
  957.         ptr->chain = Record[20];
  958.         ptr->init = (int)ReadValue(21,4);
  959.         ptr->term = (int)ReadValue(32,4);
  960.  
  961.     } else if( !strncmp("CONE",Record+1,4) )
  962.     {   if( !model ) continue;
  963.  
  964.         if( (srcatm = (int)ReadValue(7,5)) )
  965.         {   if( (dstatm = (int)ReadValue(12,5)) )
  966.             if( dstatm > srcatm )
  967.             PDBConnectAtoms( srcatm, dstatm );
  968.  
  969.         if( !Record[17] ) continue;
  970.         if( (dstatm = (int)ReadValue(17,5)) )
  971.             if( dstatm > srcatm )
  972.             PDBConnectAtoms( srcatm, dstatm );
  973.  
  974.         if( !Record[22] ) continue;
  975.         if( (dstatm = (int)ReadValue(22,5)) )
  976.             if( dstatm > srcatm )
  977.             PDBConnectAtoms( srcatm, dstatm );
  978.  
  979.         if( !Record[27] ) continue;
  980.         if( (dstatm = (int)ReadValue(27,5)) )
  981.             if( dstatm > srcatm )
  982.             PDBConnectAtoms( srcatm, dstatm );
  983.         }
  984.  
  985.     } else if( !strncmp("COLO",Record+1,4) )
  986.     {   ProcessPDBColourMask();
  987.     } else if( !strncmp("TER",Record+1,3) )
  988.     {   if( !Record[4] || (Record[4]==' ') )
  989.         {   ConnectAtom = (void __far*)0;
  990.         CurGroup = (void __far*)0;
  991.         CurChain = (void __far*)0;
  992.         }
  993.     } else if( !strncmp("HEAD",Record+1,4) )
  994.     {   ExtractString(40,Record+11,InfoClassification);
  995.         ExtractString( 4,Record+63,InfoIdentCode);
  996.         
  997.     } else if( !strncmp("COMP",Record+1,4) )
  998.     {   if( Record[10]==' ' )  /* First COMPND record */
  999.         ExtractString(60,Record+11,InfoMoleculeName);
  1000.  
  1001.     } else if( !strncmp("CRYS",Record+1,4) )
  1002.     {   dst = InfoSpaceGroup;
  1003.         for( src=Record+56; *src && src<Record+67; src++ )
  1004.         if( *src!=' ' ) *dst++ = *src;
  1005.         *dst = 0;
  1006.  
  1007.         InfoCellA = ReadValue( 7,9)/1000.0;
  1008.         InfoCellB = ReadValue(16,9)/1000.0;
  1009.         InfoCellC = ReadValue(25,9)/1000.0;
  1010.  
  1011.         InfoCellAlpha = Deg2Rad*(ReadValue(34,7)/100.0);
  1012.         InfoCellBeta =  Deg2Rad*(ReadValue(41,7)/100.0);
  1013.         InfoCellGamma = Deg2Rad*(ReadValue(48,7)/100.0);
  1014.  
  1015.     } else if( !strncmp("ENDM",Record+1,4) )
  1016.     {   model = False;
  1017.     } else if( !strncmp("END",Record+1,3) )
  1018.         if( !Record[4] || (Record[4]==' ') )
  1019.         {   /* Treat END same as TER! */
  1020.         ConnectAtom = (void __far*)0;
  1021.         CurGroup = (void __far*)0;
  1022.         CurChain = (void __far*)0;
  1023.         }
  1024.     } while( !feof(DataFile) );
  1025.  
  1026.     if( Database )
  1027.     strcpy(InfoFileName,DataFileName);
  1028.     if( FeatList ) ProcessPDBFeatures();
  1029.     return( True );
  1030. }
  1031.  
  1032.  
  1033. static int SimpleAtomType( type )
  1034.     char *type;
  1035. {
  1036.     register int i, refno;
  1037.     auto char name[4];
  1038.  
  1039.     name[2] = name[3] = ' ';
  1040.     if( type[1] && (type[1]!=' ') )
  1041.     {   name[0] = ToUpper(type[0]);
  1042.     name[1] = ToUpper(type[1]);
  1043.     } else
  1044.     {   name[1] = ToUpper(type[0]);
  1045.     name[0] = ' ';
  1046.     }
  1047.  
  1048.     for( refno=0; refno<ElemNo; refno++ )
  1049.     if( !strncmp(ElemDesc[refno],name,4) )
  1050.         return( refno );
  1051.  
  1052.     if( ElemNo++ == MAXELEM )
  1053.     FatalDataError("Too many new atom types");
  1054.  
  1055.     for( i=0; i<4; i++ )
  1056.     ElemDesc[refno][i] = name[i];
  1057.     return( refno );
  1058. }
  1059.  
  1060.  
  1061. static void CreateMolGroup()
  1062. {
  1063.     strcpy(InfoFileName,DataFileName);
  1064.  
  1065.     CreateChain( ' ' );
  1066.     CreateGroup( 1 );
  1067.  
  1068.     CurGroup->refno = FindResNo( "MOL" );
  1069.     CurGroup->serno = 1;
  1070.     
  1071.     MinMainRes = MaxMainRes = 1;
  1072.     MinHetaRes = MaxHetaRes = 0;
  1073.     MainGroupCount = 1;
  1074. }
  1075.  
  1076.  
  1077. int LoadXYZMolecule( fp )
  1078.     FILE *fp;
  1079. {
  1080.     auto char type[12];
  1081.     auto double xpos, ypos, zpos;
  1082.     auto double charge, u, v, w;
  1083.     auto int atoms;
  1084.  
  1085.     register Atom __far *ptr;
  1086.     register char *src,*dst;
  1087.     register int i,count;
  1088.  
  1089.  
  1090.     DataFile = fp;
  1091.  
  1092.     /* Number of Atoms */
  1093.     FetchRecord();
  1094.     sscanf(Record+1,"%d",&atoms);
  1095.  
  1096.     /* Molecule (step) Description */
  1097.     FetchRecord();
  1098.     src = Record+1;
  1099.     while( *src == ' ' )
  1100.     src++;
  1101.  
  1102.     dst = InfoMoleculeName;
  1103.     for( i=0; i<78; i++ )
  1104.     if( *src ) *dst++ = *src++;
  1105.     *dst = '\0';
  1106.  
  1107.     if( atoms )
  1108.     {   CreateMolGroup();
  1109.     for( i=0; i<atoms; i++ )
  1110.     {   FetchRecord();
  1111.         ptr = CreateAtom();
  1112.         ptr->serno = i;
  1113.  
  1114.         xpos = ypos = zpos = 0.0;
  1115.         count = sscanf(Record+1,"%s %lg %lg %lg %lg %lg %lg %lg",
  1116.                type, &xpos, &ypos, &zpos, &charge, &u, &v, &w );
  1117.  
  1118.         ptr->refno = SimpleAtomType(type);
  1119.         ptr->xorg =  (Long)(250.0*xpos);
  1120.         ptr->yorg =  (Long)(250.0*ypos);
  1121.         ptr->zorg = -(Long)(250.0*zpos);
  1122.  
  1123.         if( (count==5) || (count==8) )
  1124.         {   ptr->temp = (short)(100.0*charge);
  1125.         } else ptr->temp = 0;
  1126.         ProcessAtom( ptr );
  1127.     }
  1128.     }
  1129.     return( True );
  1130. }
  1131.  
  1132.  
  1133. static int FindSybylResNo( ptr )
  1134.     char *ptr;
  1135. {
  1136.     register char *src,*dst;
  1137.     register int i, j;
  1138.     auto char name[4];
  1139.  
  1140.     src = ptr;
  1141.     dst = name;
  1142.     if( ptr[1] && (ptr[1]!='.') )
  1143.     {   *dst++ = ToUpper(*src);  src++;
  1144.     *dst++ = ToUpper(*src);  src++;
  1145.     } else
  1146.     {   *dst++ = ToUpper(*src);  src++;
  1147.     *dst++ = ' ';
  1148.     }
  1149.  
  1150.     if( *src )
  1151.     {   src++;
  1152.  
  1153.     if( *src == 'a' )
  1154.     {   *dst++ = ' ';
  1155.         *dst = ' ';
  1156.     } else if( *src == 'p' )
  1157.     {   *dst++ = '3';
  1158.         *dst = ' ';
  1159.     } else
  1160.     {   *dst++ = *src++;
  1161.         if( *src && (*src!='+') )
  1162.         {   *dst = *src;
  1163.         } else *dst = ' ';
  1164.     }
  1165.     } else
  1166.     {   *dst++ = ' ';
  1167.     *dst = ' ';
  1168.     }
  1169.  
  1170.     for( j=0; j<ElemNo; j++ )
  1171.     if( !strncmp(ElemDesc[j],name,4) )
  1172.         return( j );
  1173.  
  1174.     if( ElemNo++ == MAXELEM )
  1175.     FatalDataError("Too many new atom types");
  1176.     for( i=0; i<4; i++ ) ElemDesc[j][i] = name[i];
  1177.     return( j );
  1178. }
  1179.  
  1180.  
  1181. int LoadMol2Molecule( fp )
  1182.     FILE *fp;
  1183. {
  1184.     auto int srcatm, dstatm;
  1185.     auto int features, sets, serno;
  1186.     auto int atoms, bonds, structs;
  1187.     auto double xpos, ypos, zpos;
  1188.  
  1189.     auto char name[8];
  1190.     auto char type[4];
  1191.  
  1192.     register Atom __far *ptr;
  1193.     register char *src, *dst;
  1194.     register int i;
  1195.  
  1196.  
  1197.     DataFile = fp;
  1198.  
  1199.     while( !feof(DataFile) )
  1200.     {   FetchRecord();
  1201.     if( !Record[1] || Record[1]=='#' )
  1202.         continue;
  1203.  
  1204.     if( !strncmp("@<TRIPOS>MOLECULE",Record+1,17) )
  1205.     {   FetchRecord();  /* Molecule Name */
  1206.         src = Record+1;
  1207.         while( *src==' ' )
  1208.         src++;
  1209.  
  1210.         dst = InfoMoleculeName;
  1211.         while( (*dst++ = *src++) );
  1212.  
  1213.         FetchRecord();
  1214.         atoms = bonds = structs = features = sets = 0;
  1215.         sscanf(Record+1,"%d %d %d %d %d", &atoms, &bonds, &structs,
  1216.                           &features, &sets );
  1217.  
  1218.         FetchRecord();  /* Molecule Type  */
  1219.         FetchRecord();  /* Charge Type    */
  1220.  
  1221.     } else if( !strncmp("@<TRIPOS>ATOM",Record+1,13) )
  1222.     {   if( !atoms ) continue;
  1223.  
  1224.         CreateMolGroup();
  1225.         for( i=0; i<atoms; i++ )
  1226.         {    FetchRecord();
  1227.          ptr = CreateAtom();
  1228.  
  1229.          sscanf(Record+1,"%d %s %lg %lg %lg %s", &serno, name,
  1230.                   &xpos, &ypos, &zpos, type );
  1231.  
  1232.          ptr->refno = FindSybylResNo( type );
  1233.          ptr->serno = serno;
  1234.          /* ptr->serno = i; */
  1235.  
  1236.          ptr->xorg =  (Long)(250.0*xpos);
  1237.          ptr->yorg =  (Long)(250.0*ypos);
  1238.          ptr->zorg = -(Long)(250.0*zpos);
  1239.          ProcessAtom( ptr );
  1240.         }
  1241.  
  1242.     } else if( !strncmp("@<TRIPOS>BOND",Record+1,13) )
  1243.         for( i=0; i<bonds; i++ )
  1244.         {   FetchRecord();
  1245.         sscanf(Record+1,"%d %d %d %.2s",&serno,&srcatm,&dstatm,type);
  1246.         if( !strncmp(type,"ar",2) )
  1247.         {   ConnectAtoms(srcatm,dstatm,AromBondFlag);
  1248.         } else if( *type == '2' )
  1249.         {   ConnectAtoms(srcatm,dstatm,DoubBondFlag);
  1250.         } else /* *type == '1' */
  1251.             ConnectAtoms(srcatm,dstatm,NormBondFlag);
  1252.         }
  1253.     }
  1254.     return( True );
  1255. }
  1256.  
  1257.  
  1258.  
  1259. static int FindAlchemyResNo()
  1260. {
  1261.     register char *ptr;
  1262.     register int i, j;
  1263.     char name[4];
  1264.  
  1265.     ptr = Record+7;
  1266.     if( !isalpha(ptr[1]) )
  1267.     {   name[0] = ' ';
  1268.     for( i=0; i<3; i++ )
  1269.         name[i+1] = ToUpper(ptr[i]);
  1270.     ptr = name;
  1271.     } else
  1272.     {   for( i=0; i<4; i++ )
  1273.         ptr[i] = ToUpper(ptr[i]);
  1274.  
  1275.     for( i=0; i<MAXALCATOM; i++ )
  1276.         if( !strncmp(AlcAtomTable[i].src,ptr,4) )
  1277.         {   ptr = AlcAtomTable[i].dst;
  1278.         break;
  1279.         }
  1280.     }
  1281.  
  1282.     for( j=0; j<ElemNo; j++ )
  1283.     if( !strncmp(ElemDesc[j],ptr,4) )
  1284.         return( j );
  1285.  
  1286.     if( ElemNo++ == MAXELEM )
  1287.     FatalDataError("Too many new atom types");
  1288.     for( i=0; i<4; i++ ) ElemDesc[j][i] = ptr[i];
  1289.     return( j );
  1290. }
  1291.  
  1292.  
  1293. int LoadAlchemyMolecule( fp )
  1294.     FILE *fp;
  1295. {
  1296.     auto char type[12];
  1297.     auto int serno,srcatm,dstatm;
  1298.     register Atom __far *ptr;
  1299.     register int atoms, bonds;
  1300.     register int i;
  1301.  
  1302.     DataFile = fp;
  1303.     FetchRecord();
  1304.     atoms = (int)ReadValue(1,5);
  1305.     bonds = (int)ReadValue(14,5);
  1306.     ExtractString(38,Record+42,InfoMoleculeName);
  1307.  
  1308.     if( !atoms )
  1309.     return( False );
  1310.  
  1311.     CreateMolGroup();
  1312.     for( i=0; i<atoms; i++ )
  1313.     {   FetchRecord();
  1314.     ptr = CreateAtom();
  1315.  
  1316.     ptr->refno = FindAlchemyResNo();
  1317.     ptr->temp = (int)ReadValue(41,8);
  1318.     ptr->serno = (int)ReadValue(1,5);
  1319.     /* ptr->serno = i+1; */
  1320.  
  1321.     ptr->xorg =  ReadValue(13,7)/4;
  1322.     ptr->yorg =  ReadValue(22,7)/4;
  1323.     ptr->zorg = -ReadValue(31,7)/4;
  1324.     ProcessAtom( ptr );
  1325.     }
  1326.  
  1327.     for( i=0; i<bonds; i++ )
  1328.     {   FetchRecord();
  1329.     sscanf(Record+1,"%d %d %d %.10s",&serno,&srcatm,&dstatm,type);
  1330.  
  1331.     if( *type =='A' )             /* AROMATIC */
  1332.     {   ConnectAtoms(srcatm,dstatm,AromBondFlag);
  1333.     } else if( *type == 'D' )     /* DOUBLE */
  1334.     {   ConnectAtoms(srcatm,dstatm,DoubBondFlag);
  1335.     } else if( *type == 'T' )     /* TRIPLE */
  1336.     {   ConnectAtoms(srcatm,dstatm,TripBondFlag);
  1337.     } else /* (*type == 'S') */   /* SINGLE */
  1338.         ConnectAtoms(srcatm,dstatm,NormBondFlag);
  1339.     }
  1340.     return( True );
  1341. }
  1342.  
  1343.  
  1344. int LoadCharmmMolecule( fp )
  1345.     FILE *fp;
  1346. {
  1347.     auto char buffer[9];
  1348.     register Atom __far *ptr;
  1349.     register int refno,serno;
  1350.     register int resno,atoms;
  1351.     register int init,chain;
  1352.     register int i;
  1353.  
  1354.     DataFile = fp;
  1355.  
  1356.     do { 
  1357.     FetchRecord();
  1358.     } while( Record[1]=='*' );
  1359.     atoms = (int)ReadValue(1,5);
  1360.     if( !atoms ) return False;
  1361.  
  1362.     MinHetaRes = MaxHetaRes = 0;
  1363.     strcpy(InfoFileName,DataFileName);
  1364.     MainGroupCount = 0;
  1365.  
  1366.     chain = 0;
  1367.     init = False;
  1368.     CurChain = NULL;
  1369.     for( serno=0; serno<atoms; serno++ )
  1370.     {   FetchRecord();
  1371.  
  1372.     if( !CurChain || strncmp(Record+52,buffer,9) )
  1373.     {   for( i=0; i<9; i++ )
  1374.         buffer[i] = Record[52+i];
  1375.         CreateChain(chain+49);
  1376.         chain++;
  1377.     }
  1378.  
  1379.     resno = (int)ReadValue(6,5);
  1380.     if( !CurGroup || (CurGroup->serno!=resno) )
  1381.     {   CreateGroup( GroupPool );
  1382.         CurGroup->refno = FindResNo(Record+12);
  1383.         CurGroup->serno = resno;
  1384.  
  1385.         if( !init )
  1386.         {   MinMainRes = resno;
  1387.         MaxMainRes = resno;
  1388.         init = True;
  1389.         } else if( resno > MaxMainRes )
  1390.         {   MaxMainRes = resno;
  1391.         } else if( resno < MinMainRes )
  1392.         MinMainRes = resno;
  1393.         MainGroupCount++;
  1394.     }
  1395.  
  1396.     ptr = CreateAtom();
  1397.     for( refno=0; refno<ElemNo; refno++ )
  1398.         if( !strncmp(ElemDesc[refno],Record+16,4) )
  1399.         break;
  1400.  
  1401.     if( refno == ElemNo )
  1402.     {   if( ElemNo++ == MAXELEM )
  1403.         FatalDataError("Too many new atom types");
  1404.         for( i=0; i<4; i++ )
  1405.         ElemDesc[refno][i] = Record[16+i];
  1406.     }
  1407.  
  1408.     ptr->refno = refno;
  1409.     ptr->temp = (int)ReadValue(61,9);
  1410.     ptr->serno = (int)ReadValue(1,5);
  1411.     /* ptr->serno = serno+1; */
  1412.  
  1413.     ptr->xorg =  ReadValue(21,8)/4;
  1414.     ptr->yorg =  ReadValue(31,8)/4;
  1415.     ptr->zorg = -ReadValue(41,8)/4;
  1416.     ProcessAtom( ptr );
  1417.     }
  1418.     return( True );
  1419. }
  1420.  
  1421.  
  1422. int LoadMDLMolecule( fp )
  1423.     FILE *fp;
  1424. {
  1425.     register Bond __far *bptr;
  1426.     register Atom __far *src;
  1427.     register Atom __far *dst;
  1428.     register Atom __far *ptr;
  1429.  
  1430.     register int atoms, bonds;
  1431.     register int srcatm,dstatm;
  1432.     register int i,type,done;
  1433.     register Long dx, dy, dz;
  1434.     register Card dist2;
  1435.     register Real scale;
  1436.  
  1437.     DataFile = fp;
  1438.  
  1439.     FetchRecord(); /* Molecule Name */
  1440.     ExtractString(78,Record+1,InfoMoleculeName);
  1441.  
  1442.     FetchRecord(); /* Program Details */
  1443.     FetchRecord(); /* Comments */
  1444.  
  1445.     FetchRecord();
  1446.     atoms = (int)ReadValue(1,3);
  1447.     bonds = (int)ReadValue(4,3);
  1448.  
  1449.     if( !atoms )
  1450.     return( False );
  1451.  
  1452.     CreateMolGroup();
  1453.     for( i=1; i<=atoms; i++ )
  1454.     {   FetchRecord();
  1455.     ptr = CreateAtom();
  1456.     ptr->refno = SimpleAtomType(Record+32);
  1457.  
  1458.     switch( (int)ReadValue(37,3) )
  1459.     {   case(1):  ptr->temp =  300;  break;
  1460.         case(2):  ptr->temp =  200;  break;
  1461.         case(3):  ptr->temp =  100;  break;
  1462.         case(5):  ptr->temp = -100;  break;
  1463.         case(6):  ptr->temp = -200;  break;
  1464.         case(7):  ptr->temp = -300;  break;
  1465.         default:  ptr->temp = 0;
  1466.     }
  1467.     ptr->serno = i;
  1468.  
  1469.     ptr->xorg =  ReadValue(1,10)/40;
  1470.     ptr->yorg =  ReadValue(11,10)/40;
  1471.     ptr->zorg = -ReadValue(21,10)/40;
  1472.     ProcessAtom( ptr );
  1473.     }
  1474.  
  1475.     for( i=0; i<bonds; i++ )
  1476.     {   FetchRecord();
  1477.     srcatm = (int)ReadValue(1,3);
  1478.     dstatm = (int)ReadValue(4,3);
  1479.     type = (int)ReadValue(7,3);
  1480.  
  1481.     if( type==2 )                 /* DOUBLE */
  1482.     {   ConnectAtoms(srcatm,dstatm,DoubBondFlag);
  1483.     } else if( type==3 )          /* TRIPLE */
  1484.     {   ConnectAtoms(srcatm,dstatm,TripBondFlag);
  1485.     } else if( type==4 )          /* AROMATIC */
  1486.     {   ConnectAtoms(srcatm,dstatm,AromBondFlag);
  1487.     } else                        /* SINGLE */
  1488.         ConnectAtoms(srcatm,dstatm,NormBondFlag);
  1489.     }
  1490.  
  1491.     done = False;
  1492.     for( bptr=CurMolecule->blist; bptr; bptr=bptr->bnext )
  1493.     if( bptr->flag & NormBondFlag )
  1494.     {   src = bptr->srcatom;
  1495.         dst = bptr->dstatom;
  1496.         if( (src->refno==2) && (dst->refno==2) )
  1497.         {   dx = dst->xorg - src->xorg;
  1498.         dy = dst->yorg - src->yorg;
  1499.         dz = dst->zorg - src->zorg;
  1500.         if( dx || dy || dz )
  1501.         {   dist2 = dx*dx + dy*dy + dz*dz;
  1502.             scale = 385.0/sqrt(dist2);
  1503.             break;
  1504.         }
  1505.         }
  1506.     }
  1507.  
  1508.     if( bptr )
  1509.     {   for( ptr=CurGroup->alist; ptr; ptr=ptr->anext )
  1510.     {   ptr->xorg = (Long)(ptr->xorg*scale);
  1511.         ptr->yorg = (Long)(ptr->yorg*scale);
  1512.         ptr->zorg = (Long)(ptr->zorg*scale);
  1513.     }
  1514.     MinX = (Long)(MinX*scale);  MaxX = (Long)(MaxX*scale);
  1515.     MinY = (Long)(MinY*scale);  MaxY = (Long)(MaxY*scale);
  1516.     MinZ = (Long)(MinZ*scale);  MaxZ = (Long)(MaxZ*scale);
  1517.     }
  1518.     return( True );
  1519. }
  1520.  
  1521.  
  1522.  
  1523. int SavePDBMolecule( filename )
  1524.     char *filename;
  1525. {
  1526.     register double x, y, z;
  1527.     register Group __far *prev;
  1528.     register Chain __far *chain;
  1529.     register Group __far *group;
  1530.     register Atom __far *aptr;
  1531.     register char *ptr;
  1532.     register int count;
  1533.     register char ch;
  1534.     register int i;
  1535.  
  1536.     if( !Database )
  1537.     return( False );
  1538.  
  1539.     DataFile = fopen( filename, "w" );
  1540.     if( !DataFile )
  1541.     {   if( CommandActive )
  1542.         WriteChar('\n');
  1543.     WriteString("Error: Unable to create file!\n\n");
  1544.     CommandActive=False;
  1545.     return( False );
  1546.     }
  1547.  
  1548.     if( *InfoClassification || *InfoIdentCode )
  1549.     {   fputs("HEADER    ",DataFile);
  1550.  
  1551.     ptr = InfoClassification;
  1552.     for( i=11; i<=50; i++ )
  1553.         putc( (*ptr ? *ptr++ : ' '), DataFile );
  1554.     fprintf(DataFile,"13-JUL-93   %.4s\n",InfoIdentCode);
  1555.     }
  1556.  
  1557.     if( *InfoMoleculeName )
  1558.     fprintf(DataFile,"COMPND    %.60s\n",InfoMoleculeName);
  1559.  
  1560.     prev = (void __far*)0;
  1561.  
  1562.     count = 1;
  1563.     ForEachAtom
  1564.     if( aptr->flag&SelectFlag )
  1565.     {   if( prev && (chain->ident!=ch) )
  1566.         fprintf( DataFile, "TER   %5d      %.3s %c%4d \n", 
  1567.              count++, Residue[prev->refno], ch, prev->serno);
  1568.  
  1569.         if( aptr->flag&HeteroFlag )
  1570.         {      fputs("HETATM",DataFile);
  1571.         } else fputs("ATOM  ",DataFile);
  1572.         fprintf( DataFile, "%5d %.4s %.3s %c%4d    ",
  1573.              count++, ElemDesc[aptr->refno], Residue[group->refno],
  1574.              chain->ident, group->serno );
  1575.  
  1576.         x = (double)aptr->xorg/250.0;
  1577.         y = (double)aptr->yorg/250.0;
  1578.         z = (double)aptr->zorg/250.0;
  1579.  
  1580. #ifdef INVERT
  1581.         fprintf(DataFile,"%8.3f%8.3f%8.3f",x,-y,-z);
  1582. #else
  1583.         fprintf(DataFile,"%8.3f%8.3f%8.3f",x,y,-z);
  1584. #endif
  1585.         fprintf(DataFile,"  1.00%6.2f\n",aptr->temp/100.0);
  1586.         
  1587.         ch = chain->ident;
  1588.         prev = group;
  1589.     }
  1590.  
  1591.     if( prev )
  1592.     fprintf( DataFile, "TER   %5d      %.3s %c%4d \n", 
  1593.          count, Residue[prev->refno], ch, prev->serno);
  1594.  
  1595.     fputs("END   \n",DataFile);
  1596.     fclose( DataFile );
  1597. #ifdef APPLEMAC
  1598.     SetFileInfo(filename,'RSML','TEXT',131);
  1599. #endif
  1600.     return( True );
  1601. }
  1602.  
  1603. int SaveXYZMolecule( filename )
  1604.     char *filename;
  1605. {
  1606.     return( True );
  1607. }
  1608.  
  1609.  
  1610. int SaveCIFMolecule( filename )
  1611.     char *filename;
  1612. {
  1613.     if( !filename )
  1614.         return( False );
  1615.     return( True );
  1616. }
  1617.  
  1618.  
  1619. int SaveAlchemyMolecule( filename )
  1620.     char *filename;
  1621. {
  1622.     register Real x, y, z;
  1623.     register float xpos, ypos, zpos;
  1624.     register Chain __far *chain;
  1625.     register Group __far *group;
  1626.     register Atom __far *aptr;
  1627.     register Bond __far *bptr;
  1628.     register char *ptr;
  1629.     register int atomno;
  1630.     register int bondno;
  1631.     register int num;
  1632.  
  1633.     if( !Database )
  1634.     return( False );
  1635.  
  1636.     DataFile = fopen( filename, "w" );
  1637.     if( !DataFile )
  1638.     {   if( CommandActive )
  1639.         WriteChar('\n');
  1640.     WriteString("Error: Unable to create file!\n\n");
  1641.     CommandActive=False;
  1642.     return( False );
  1643.     }
  1644.  
  1645.     atomno = 0;
  1646.     ForEachAtom
  1647.     if( aptr->flag & SelectFlag )
  1648.     {   aptr->mbox = 0;
  1649.         atomno++;
  1650.     }
  1651.  
  1652.     bondno = 0;
  1653.     ForEachBond
  1654.     if( (bptr->srcatom->flag&bptr->dstatom->flag) & SelectFlag )
  1655.     {   if( bptr->flag&AromBondFlag )
  1656.         {   bptr->srcatom->mbox = -1;
  1657.         bptr->dstatom->mbox = -1;
  1658.         } else if( !(bptr->flag&HydrBondFlag) )
  1659.         {   num = (bptr->flag&DoubBondFlag)? 2 : 1;
  1660.         if( bptr->srcatom->mbox>0 ) 
  1661.             bptr->srcatom->mbox += num;
  1662.         if( bptr->dstatom->mbox>0 ) 
  1663.             bptr->dstatom->mbox += num;
  1664.         }
  1665.         bondno++;
  1666.     }
  1667.  
  1668.     fprintf(DataFile,"%5d ATOMS, %5d BONDS, ",atomno,bondno);
  1669.     fprintf(DataFile,"    0 CHARGES, %s\n", InfoMoleculeName );
  1670.  
  1671.     atomno = 1;
  1672.     ForEachAtom
  1673.     if( aptr->flag & SelectFlag )
  1674.     {   aptr->mbox = atomno;
  1675.         fprintf(DataFile,"%5d ",atomno++);
  1676.  
  1677.         switch( GetElemNumber(aptr) )
  1678.         {   case( 6 ):  if( aptr->mbox == -1 )
  1679.                 {   ptr = "CAR ";
  1680.                 } else if( aptr->mbox == 1 )
  1681.                 {   ptr = "C3  ";
  1682.                 } else ptr = "C2  ";
  1683.                 fputs( ptr, DataFile );
  1684.                 break;
  1685.  
  1686.         case( 7 ):  if( aptr->mbox == -1 )
  1687.                 {   ptr = "NAR ";
  1688.                 } else ptr = "N2  ";
  1689.                 fputs( ptr, DataFile );
  1690.                 break;
  1691.  
  1692.         case( 8 ):  if( aptr->mbox == 2 )
  1693.                 {   ptr = "O2  ";
  1694.                 } else ptr = "O3  ";
  1695.                 fputs( ptr, DataFile );
  1696.                 break;
  1697.  
  1698.         case( 1 ):  fputs( "H   ", DataFile );  break;
  1699.  
  1700.         default:    ptr = ElemDesc[aptr->refno];
  1701.                 if( *ptr==' ' )
  1702.                 {   fprintf(DataFile,"%.3s ",ptr+1);
  1703.                 } else fprintf(DataFile,"%.4s",ptr);
  1704.         }
  1705.  
  1706.         x = aptr->xorg/250.0;
  1707.         y = aptr->yorg/250.0;
  1708.         z = aptr->zorg/250.0;
  1709.  
  1710.         /* Apply Current Viewpoint Rotation Matrix */
  1711.         xpos = (float)(x*RotX[0] + y*RotX[1] + z*RotX[2]);
  1712.         ypos = (float)(x*RotY[0] + y*RotY[1] + z*RotY[2]);
  1713.         zpos = (float)(x*RotZ[0] + y*RotZ[1] + z*RotZ[2]);
  1714.  
  1715. #ifdef INVERT
  1716.         fprintf(DataFile,"  %8.4f %8.4f %8.4f",xpos,-ypos,-zpos);
  1717. #else
  1718.         fprintf(DataFile,"  %8.4f %8.4f %8.4f",xpos,ypos,-zpos);
  1719. #endif
  1720.         fprintf(DataFile,"    %7.4f\n",aptr->temp/1000.0);
  1721.     }
  1722.  
  1723.     bondno = 1;
  1724.     ForEachBond
  1725.     if( (bptr->srcatom->flag&bptr->dstatom->flag) & SelectFlag )
  1726.     {   fprintf(DataFile,"%5d %5d %5d  ", bondno++,
  1727.             bptr->srcatom->mbox, bptr->dstatom->mbox );
  1728.         if( bptr->flag & AromBondFlag )
  1729.         {   ptr = "AROMATIC\n";
  1730.         } else if( bptr->flag & TripBondFlag )
  1731.         {   ptr = "TRIPLE\n";
  1732.         } else if( bptr->flag & DoubBondFlag )
  1733.         {   ptr = "DOUBLE\n";
  1734.         } else ptr = "SINGLE\n";
  1735.         fputs( ptr, DataFile );
  1736.     }
  1737.  
  1738.     ForEachAtom
  1739.         if( aptr->flag & SelectFlag )
  1740.             aptr->mbox = 0;
  1741.     fclose( DataFile );
  1742. #ifdef APPLEMAC
  1743.     SetFileInfo(filename,'RSML','TEXT',131);
  1744. #endif
  1745.     return( True );
  1746. }
  1747.  
  1748.  
  1749. static void TestBonded( sptr, dptr, flag )
  1750.     Atom __far *sptr, __far *dptr; 
  1751.     int flag;
  1752. {
  1753.     register Bond __far *bptr;
  1754.     register Long dx, dy, dz;
  1755.     register Long max, dist;
  1756.  
  1757.     if( flag )
  1758.     {    /* Sum of covalent radii with 0.5A tolerance */
  1759.          dist = Element[GetElemNumber(sptr)].covalrad + 
  1760.                 Element[GetElemNumber(dptr)].covalrad + 125;
  1761.          max = dist*dist;  
  1762.     } else 
  1763.     {    /* Fast Bio-Macromolecule Bonding Calculation */
  1764.          if( (sptr->flag|dptr->flag) & HydrogenFlag )
  1765.      {      max = MaxHBondDist;
  1766.          } else max = MaxBondDist;
  1767.     }
  1768.  
  1769.     dx = sptr->xorg-dptr->xorg;   if( (dist=dx*dx)>max ) return;
  1770.     dy = sptr->yorg-dptr->yorg;   if( (dist+=dy*dy)>max ) return;
  1771.     dz = sptr->zorg-dptr->zorg;   if( (dist+=dz*dz)>max ) return;
  1772.  
  1773.     if( dist > MinBondDist )
  1774.     {   /* Reset Non-bonded flags! */
  1775.     sptr->flag &= ~NonBondFlag;
  1776.     dptr->flag &= ~NonBondFlag;
  1777.  
  1778.         if( (sptr->flag|dptr->flag) & HydrogenFlag )
  1779.         {      flag = HydrBondFlag;
  1780.         } else flag = NormBondFlag;
  1781.     bptr = ProcessBond(sptr,dptr,flag);
  1782.     bptr->bnext = CurMolecule->blist;
  1783.     CurMolecule->blist = bptr;
  1784.     InfoBondCount++;
  1785.     }
  1786. }
  1787.  
  1788.  
  1789. static void ReclaimHBonds( ptr )
  1790.     HBond __far *ptr;
  1791. {
  1792.     register HBond __far *temp;
  1793.  
  1794.     if( (temp = ptr) )
  1795.     {   while( temp->hnext )
  1796.         temp=temp->hnext;
  1797.     temp->hnext = FreeHBond;
  1798.     FreeHBond = ptr;
  1799.     }
  1800. }
  1801.  
  1802.  
  1803. static void ReclaimBonds( ptr )
  1804.     Bond __far *ptr;
  1805. {
  1806.     register Bond __far *temp;
  1807.  
  1808.     if( (temp = ptr) )
  1809.     {   while( temp->bnext )
  1810.         temp=temp->bnext;
  1811.     temp->bnext = FreeBond;
  1812.     FreeBond = ptr;
  1813.     }
  1814. }
  1815.  
  1816.  
  1817. void CreateMoleculeBonds( info, flag )
  1818.     int info, flag;
  1819. {
  1820.     register int i, x, y, z;
  1821.     register Long tx, ty, tz;
  1822.     register Long mx, my, mz; 
  1823.     register Long dx, dy, dz;
  1824.     register int lx, ly, lz, hx, hy, hz;
  1825.     register Atom __far *aptr, __far *dptr;
  1826.     register Chain __far *chain;
  1827.     register Group __far *group;
  1828.     char buffer[40];
  1829.  
  1830.  
  1831.     if( !Database ) 
  1832.     return;
  1833.  
  1834.     dx = (MaxX-MinX)+1;
  1835.     dy = (MaxY-MinY)+1;
  1836.     dz = (MaxZ-MinZ)+1;
  1837.  
  1838.     InfoBondCount = 0;
  1839.     ReclaimBonds( CurMolecule->blist );
  1840.     CurMolecule->blist = (void __far*)0;
  1841.     ResetVoxelData();
  1842.  
  1843.     for( chain=Database->clist; chain; chain=chain->cnext )
  1844.     {   ResetVoxelData();
  1845.     for( group=chain->glist; group; group=group->gnext )
  1846.         for( aptr=group->alist; aptr; aptr=aptr->anext )
  1847.         {   /* Initially non-bonded! */
  1848.         aptr->flag |= NonBondFlag;
  1849.  
  1850.         mx = aptr->xorg-MinX;
  1851.         my = aptr->yorg-MinY;
  1852.         mz = aptr->zorg-MinZ;
  1853.  
  1854.         tx = mx-AbsMaxBondDist;  
  1855.         ty = my-AbsMaxBondDist;  
  1856.         tz = mz-AbsMaxBondDist;  
  1857.  
  1858.         lx = (tx>0)? (int)((VOXORDER*tx)/dx) : 0;
  1859.         ly = (ty>0)? (int)((VOXORDER*ty)/dy) : 0;
  1860.         lz = (tz>0)? (int)((VOXORDER*tz)/dz) : 0;
  1861.  
  1862.         tx = mx+AbsMaxBondDist;  
  1863.         ty = my+AbsMaxBondDist;  
  1864.         tz = mz+AbsMaxBondDist;
  1865.  
  1866.         hx = (tx<dx)? (int)((VOXORDER*tx)/dx) : VOXORDER-1;
  1867.         hy = (ty<dy)? (int)((VOXORDER*ty)/dy) : VOXORDER-1;
  1868.         hz = (tz<dz)? (int)((VOXORDER*tz)/dz) : VOXORDER-1;
  1869.     
  1870.         for( x=lx; x<=hx; x++ )
  1871.         {   i = VOXORDER2*x + VOXORDER*ly;
  1872.             for( y=ly; y<=hy; y++ )
  1873.             {   for( z=lz; z<=hz; z++ )
  1874.                 if( (dptr = (Atom __far*)HashTable[i+z]) )
  1875.                 do { TestBonded(aptr,dptr,flag);
  1876.                 } while( (dptr = dptr->next) );
  1877.             i += VOXORDER;
  1878.             }
  1879.         }
  1880.         
  1881.         x = (int)((VOXORDER*mx)/dx);
  1882.         y = (int)((VOXORDER*my)/dy);
  1883.         z = (int)((VOXORDER*mz)/dz);
  1884.  
  1885.         i = VOXORDER2*x + VOXORDER*y + z;
  1886.         aptr->next = (Atom __far*)HashTable[i];
  1887.         HashTable[i] = (void __far*)aptr;
  1888.         }
  1889.     VoxelsClean = False;
  1890.     }
  1891.  
  1892.     if( info )
  1893.     {   if( CommandActive )
  1894.         WriteChar('\n');
  1895.     CommandActive=False;
  1896.     sprintf(buffer,"Number of Bonds ..... %ld\n\n",InfoBondCount);
  1897.     WriteString(buffer);
  1898.     }
  1899. }
  1900.  
  1901.  
  1902. Atom __far *FindGroupAtom( group, n )
  1903.     Group __far *group;  Byte n;
  1904. {
  1905.     register Atom __far *ptr;
  1906.  
  1907.     ptr = group->alist;
  1908.     while( ptr )
  1909.     {   if( ptr->refno == n )
  1910.         return( ptr );
  1911.     ptr = ptr->anext;
  1912.     }
  1913.     return( (Atom __far*)0 );
  1914. }
  1915.  
  1916.  
  1917. void TestDisulphideBridge( group1, group2, cys1 )
  1918.     Group __far *group1, __far *group2;  
  1919.     Atom __far *cys1;
  1920. {
  1921.     register HBond __far *ptr;
  1922.     register Atom __far *cys2;
  1923.     register int dx, dy, dz;
  1924.     register Long max,dist;
  1925.  
  1926.     if( !(cys2=FindGroupAtom(group2,20)) )
  1927.     return;
  1928.  
  1929.     max = (Long)750*750;
  1930.     dx = (int)(cys1->xorg-cys2->xorg);   if( (dist=(Long)dx*dx)>max ) return;
  1931.     dy = (int)(cys1->yorg-cys2->yorg);   if( (dist+=(Long)dy*dy)>max ) return;
  1932.     dz = (int)(cys1->zorg-cys2->zorg);   if( (dist+=(Long)dz*dz)>max ) return;
  1933.  
  1934.     if( !(ptr = FreeHBond) )
  1935.     {   MemSize += sizeof(HBond);
  1936.     ptr = (HBond __far*)_fmalloc(sizeof(HBond));
  1937.     if( !ptr ) FatalDataError("Memory allocation failed");
  1938.     RegisterAlloc( ptr );
  1939.     } else FreeHBond = ptr->hnext;
  1940.  
  1941.     ptr->dst = cys1;
  1942.     if( !(ptr->dstCA=FindGroupAtom(group1,1)) )
  1943.     ptr->dstCA = cys1;
  1944.  
  1945.     ptr->src = cys2;
  1946.     if( !(ptr->srcCA=FindGroupAtom(group2,1)) )
  1947.     ptr->srcCA = cys2;
  1948.  
  1949.     ptr->offset = 0;
  1950.     ptr->energy = 0;
  1951.     ptr->flag = 0;
  1952.     ptr->col = 0;
  1953.  
  1954.     ptr->hnext = CurMolecule->slist;
  1955.     CurMolecule->slist = ptr;
  1956.  
  1957.     group1->flag |= CystineFlag;
  1958.     group2->flag |= CystineFlag;
  1959.     InfoSSBondCount++;
  1960. }
  1961.  
  1962.  
  1963. void FindDisulphideBridges()
  1964. {
  1965.     register Chain __far *chn1;
  1966.     register Chain __far *chn2;
  1967.     register Group __far *group1;
  1968.     register Group __far *group2;
  1969.     register Atom __far *cys;
  1970.     char buffer[40];
  1971.  
  1972.     if( !Database ) return;
  1973.     ReclaimHBonds( CurMolecule->slist );
  1974.     InfoSSBondCount = 0;
  1975.  
  1976.     for(chn1=Database->clist;chn1;chn1=chn1->cnext)
  1977.     for(group1=chn1->glist;group1;group1=group1->gnext)
  1978.         if( IsCysteine(group1->refno) && (cys=FindGroupAtom(group1,20)) )
  1979.         {   for(group2=group1->gnext;group2;group2=group2->gnext)
  1980.             if( IsCysteine(group2->refno) )
  1981.             TestDisulphideBridge(group1,group2,cys);
  1982.  
  1983.         for(chn2=chn1->cnext;chn2;chn2=chn2->cnext)
  1984.             for(group2=chn2->glist;group2;group2=group2->gnext)
  1985.             if( IsCysteine(group2->refno) )
  1986.                 TestDisulphideBridge(group1,group2,cys);
  1987.         }
  1988.  
  1989.     if( CommandActive )
  1990.     WriteChar('\n');
  1991.     CommandActive=False;
  1992.     
  1993.     sprintf(buffer,"Number of Bridges ... %d\n\n",InfoSSBondCount);
  1994.     WriteString(buffer);
  1995. }
  1996.  
  1997.  
  1998. #ifdef FUNCPROTO
  1999. static void CreateHydrogenBond( Atom __far*, Atom __far*,
  2000.                 Atom __far*, Atom __far*, int, int );
  2001. static int IsHBonded( Atom __far*, Atom __far*, HBond __far* );
  2002. static void TestLadder( Chain __far* );
  2003. #endif
  2004.  
  2005.  
  2006. static void CreateHydrogenBond( srcCA, dstCA, src, dst, energy, offset )
  2007.     Atom __far *srcCA, __far *dstCA;
  2008.     Atom __far *src, __far *dst;
  2009.     int energy, offset;
  2010. {
  2011.     register HBond __far *ptr;
  2012.     register int i,flag;
  2013.  
  2014.     if( !(ptr = FreeHBond) )
  2015.     {   MemSize += HBondPool*sizeof(HBond);
  2016.     ptr = (HBond __far *)_fmalloc( HBondPool*sizeof(HBond) );
  2017.     if( !ptr ) FatalDataError("Memory allocation failed");
  2018.     RegisterAlloc( ptr );
  2019.     for( i=1; i<HBondPool; i++ )
  2020.     {   ptr->hnext = FreeHBond;
  2021.         FreeHBond = ptr++;
  2022.     } 
  2023.     } else FreeHBond = ptr->hnext;
  2024.  
  2025.     if( (offset>=-128) && (offset<127) )
  2026.     {   ptr->offset = (Char)offset;
  2027.     } else ptr->offset = 0;
  2028.  
  2029.     flag = ZoneBoth? src->flag&dst->flag : src->flag|dst->flag;
  2030.     ptr->flag = flag & SelectFlag;
  2031.  
  2032.     ptr->src = src;
  2033.     ptr->dst = dst;
  2034.     ptr->srcCA = srcCA;
  2035.     ptr->dstCA = dstCA;
  2036.     ptr->energy = energy;
  2037.     ptr->col = 0;
  2038.  
  2039.     *CurHBond = ptr;
  2040.     ptr->hnext = (void __far*)0;
  2041.     CurHBond = &ptr->hnext;
  2042.     InfoHBondCount++;
  2043. }
  2044.  
  2045. /* Coupling constant for Electrostatic Energy   */
  2046. /* QConst = -332 * 0.42 * 0.2 * 1000.0 [*250.0] */
  2047. #define QConst (-6972000.0)
  2048. #define MaxHDist ((Long)2250*2250)
  2049. #define MinHDist ((Long)125*125)
  2050.  
  2051.  
  2052. /* Protein Donor Atom Coordinates */
  2053. static int hxorg,hyorg,hzorg;
  2054. static int nxorg,nyorg,nzorg;
  2055. static Atom __far *best1CA;
  2056. static Atom __far *best2CA;
  2057. static Atom __far *best1;
  2058. static Atom __far *best2;
  2059. static Atom __far *optr;
  2060. static int res1,res2;
  2061. static int off1,off2;
  2062.  
  2063.  
  2064. static int CalculateBondEnergy( group )
  2065.     Group __far *group;
  2066. {
  2067.     register double dho,dhc;
  2068.     register double dnc,dno;
  2069.  
  2070.     register Atom __far *cptr;
  2071.     register Long dx,dy,dz;
  2072.     register Long dist;
  2073.     register int result;
  2074.  
  2075.     if( !(cptr=FindGroupAtom(group,2)) )  return(0);
  2076.     if( !(optr=FindGroupAtom(group,3)) )  return(0);
  2077.  
  2078.     dx = hxorg-optr->xorg;  
  2079.     dy = hyorg-optr->yorg;  
  2080.     dz = hzorg-optr->zorg;
  2081.     dist = dx*dx+dy*dy+dz*dz;
  2082.     if( dist < MinHDist ) 
  2083.     return( -9900 );
  2084.     dho = sqrt((double)dist);
  2085.  
  2086.     dx = hxorg-cptr->xorg;  
  2087.     dy = hyorg-cptr->yorg;  
  2088.     dz = hzorg-cptr->zorg;
  2089.     dist = dx*dx+dy*dy+dz*dz;
  2090.     if( dist < MinHDist ) 
  2091.     return( -9900 );
  2092.     dhc = sqrt((double)dist);
  2093.  
  2094.     dx = nxorg-cptr->xorg;  
  2095.     dy = nyorg-cptr->yorg;  
  2096.     dz = nzorg-cptr->zorg;
  2097.     dist = dx*dx+dy*dy+dz*dz;
  2098.     if( dist < MinHDist ) 
  2099.     return( -9900 );
  2100.     dnc = sqrt((double)dist);
  2101.  
  2102.     dx = nxorg-optr->xorg;  
  2103.     dy = nyorg-optr->yorg;  
  2104.     dz = nzorg-optr->zorg;
  2105.     dist = dx*dx+dy*dy+dz*dz;
  2106.     if( dist < MinHDist ) 
  2107.     return( -9900 );
  2108.     dno = sqrt((double)dist);
  2109.  
  2110.     result = (int)(QConst/dho - QConst/dhc + QConst/dnc - QConst/dno);
  2111.  
  2112.     if( result<-9900 ) 
  2113.     {   return( -9900 );
  2114.     } else if( result>-500 ) 
  2115.     return( 0 );
  2116.  
  2117.     return( result );
  2118.  
  2119. }
  2120.  
  2121.  
  2122. void CalcHydrogenBonds()
  2123. {
  2124.     register int energy, offset, refno;
  2125.     register Chain __far *chn1, __far *chn2;
  2126.     register Group __far *group1;
  2127.     register Group __far *group2;
  2128.     register Group __far *best;
  2129.     register Atom __far *ca1;
  2130.     register Atom __far *ca2;
  2131.     register Atom __far *pc1;
  2132.     register Atom __far *po1;
  2133.     register Atom __far *n1;
  2134.     register Long max,dist;
  2135.     register int pos1,pos2;
  2136.     register int dx,dy,dz;
  2137.     register double dco;
  2138.     char buffer[40];
  2139.  
  2140.     if( !Database ) return;
  2141.     ReclaimHBonds( CurMolecule->hlist );
  2142.     CurMolecule->hlist = (void __far*)0;
  2143.     CurHBond = &CurMolecule->hlist;
  2144.     InfoHBondCount = 0;
  2145.  
  2146.     if( MainAtomCount > 10000 )
  2147.     {   if( CommandActive )
  2148.         WriteChar('\n');
  2149.     WriteString("Please wait... ");
  2150.     CommandActive=True;
  2151.     }
  2152.  
  2153.     for(chn1=Database->clist;chn1;chn1=chn1->cnext)
  2154.     {   if( !chn1->glist ) continue;
  2155.  
  2156.     if( IsProtein(chn1->glist->refno) )
  2157.     {   pc1 = po1 = (void __far*)0;
  2158.         pos1 = 0;
  2159.         for(group1=chn1->glist;group1;group1=group1->gnext)
  2160.         {   pos1++;
  2161.         if( pc1 && po1 )
  2162.         {   dx = (int)(pc1->xorg - po1->xorg);
  2163.             dy = (int)(pc1->yorg - po1->yorg);
  2164.             dz = (int)(pc1->zorg - po1->zorg);
  2165.         } else
  2166.         {   pc1 = FindGroupAtom(group1,2);
  2167.             po1 = FindGroupAtom(group1,3);
  2168.             continue;
  2169.         }
  2170.  
  2171.         pc1 = FindGroupAtom(group1,2);
  2172.         po1 = FindGroupAtom(group1,3);
  2173.  
  2174.         if( !IsAmino(group1->refno) || IsProline(group1->refno) )
  2175.             continue;
  2176.  
  2177.         if( !(ca1=FindGroupAtom(group1,1)) ) continue;
  2178.         if( !(n1=FindGroupAtom(group1,0)) )  continue;
  2179.  
  2180.         dist = (Long)dx*dx + (Long)dy*dy + (Long)dz*dz;
  2181.         dco = sqrt( (double)dist )/250.0;
  2182.  
  2183.         nxorg = (int)n1->xorg;   hxorg = nxorg + (int)(dx/dco);
  2184.         nyorg = (int)n1->yorg;   hyorg = nyorg + (int)(dy/dco);
  2185.         nzorg = (int)n1->zorg;   hzorg = nzorg + (int)(dz/dco);
  2186.         res1 = res2 = 0;
  2187.  
  2188.         /* Only Hydrogen Bond within a single chain!       */
  2189.         /* for(chn2=Database->clist;chn2;chn2=chn2->cnext) */
  2190.  
  2191.         chn2 = chn1;
  2192.         {   /* Only consider non-empty peptide chains! */
  2193.             if( !chn2->glist || !IsProtein(chn2->glist->refno) )
  2194.             continue;
  2195.  
  2196.             pos2 = 0;
  2197.             for(group2=chn2->glist;group2;group2=group2->gnext)
  2198.             {   pos2++;
  2199.             if( (group2==group1) || (group2->gnext==group1) )
  2200.                 continue;
  2201.  
  2202.             if( !IsAmino(group2->refno) ) 
  2203.                 continue;
  2204.             if( !(ca2=FindGroupAtom(group2,1)) ) 
  2205.                 continue;
  2206.  
  2207.             dx = (int)(ca1->xorg-ca2->xorg);
  2208.             if( (dist=(Long)dx*dx) > MaxHDist )
  2209.                 continue;
  2210.  
  2211.             dy = (int)(ca1->yorg-ca2->yorg);
  2212.             if( (dist+=(Long)dy*dy) > MaxHDist )
  2213.                 continue;
  2214.  
  2215.             dz = (int)(ca1->zorg-ca2->zorg);
  2216.             if( (dist+=(Long)dz*dz) > MaxHDist )
  2217.                 continue;
  2218.  
  2219.             if( (energy = CalculateBondEnergy(group2)) )
  2220.             {   if( chn1 == chn2 )
  2221.                 {   offset = pos1 - pos2;
  2222.                 } else offset = 0;
  2223.  
  2224.                 if( energy<res1 )
  2225.                 {   best2CA = best1CA;  best1CA = ca2;
  2226.                 best2 = best1;      best1 = optr;
  2227.                 res2 = res1;        res1 = energy;
  2228.                 off2 = off1;        off1 = offset;
  2229.                 } else if( energy<res2 )
  2230.                 {   best2CA = ca2;
  2231.                 best2 = optr;
  2232.                 res2 = energy;
  2233.                 off2 = offset;
  2234.                 }
  2235.             }
  2236.             }  /* group2 */
  2237.         }      /* chn2 */
  2238.  
  2239.         if( res1 ) 
  2240.         {   if( res2 ) 
  2241.             CreateHydrogenBond(ca1,best2CA,n1,best2,res2,off2);
  2242.             CreateHydrogenBond(ca1,best1CA,n1,best1,res1,off1);
  2243.         }
  2244.         }
  2245.     } else if( IsNucleo(chn1->glist->refno) )
  2246.         for(group1=chn1->glist;group1;group1=group1->gnext)
  2247.         {   if( !IsPurine(group1->refno) ) continue;
  2248.         /* Find N1 of Purine Group */
  2249.         if( !(n1=FindGroupAtom(group1,21)) )
  2250.             continue;
  2251.  
  2252.         /* Maximum N1-N3 distance 5A */
  2253.         refno = NucleicCompl(group1->refno);
  2254.         max = (Long)1250*1250;
  2255.         best = (void __far*)0;
  2256.  
  2257.         for(chn2=Database->clist;chn2;chn2=chn2->cnext)
  2258.         {   /* Only consider non-empty peptide chains! */
  2259.             if( (chn1==chn2) || !chn2->glist || 
  2260.             !IsNucleo(chn2->glist->refno) )
  2261.             continue;
  2262.  
  2263.             for(group2=chn2->glist;group2;group2=group2->gnext)
  2264.             if( group2->refno == (Byte)refno )
  2265.             {   /* Find N3 of Pyramidine Group */
  2266.                 if( !(ca1=FindGroupAtom(group2,23)) )
  2267.                 continue;
  2268.  
  2269.                 dx = (int)(ca1->xorg - n1->xorg);
  2270.                 if( (dist=(Long)dx*dx) >= max ) 
  2271.                 continue;
  2272.  
  2273.                 dy = (int)(ca1->yorg - n1->yorg);
  2274.                 if( (dist+=(Long)dy*dy) >= max ) 
  2275.                 continue;
  2276.  
  2277.                 dz = (int)(ca1->zorg - n1->zorg);
  2278.                 if( (dist+=(Long)dz*dz) >= max )
  2279.                 continue;
  2280.  
  2281.                 best1 = ca1;
  2282.                 best = group2;
  2283.                 max = dist;
  2284.             }
  2285.         }
  2286.  
  2287.         if( best )
  2288.         {   /* Find the sugar phosphorous atoms */
  2289.             ca1 = FindGroupAtom( group1, 7 );
  2290.             ca2 = FindGroupAtom( best, 7 );
  2291.  
  2292.             CreateHydrogenBond( ca1, ca2, n1, best1, 0, 0 );
  2293.             if( IsGuanine(group1->refno) )
  2294.             {   /* Guanine-Cytosine */
  2295.             if( (ca1=FindGroupAtom(group1,22)) &&  /* G.N2 */
  2296.                 (ca2=FindGroupAtom(best,26)) )     /* C.O2 */
  2297.                 CreateHydrogenBond( (void __far*)0, (void __far*)0,
  2298.                         ca1, ca2, 0, 0 );
  2299.  
  2300.             if( (ca1=FindGroupAtom(group1,28)) &&  /* G.O6 */
  2301.                 (ca2=FindGroupAtom(best,24)) )     /* C.N4 */
  2302.                 CreateHydrogenBond( (void __far*)0, (void __far*)0,
  2303.                         ca1, ca2, 0, 0 );
  2304.  
  2305.             } else /* Adenine-Thymine */
  2306.             if( (ca1=FindGroupAtom(group1,25)) &&  /* A.N6 */
  2307.                 (ca2=FindGroupAtom(best,27)) )     /* T.O4 */
  2308.                 CreateHydrogenBond( (void __far*)0, (void __far*)0,
  2309.                         ca1, ca2, 0, 0 );
  2310.         }
  2311.         }
  2312.     }
  2313.  
  2314.     if( CommandActive )
  2315.     WriteChar('\n');
  2316.     CommandActive=False;
  2317.     
  2318.     sprintf(buffer,"Number of H-Bonds ... %d\n",InfoHBondCount);
  2319.     WriteString(buffer);
  2320. }
  2321.  
  2322.  
  2323. static int IsHBonded( src, dst, ptr )
  2324.     Atom __far *src, __far *dst;
  2325.     HBond __far *ptr;
  2326. {
  2327.     while( ptr && (ptr->srcCA==src) )
  2328.     if( ptr->dstCA == dst )
  2329.     {   return( True );
  2330.     } else ptr=ptr->hnext;
  2331.     return( False );
  2332. }
  2333.  
  2334.  
  2335. static void FindAlphaHelix( pitch, flag )
  2336.     int pitch, flag;
  2337. {
  2338.     register HBond __far *hbond;
  2339.     register Chain __far *chain;
  2340.     register Group __far *group;
  2341.     register Group __far *first;
  2342.     register Group __far *ptr;
  2343.     register Atom __far *srcCA;
  2344.     register Atom __far *dstCA;
  2345.     register int res,dist,prev;
  2346.  
  2347.     /* Protein chains only! */
  2348.     hbond = Database->hlist;
  2349.     for( chain=Database->clist; chain; chain=chain->cnext )
  2350.     if( (first=chain->glist) && IsProtein(first->refno) )
  2351.     {   prev = False; dist = 0;
  2352.     for( group=chain->glist; hbond && group; group=group->gnext )
  2353.     {   if( IsAmino(group->refno) && (srcCA=FindGroupAtom(group,1)) )
  2354.         {   if( dist==pitch )
  2355.         {   res = False;
  2356.             dstCA=FindGroupAtom(first,1);
  2357.  
  2358.             while( hbond && hbond->srcCA == srcCA )
  2359.             {   if( hbond->dstCA==dstCA ) res=True;
  2360.             hbond = hbond->hnext;
  2361.             }
  2362.  
  2363.             if( res )
  2364.             {   if( prev )
  2365.             {   if( !(first->struc&HelixFlag) ) 
  2366.                 InfoHelixCount++;
  2367.  
  2368.                 ptr = first;
  2369.                 do {
  2370.                 ptr->struc |= flag;
  2371.                 ptr = ptr->gnext;
  2372.                 } while( ptr != group );
  2373.             } else prev = True;
  2374.             } else prev = False;
  2375.         } else while( hbond && hbond->srcCA==srcCA )
  2376.             hbond = hbond->hnext;
  2377.         } else prev = False;
  2378.  
  2379.         if( group->struc&HelixFlag )
  2380.         {   first = group; prev = False; dist = 1;
  2381.         } else if( dist==pitch )
  2382.         {   first = first->gnext;
  2383.         } else dist++;
  2384.     }
  2385.     } else if( first && IsNucleo(first->refno) )
  2386.     while( hbond && !IsAminoBackbone(hbond->src->refno) )
  2387.         hbond = hbond->hnext;
  2388. }
  2389.  
  2390.  
  2391. static Atom __far *cprevi, __far *ccurri, __far *cnexti;
  2392. static HBond __far *hcurri, __far *hnexti;
  2393. static Group __far *curri, __far *nexti;
  2394.  
  2395.  
  2396.  
  2397. static void TestLadder( chain )
  2398.     Chain __far *chain;
  2399. {
  2400.     register Atom __far *cprevj, __far *ccurrj, __far *cnextj;
  2401.     register HBond __far *hcurrj, __far *hnextj;
  2402.     register Group __far *currj, __far *nextj;
  2403.     register int count, result, found;
  2404.  
  2405.     /* Already part of atleast one ladder */
  2406.     found = curri->flag & SheetFlag;
  2407.     nextj = nexti->gnext;
  2408.  
  2409.     hnextj = hnexti;
  2410.     while( hnextj && hnextj->srcCA==cnexti )
  2411.     hnextj = hnextj->hnext;
  2412.  
  2413.     while( True )
  2414.     {   if( nextj )
  2415.         if( IsProtein(chain->glist->refno) )
  2416.         {   count = 1;
  2417.         do {
  2418.             cnextj = FindGroupAtom(nextj,1);
  2419.             if( count == 3 )
  2420.             {   if( IsHBonded(cnexti,ccurrj,hnexti) &&
  2421.                 IsHBonded(ccurrj,cprevi,hcurrj) )
  2422.             {   result = ParaLadder;
  2423.             } else if( IsHBonded(cnextj,ccurri,hnextj) &&
  2424.                    IsHBonded(ccurri,cprevj,hcurri) )
  2425.             {   result = ParaLadder;
  2426.             } else if( IsHBonded(cnexti,cprevj,hnexti) &&
  2427.                    IsHBonded(cnextj,cprevi,hnextj) )
  2428.             {   result = AntiLadder;
  2429.             } else if( IsHBonded(ccurrj,ccurri,hcurrj) &&
  2430.                    IsHBonded(ccurri,ccurrj,hcurri) )
  2431.             {   result = AntiLadder;
  2432.             } else result = NoLadder;
  2433.  
  2434.             if( result )
  2435.             {   curri->struc |= SheetFlag;
  2436.                 currj->struc |= SheetFlag;
  2437.                 if( found ) return;
  2438.                 found = True;
  2439.             }
  2440.             } else count++;
  2441.  
  2442.             cprevj = ccurrj; ccurrj = cnextj; 
  2443.             currj = nextj;   hcurrj = hnextj;
  2444.  
  2445.             while( hnextj && hnextj->srcCA==cnextj )
  2446.             hnextj = hnextj->hnext;
  2447.         } while( (nextj = nextj->gnext) );
  2448.  
  2449.         } else if( IsNucleo(chain->glist->refno) )
  2450.         while( hnextj && !IsAminoBackbone(hnextj->src->refno) )
  2451.             hnextj = hnextj->hnext;
  2452.  
  2453.     if( (chain = chain->cnext) ) 
  2454.     {   nextj = chain->glist;
  2455.     } else return;
  2456.     }
  2457. }
  2458.  
  2459.  
  2460. static void FindBetaSheets()
  2461. {
  2462.     register Chain __far *chain;
  2463.     register int ladder;
  2464.     register int count;
  2465.  
  2466.     hnexti = Database->hlist;
  2467.     for( chain=Database->clist; chain; chain=chain->cnext )
  2468.     if( (nexti = chain->glist) )
  2469.         if( IsProtein(nexti->refno) )
  2470.         {   count = 1;
  2471.         ladder = False;
  2472.         do {
  2473.             cnexti = FindGroupAtom(nexti,1);
  2474.  
  2475.             if( count == 3 )
  2476.             {   TestLadder( chain );
  2477.             if( curri->struc & SheetFlag )
  2478.             {   if( !ladder )
  2479.                 {   InfoLadderCount++;
  2480.                 ladder = True;
  2481.                 }
  2482.             } else ladder = False;
  2483.             } else count++;
  2484.  
  2485.             cprevi = ccurri; ccurri = cnexti; 
  2486.             curri = nexti;   hcurri = hnexti;
  2487.             while( hnexti && hnexti->srcCA==cnexti )
  2488.             hnexti = hnexti->hnext;
  2489.         } while( (nexti = nexti->gnext) );
  2490.  
  2491.         } else if( IsNucleo(nexti->refno) )
  2492.         while( hnexti && !IsAminoBackbone(hnexti->src->refno) )
  2493.             hnexti = hnexti->hnext;
  2494. }
  2495.  
  2496.  
  2497. static void FindTurnStructure()
  2498. {
  2499.     static Atom __far *aptr[5];
  2500.     register Chain __far *chain;
  2501.     register Group __far *group;
  2502.     register Group __far *prev;
  2503.     register Atom __far *ptr;
  2504.     register Long ux,uy,uz,mu;
  2505.     register Long vx,vy,vz,mv;
  2506.     register int i,found,len;
  2507.     register Real CosKappa;
  2508.  
  2509.     for( chain=Database->clist; chain; chain=chain->cnext )
  2510.     if( chain->glist && IsProtein(chain->glist->refno) )
  2511.     {   len = 0;  found = False;
  2512.         for( group=chain->glist; group; group=group->gnext )
  2513.         {    ptr = FindGroupAtom(group,1);
  2514.          if( ptr && (ptr->flag&BreakFlag) )
  2515.          {   found = False;
  2516.              len = 0;
  2517.          } else if( len==5 )
  2518.          {   for( i=0; i<4; i++ )
  2519.              aptr[i] = aptr[i+1];
  2520.              len = 4;
  2521.          } else if( len==2 )
  2522.              prev = group;
  2523.  
  2524.          aptr[len++] = ptr;
  2525.          if( len==5 ) 
  2526.          {   if( !(prev->struc&(HelixFlag|SheetFlag)) &&
  2527.              aptr[0] && aptr[2] && aptr[4] )
  2528.              {   ux = aptr[2]->xorg - aptr[0]->xorg;
  2529.              uy = aptr[2]->yorg - aptr[0]->yorg;
  2530.              uz = aptr[2]->zorg - aptr[0]->zorg;
  2531.  
  2532.              vx = aptr[4]->xorg - aptr[2]->xorg;
  2533.              vy = aptr[4]->yorg - aptr[2]->yorg;
  2534.              vz = aptr[4]->zorg - aptr[2]->zorg;
  2535.  
  2536.              mu = ux*ux + uy*uy + uz*uz;
  2537.              mv = vx*vx + vz*vz + vy*vy;
  2538.              if( mu && mv )
  2539.              {   CosKappa = (Real)(ux*vx + uy*vy + uz*vz);
  2540.                  CosKappa /= sqrt( (Real)mu*mv );
  2541.                  if( CosKappa<Cos70Deg )
  2542.                  {   if( !found )
  2543.                      InfoTurnCount++;
  2544.                  prev->struc |= TurnFlag;
  2545.                  }
  2546.              }
  2547.              }
  2548.              found = prev->struc&TurnFlag;
  2549.              prev = prev->gnext;
  2550.          } /* len==5 */
  2551.         }
  2552.     }
  2553. }
  2554.  
  2555. void DetermineStructure()
  2556. {
  2557.     register Chain __far *chain;
  2558.     register Group __far *group;
  2559.     char buffer[40];
  2560.  
  2561.     if( !Database )
  2562.     return;
  2563.  
  2564.     if( InfoHBondCount<0 )
  2565.     CalcHydrogenBonds();
  2566.  
  2567.     if( InfoHelixCount>=0 )
  2568.     for( chain=Database->clist; chain; chain=chain->cnext )
  2569.         for( group=chain->glist; group; group=group->gnext )
  2570.         group->struc = 0;
  2571.  
  2572.     InfoStrucSource = SourceCalc;
  2573.     InfoLadderCount = 0;
  2574.     InfoHelixCount = 0;
  2575.     InfoTurnCount = 0;
  2576.  
  2577.     if( InfoHBondCount )
  2578.     {   FindAlphaHelix(4,Helix4Flag);
  2579.     FindBetaSheets();
  2580.     FindAlphaHelix(3,Helix3Flag);
  2581.     FindAlphaHelix(5,Helix5Flag);
  2582.     FindTurnStructure();
  2583.     }
  2584.  
  2585.     if( CommandActive )
  2586.     WriteChar('\n');
  2587.     CommandActive=False;
  2588.  
  2589.     sprintf(buffer,"Number of Helices ... %d\n",InfoHelixCount);
  2590.     WriteString(buffer);
  2591.     sprintf(buffer,"Number of Strands ... %d\n",InfoLadderCount);
  2592.     WriteString(buffer);
  2593.     sprintf(buffer,"Number of Turns ..... %d\n",InfoTurnCount);
  2594.     WriteString(buffer);
  2595. }
  2596.  
  2597.  
  2598. void RenumberMolecule( start )
  2599.     int start;
  2600. {
  2601.     register Chain __far *chain;
  2602.     register Group __far *group;
  2603.     register int hinit, minit;
  2604.     register int resno;
  2605.  
  2606.     if( !Database )
  2607.     return;
  2608.  
  2609.     hinit = minit = False;
  2610.     for( chain=Database->clist; chain; chain=chain->cnext )
  2611.     {   resno = start;
  2612.     for( group=chain->glist; group; group=group->gnext )
  2613.     {   if( group->alist->flag & HeteroFlag )
  2614.         {   if( hinit )
  2615.         {   if( resno > MaxHetaRes )
  2616.             {   MaxHetaRes = resno;
  2617.             } else if( resno < MinHetaRes )
  2618.             MinHetaRes = resno;
  2619.         } else MinHetaRes = MaxHetaRes = resno;
  2620.         hinit = True;
  2621.         } else
  2622.         {   if( minit )
  2623.         {   if( resno > MaxMainRes )
  2624.             {   MaxMainRes = resno;
  2625.             } else if( resno < MinMainRes )
  2626.             MinMainRes = resno;
  2627.         } else MinMainRes = MaxMainRes = resno;
  2628.         minit = True;
  2629.         }
  2630.         group->serno = resno++;
  2631.     }
  2632.     }
  2633. }
  2634.  
  2635.  
  2636. static void ReclaimAtoms( ptr )
  2637.     Atom __far *ptr;
  2638. {
  2639.     register Atom __far *temp;
  2640.  
  2641.     if( (temp = ptr) )
  2642.     {   while( temp->anext )
  2643.         temp=temp->anext;
  2644.     temp->anext = FreeAtom;
  2645.     FreeAtom = ptr;
  2646.     }
  2647. }
  2648.  
  2649. static void ResetDatabase()
  2650. {
  2651.     
  2652.     Database = CurMolecule = (void __far*)0;
  2653.     MainGroupCount = HetaGroupCount = 0;
  2654.     InfoChainCount = HetaAtomCount = 0;
  2655.     MainAtomCount = InfoBondCount = 0;  
  2656.     SelectCount = 0;
  2657.  
  2658.     InfoStrucSource = SourceNone;
  2659.     InfoSSBondCount = InfoHBondCount = -1;
  2660.     InfoHelixCount = InfoLadderCount = -1;
  2661.     InfoTurnCount = -1;
  2662.  
  2663.     CurGroup = (void __far*)0;
  2664.     CurChain = (void __far*)0;
  2665.     CurAtom = (void __far*)0;
  2666.  
  2667.     MinX = MinY = MinZ = 0;
  2668.     MaxX = MaxY = MaxZ = 0;
  2669.  
  2670.     MinMainTemp = MaxMainTemp = 0;
  2671.     MinHetaTemp = MaxHetaTemp = 0;
  2672.     MinMainRes = MaxMainRes = 0;
  2673.     MinHetaRes = MaxHetaRes = 0;
  2674.  
  2675.     *InfoMoleculeName = 0;
  2676.     *InfoClassification = 0;
  2677.     *InfoIdentCode = 0;
  2678.     *InfoSpaceGroup = 0;
  2679.     *InfoFileName = 0;
  2680.  
  2681.     VoxelsClean = False;
  2682.     HMinMaxFlag = False;
  2683.     MMinMaxFlag = False;
  2684.     HasHydrogen = False;
  2685.     ElemNo = MINELEM;
  2686.     ResNo = MINRES;
  2687.     MaskCount = 0;
  2688. }
  2689.  
  2690.  
  2691. void DestroyDatabase()
  2692. {
  2693.     register void __far *temp;
  2694.     register Group __far *gptr;
  2695.  
  2696.     if( Database )
  2697.     {   ReclaimHBonds( Database->slist );
  2698.     ReclaimHBonds( Database->hlist );
  2699.     ReclaimBonds( Database->blist );
  2700.  
  2701.     while( Database->clist )
  2702.     {   ReclaimBonds(Database->clist->blist);
  2703.         if( (gptr = Database->clist->glist) )
  2704.         {   ReclaimAtoms(gptr->alist);
  2705.         while( gptr->gnext )
  2706.         {   gptr = gptr->gnext;
  2707.             ReclaimAtoms(gptr->alist);
  2708.         }
  2709.         gptr->gnext = FreeGroup;
  2710.         FreeGroup = Database->clist->glist;
  2711.         }
  2712.         temp = Database->clist->cnext;
  2713.         Database->clist->cnext = FreeChain;
  2714.         FreeChain = Database->clist;
  2715.         Database->clist = temp;
  2716.     }
  2717.  
  2718.     FreeMolecule = Database;
  2719.     Database = (void __far*)0;
  2720.     }
  2721.     ResetDatabase();
  2722. }
  2723.  
  2724.  
  2725. void PurgeDatabase()
  2726. {
  2727. #ifdef APPLEMAC
  2728.     register AllocRef *ptr;
  2729.     register AllocRef *tmp;
  2730.     register int i;
  2731.     
  2732.     /* Avoid Memory Leaks */
  2733.     for( ptr=AllocList; ptr; ptr=tmp )
  2734.     {   for( i=0; i<ptr->count; i++ )
  2735.         _ffree( ptr->data[i] );
  2736.     tmp = ptr->next;
  2737.     _ffree( ptr );
  2738.     }
  2739. #endif
  2740. }
  2741.  
  2742.  
  2743. void InitialiseDatabase()
  2744. {
  2745.     FreeMolecule = (void __far*)0;
  2746.     FreeHBond = (void __far*)0;
  2747.     FreeChain = (void __far*)0;
  2748.     FreeGroup = (void __far*)0;
  2749.     FreeAtom = (void __far*)0;
  2750.     FreeBond = (void __far*)0;
  2751.  
  2752. #ifdef APPLEMAC
  2753.     AllocList = (void*)0;
  2754. #endif
  2755.  
  2756.     ResetDatabase();
  2757. }
  2758.  
  2759.