This is molecule.c in view mode; [Download] [Up]
/* molecule.c
* RasMol2 Molecular Graphics
* Roger Sayle, August 1995
* Version 2.6
*/
#include "rasmol.h"
#ifdef IBMPC
#include <windows.h>
#include <malloc.h>
#endif
#ifdef APPLEMAC
#include <Types.h>
#endif
#ifndef sun386
#include <stdlib.h>
#endif
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include <math.h>
#define MOLECULE
#include "molecule.h"
#include "command.h"
#include "abstree.h"
#include "transfor.h"
#include "render.h"
#define HBondPool 32
#define BondPool 32
#define AtomPool 32
#define NoLadder 0x00
#define ParaLadder 0x01
#define AntiLadder 0x02
#define Cos70Deg 0.34202014332567
#define MaxHBondDist ((Long)300*300)
#define MaxBondDist ((Long)475*475)
#define MinBondDist ((Long)100*100)
#define AbsMaxBondDist 600
#ifdef APPLEMAC
#define AllocSize 256
typedef struct _AllocRef {
struct _AllocRef *next;
void *data[AllocSize];
int count;
} AllocRef;
static AllocRef *AllocList;
#endif
typedef struct {
char name[4];
int code;
} SynonymTable;
#define RESSYNMAX 16
static SynonymTable ResSynonym[RESSYNMAX] = {
{ "ADE", 24 }, /* A : Adenosine */
{ "CPR", 11 }, /* PRO : Cis-proline */
{ "CSH", 17 }, /* CYS : Cystine */
{ "CSM", 17 }, /* CYS : Cystine */
{ "CYH", 17 }, /* CYS : Cystine */
{ "CYT", 25 }, /* C : Cytosine */
{ "D2O", 47 }, /* DOD : Heavy Water */
{ "GUA", 26 }, /* G : Guanosine */
{ "H2O", 46 }, /* HOH : Solvent */
{ "SOL", 46 }, /* HOH : Solvent */
{ "SUL", 48 }, /* SO4 : Sulphate */
{ "THY", 27 }, /* T : Thymidine */
{ "TIP", 46 }, /* HOH : Water */
{ "TRY", 20 }, /* TRP : Tryptophan */
{ "URI", 28 }, /* U : Uridine */
{ "WAT", 46 } /* HOH : Water */
};
static Molecule __far *FreeMolecule;
static HBond __far *FreeHBond;
static Chain __far *FreeChain;
static Group __far *FreeGroup;
static Atom __far *FreeAtom;
static Bond __far *FreeBond;
static IntCoord __far *IntPrev;
static HBond __far * __far *CurHBond;
static int MemSize;
/* Macros for commonly used loops */
#define ForEachAtom for(chain=Database->clist;chain;chain=chain->cnext) \
for(group=chain->glist;group;group=group->gnext) \
for(aptr=group->alist;aptr;aptr=aptr->anext)
#define ForEachBond for(bptr=Database->blist;bptr;bptr=bptr->bnext)
/* Forward Reference */
void DestroyDatabase();
#ifdef APPLEMAC
/* External RasMac Function Declaration! */
void SetFileInfo( char*, OSType, OSType, short );
#endif
static void FatalDataError(ptr)
char *ptr;
{
char buffer[80];
sprintf(buffer,"Database Error: %s!",ptr);
RasMolFatalExit(buffer);
}
void DescribeMolecule()
{
char buffer[40];
if( CommandActive )
WriteChar('\n');
CommandActive=False;
if( *InfoMoleculeName )
{ WriteString("Molecule name ....... ");
WriteString(InfoMoleculeName);
WriteChar('\n');
}
if( *InfoClassification )
{ WriteString("Classification ...... ");
WriteString(InfoClassification);
WriteChar('\n');
}
if( Database && (MainGroupCount>1) )
{ WriteString("Secondary Structure . ");
if( InfoStrucSrc==SourceNone )
{ WriteString("No Assignment\n");
} else if( InfoStrucSrc==SourcePDB )
{ WriteString("PDB Data Records\n");
} else WriteString("Calculated\n");
}
if( *InfoIdentCode )
{ WriteString("Brookhaven Code ..... ");
WriteString(InfoIdentCode);
WriteChar('\n');
}
if( InfoChainCount>1 )
{ sprintf(buffer,"Number of Chains .... %d\n",InfoChainCount);
WriteString(buffer);
}
sprintf(buffer,"Number of Groups .... %d",MainGroupCount);
WriteString(buffer);
if( HetaAtomCount )
{ sprintf(buffer," (%d)\n",HetaGroupCount);
WriteString(buffer);
} else WriteChar('\n');
sprintf(buffer,"Number of Atoms ..... %ld",(long)MainAtomCount);
WriteString(buffer);
if( HetaAtomCount )
{ sprintf(buffer," (%d)\n",HetaAtomCount);
WriteString(buffer);
} else WriteChar('\n');
if( InfoBondCount )
{ sprintf(buffer,"Number of Bonds ..... %ld\n",(long)InfoBondCount);
WriteString(buffer);
}
if( InfoSSBondCount != -1 )
{ WriteString("Number of Bridges ... ");
sprintf(buffer,"%d\n\n",InfoSSBondCount);
WriteString(buffer);
}
if( InfoHBondCount != -1 )
{ WriteString("Number of H-Bonds ... ");
sprintf(buffer,"%d\n",InfoHBondCount);
WriteString(buffer);
}
if( InfoHelixCount != -1 )
{ WriteString("Number of Helices ... ");
sprintf(buffer,"%d\n",InfoHelixCount);
WriteString(buffer);
WriteString("Number of Strands ... ");
sprintf(buffer,"%d\n",InfoLadderCount);
WriteString(buffer);
WriteString("Number of Turns ..... ");
sprintf(buffer,"%d\n",InfoTurnCount);
WriteString(buffer);
}
}
#ifdef APPLEMAC
/* Avoid System Memory Leaks! */
static void RegisterAlloc( data )
void *data;
{
register AllocRef *ptr;
if( !AllocList || (AllocList->count==AllocSize) )
{ ptr = (AllocRef *)_fmalloc( sizeof(AllocRef) );
if( !ptr ) FatalDataError("Memory allocation failed");
ptr->next = AllocList;
ptr->data[0] = data;
ptr->count = 1;
AllocList = ptr;
} else AllocList->data[AllocList->count++] = data;
}
#else
#define RegisterAlloc(x)
#endif
/*==================================*/
/* Group & Chain Handling Functions */
/*==================================*/
void CreateChain( ident )
int ident;
{
register Chain __far *prev;
if( !CurMolecule )
{ if( !(CurMolecule = FreeMolecule) )
{ MemSize += sizeof(Molecule);
CurMolecule = (Molecule __far *)_fmalloc(sizeof(Molecule));
if( !CurMolecule ) FatalDataError("Memory allocation failed");
RegisterAlloc( CurMolecule );
} else FreeMolecule = (void __far*)0;
CurChain = (void __far*)0;
CurMolecule->slist = (void __far*)0;
CurMolecule->hlist = (void __far*)0;
CurMolecule->blist = (void __far*)0;
CurMolecule->clist = (void __far*)0;
Database = CurMolecule;
}
/* Handle chain breaks! */
if( !(prev=CurChain) )
if( (prev=CurMolecule->clist) )
while( prev->cnext )
prev = prev->cnext;
if( !(CurChain = FreeChain) )
{ MemSize += sizeof(Chain);
CurChain = (Chain __far *)_fmalloc(sizeof(Chain));
if( !CurChain ) FatalDataError("Memory allocation failed");
RegisterAlloc( CurChain );
} else FreeChain = FreeChain->cnext;
if( prev )
{ prev->cnext = CurChain;
} else CurMolecule->clist = CurChain;
CurChain->cnext = (void __far*)0;
CurChain->ident = ident;
CurChain->model = NMRModel;
CurChain->glist = (void __far*)0;
CurChain->blist = (void __far*)0;
ConnectAtom = (void __far*)0;
CurGroup = (void __far*)0;
InfoChainCount++;
}
void CreateGroup( pool )
int pool;
{
register Group __far *ptr;
register int i;
if( !(ptr = FreeGroup) )
{ MemSize += pool*sizeof(Group);
ptr = (Group __far *)_fmalloc( pool*sizeof(Group) );
if( !ptr ) FatalDataError("Memory allocation failed");
RegisterAlloc( ptr );
for( i=1; i<pool; i++ )
{ ptr->gnext = FreeGroup;
FreeGroup = ptr++;
}
} else FreeGroup = ptr->gnext;
if( CurGroup )
{ ptr->gnext = CurGroup->gnext;
CurGroup->gnext = ptr;
} else
{ ptr->gnext = CurChain->glist;
CurChain->glist = ptr;
}
CurGroup = ptr;
CurAtom = (void __far*)0;
ptr->alist = (void __far*)0;
ptr->insert = ' ';
ptr->struc = 0;
ptr->flag = 0;
ptr->col1 = 0;
ptr->col2 = 0;
}
int FindResNo( ptr )
char *ptr;
{
register int hi,lo;
register int refno;
register int flag;
register int mid;
for( refno=0; refno<ResNo; refno++ )
if( !strncmp(Residue[refno],ptr,3) )
return( refno );
lo = 0;
hi = RESSYNMAX;
while( lo < hi )
{ mid = (hi+lo)>>1;
flag = strncmp(ResSynonym[mid].name,ptr,3);
if( !flag ) return( ResSynonym[mid].code );
/* Binary Search */
if( flag<0 )
{ lo = mid+1;
} else hi = mid;
}
if( ResNo++ == MAXRES )
FatalDataError("Too many new residues");
Residue[refno][0] = *ptr++;
Residue[refno][1] = *ptr++;
Residue[refno][2] = *ptr;
return( refno );
}
void ProcessGroup( heta )
int heta;
{
register int serno;
serno = CurGroup->serno;
if( IsSolvent(CurGroup->refno) )
heta = True;
if( heta )
{ HetaGroupCount++;
if( HMinMaxFlag )
{ if( serno > MaxHetaRes )
{ MaxHetaRes = serno;
} else if( serno < MinHetaRes )
MinHetaRes = serno;
} else MinHetaRes = MaxHetaRes = serno;
} else
{ MainGroupCount++;
if( MMinMaxFlag )
{ if( serno > MaxMainRes )
{ MaxMainRes = serno;
} else if( serno < MinMainRes )
MinMainRes = serno;
} else MinMainRes = MaxMainRes = serno;
}
}
void CreateMolGroup()
{
strcpy(InfoFileName,DataFileName);
CreateChain( ' ' );
CreateGroup( 1 );
CurGroup->refno = FindResNo( "MOL" );
CurGroup->serno = 1;
MinMainRes = MaxMainRes = 1;
MinHetaRes = MaxHetaRes = 0;
MainGroupCount = 1;
}
/*=========================*/
/* Atom Handling Functions */
/*=========================*/
Atom __far *CreateAtom()
{
register Atom __far *ptr;
register int i;
if( !(ptr = FreeAtom) )
{ MemSize += AtomPool*sizeof(Atom);
ptr = (Atom __far *)_fmalloc( AtomPool*sizeof(Atom) );
if( !ptr ) FatalDataError("Memory allocation failed");
RegisterAlloc( ptr );
for( i=1; i<AtomPool; i++ )
{ ptr->anext = FreeAtom;
FreeAtom = ptr++;
}
} else FreeAtom = ptr->anext;
if( CurAtom )
{ ptr->anext = CurAtom->anext;
CurAtom->anext = ptr;
} else
{ ptr->anext = CurGroup->alist;
CurGroup->alist = ptr;
}
CurAtom = ptr;
SelectCount++;
ptr->flag = SelectFlag | NonBondFlag;
ptr->label = (void*)0;
ptr->radius = 375;
ptr->altl = ' ';
ptr->mbox = 0;
ptr->col = 0;
return( ptr );
}
void ProcessAtom( ptr )
Atom __far *ptr;
{
ptr->elemno = GetElemNumber(CurGroup,ptr);
if( ptr->elemno == 1 )
{ ptr->flag |= HydrogenFlag;
HasHydrogen = True;
}
if( !IsSolvent(CurGroup->refno) )
{ if( !(ptr->flag&(HydrogenFlag|HeteroFlag)) )
ptr->flag |= NormAtomFlag;
} else ptr->flag |= HeteroFlag;
#ifdef INVERT
ptr->yorg = -ptr->yorg;
#endif
if( HMinMaxFlag || MMinMaxFlag )
{ if( ptr->xorg < MinX )
{ MinX = ptr->xorg;
} else if( ptr->xorg > MaxX )
MaxX = ptr->xorg;
if( ptr->yorg < MinY )
{ MinY = ptr->yorg;
} else if( ptr->yorg > MaxY )
MaxY = ptr->yorg;
if( ptr->zorg < MinZ )
{ MinZ = ptr->zorg;
} else if( ptr->zorg > MaxZ )
MaxZ = ptr->zorg;
} else
{ MinX = MaxX = ptr->xorg;
MinY = MaxY = ptr->yorg;
MinZ = MaxZ = ptr->zorg;
}
if( ptr->flag & HeteroFlag )
{ if( HMinMaxFlag )
{ if( ptr->temp < MinHetaTemp )
{ MinHetaTemp = ptr->temp;
} else if( ptr->temp > MaxHetaTemp )
MaxHetaTemp = ptr->temp;
} else MinHetaTemp = MaxHetaTemp = ptr->temp;
HMinMaxFlag = True;
HetaAtomCount++;
} else
{ if( MMinMaxFlag )
{ if( ptr->temp < MinMainTemp )
{ MinMainTemp = ptr->temp;
} else if( ptr->temp > MaxMainTemp )
MaxMainTemp = ptr->temp;
} else MinMainTemp = MaxMainTemp = ptr->temp;
MMinMaxFlag = True;
MainAtomCount++;
}
}
Atom __far *FindGroupAtom( group, n )
Group __far *group; Byte n;
{
register Atom __far *ptr;
for( ptr=group->alist; ptr; ptr=ptr->anext )
if( ptr->refno == n ) return( ptr );
return( (Atom __far*)0 );
}
int NewAtomType( ptr )
char *ptr;
{
register int refno;
register int i;
for( refno=0; refno<ElemNo; refno++ )
if( !strncmp(ElemDesc[refno],ptr,4) )
return(refno);
if( ElemNo++ == MAXELEM )
FatalDataError("Too many new atom types");
for( i=0; i<4; i++ )
ElemDesc[refno][i] = ptr[i];
return( refno );
}
int SimpleAtomType( type )
char *type;
{
char name[4];
name[2] = name[3] = ' ';
if( type[1] && (type[1]!=' ') )
{ name[0] = ToUpper(type[0]);
name[1] = ToUpper(type[1]);
} else
{ name[1] = ToUpper(type[0]);
name[0] = ' ';
}
return( NewAtomType(name) );
}
int ComplexAtomType( ptr )
char *ptr;
{
static char name[4];
register int i;
if( isdigit(ptr[1]) )
{ name[1] = ToUpper(*ptr);
name[2] = name[3] = ' ';
name[0] = ' ';
} else for( i=0; i<4; i++ )
name[i] = ToUpper(ptr[i]);
/* Handle Unconventional Naming */
if( IsProtein(CurGroup->refno) )
{ if( name[0]=='H' )
{ name[0]=' ';
name[1]='H';
}
} else if( IsNucleo(CurGroup->refno) )
{ if( name[3]=='\'' )
name[3] = '*';
if( (name[1]=='O') && (name[2]=='P') )
{ if( !strncmp(name," OP1",4) ||
!strncmp(name,"1OP ",4) )
return( 8 );
if( !strncmp(name," OP2",4) ||
!strncmp(name,"2OP ",4) )
return( 9 );
}
}
return( NewAtomType(name) );
}
/*===============================*/
/* Z-Matrix Conversion Functions */
/*===============================*/
#ifdef FUNCPROTO
static IntCoord __far* GetInternalCoord( int );
#endif
void InitInternalCoords()
{
IntList = (IntCoord __far*)0;
IntPrev = (IntCoord __far*)0;
}
IntCoord __far* AllocInternalCoord()
{
register IntCoord __far *ptr;
ptr = (IntCoord __far*)_fmalloc(sizeof(IntCoord));
if( !ptr ) FatalDataError("Memory allocation failed");
ptr->inext = (IntCoord __far*)0;
if( IntPrev )
{ IntPrev->inext = ptr;
} else IntList = ptr;
IntPrev = ptr;
return( ptr );
}
static IntCoord __far* GetInternalCoord( index )
int index;
{
register IntCoord __far *ptr;
ptr = IntList;
while( (index>1) && ptr->inext )
{ ptr = ptr->inext;
index--;
}
return( ptr );
}
void FreeInternalCoords()
{
register IntCoord __far *ptr;
while( (ptr = IntList) )
{ IntList = ptr->inext;
_ffree( ptr );
}
}
int ConvertInternal2Cartesian()
{
register IntCoord __far *ptr;
register IntCoord __far *na;
register IntCoord __far *nb;
register IntCoord __far *nc;
register double cosph,sinph,costh,sinth,coskh,sinkh;
register double cosa,sina,cosd,sind;
register double dist,angle,dihed;
register double xpd,ypd,zpd,xqd,yqd,zqd;
register double xa,ya,za,xb,yb,zb;
register double rbc,xyb,yza,temp;
register double xpa,ypa,zqa;
register double xd,yd,zd;
register int flag;
/* Atom #1 */
ptr = IntList;
ptr->dist = 0.0;
ptr->angle = 0.0;
ptr->dihed = 0.0;
if( !(ptr=ptr->inext) )
return( True );
/* Atom #2 */
ptr->angle = 0.0;
ptr->dihed = 0.0;
if( !(ptr=ptr->inext) )
return( True );
/* Atom #3 */
dist = ptr->dist;
angle = Deg2Rad*ptr->angle;
cosa = cos(angle);
sina = sin(angle);
if( ptr->na == 1 )
{ na = IntList;
ptr->dist = na->dist + cosa*dist;
} else /* ptr->na == 2 */
{ na = IntList->inext;
ptr->dist = na->dist - cosa*dist;
}
ptr->angle = sina*dist;
ptr->dihed = 0.0;
while( ptr=ptr->inext )
{ dist = ptr->dist;
angle = Deg2Rad*ptr->angle;
dihed = Deg2Rad*ptr->dihed;
/* Optimise this access?? */
na = GetInternalCoord(ptr->na);
nb = GetInternalCoord(ptr->nb);
nc = GetInternalCoord(ptr->nc);
xb = nb->dist - na->dist;
yb = nb->angle - na->angle;
zb = nb->dihed - na->dihed;
rbc = xb*xb + yb*yb + zb*zb;
if( rbc < 0.0001 )
return( False );
rbc= 1.0/sqrt(rbc);
cosa = cos(angle);
sina = sin(angle);
if( fabs(cosa) >= 0.999999 )
{ /* Colinear */
temp = dist*rbc*cosa;
ptr->dist = na->dist + temp*xb;
ptr->angle = na->angle + temp*yb;
ptr->dihed = na->dihed + temp*zb;
} else
{ xa = nc->dist - na->dist;
ya = nc->angle - na->angle;
za = nc->dihed - na->dihed;
sind = -sin(dihed);
cosd = cos(dihed);
xd = dist*cosa;
yd = dist*sina*cosd;
zd = dist*sina*sind;
xyb = sqrt(xb*xb + yb*yb);
if( xyb < 0.1 )
{ /* Rotate about y-axis! */
temp = za; za = -xa; xa = temp;
temp = zb; zb = -xb; xb = temp;
xyb = sqrt(xb*xb + yb*yb);
flag = True;
} else flag = False;
costh = xb/xyb;
sinth = yb/xyb;
xpa = costh*xa + sinth*ya;
ypa = costh*ya - sinth*xa;
sinph = zb*rbc;
cosph = sqrt(1.0 - sinph*sinph);
zqa = cosph*za - sinph*xpa;
yza = sqrt(ypa*ypa + zqa*zqa);
if( yza > 1.0E-10 )
{ coskh = ypa/yza;
sinkh = zqa/yza;
ypd = coskh*yd - sinkh*zd;
zpd = coskh*zd + sinkh*yd;
} else
{ /* coskh = 1.0; */
/* sinkh = 0.0; */
ypd = yd;
zpd = zd;
}
xpd = cosph*xd - sinph*zpd;
zqd = cosph*zpd + sinph*xd;
xqd = costh*xpd - sinth*ypd;
yqd = costh*ypd + sinth*xpd;
if( flag )
{ /* Rotate about y-axis! */
ptr->dist = na->dist - zqd;
ptr->angle = na->angle + yqd;
ptr->dihed = na->dihed + xqd;
} else
{ ptr->dist = na->dist + xqd;
ptr->angle = na->angle + yqd;
ptr->dihed = na->dihed + zqd;
}
}
}
return( True );
}
/*=========================*/
/* Bond Handling Functions */
/*=========================*/
#ifdef FUNCPROTO
Bond __far *ProcessBond( Atom __far*, Atom __far*, int );
static void CreateHydrogenBond( Atom __far*, Atom __far*,
Atom __far*, Atom __far*, int, int );
#endif
Bond __far *ProcessBond( src, dst, flag )
Atom __far *src, __far *dst;
int flag;
{
register Bond __far *ptr;
register int i;
if( flag & (DoubBondFlag|TripBondFlag) )
DrawDoubleBonds = True;
if( !(ptr = FreeBond) )
{ MemSize += BondPool*sizeof(Bond);
ptr = (Bond __far *)_fmalloc( BondPool*sizeof(Bond) );
if( !ptr ) FatalDataError("Memory allocation failed");
RegisterAlloc( ptr );
for( i=1; i<BondPool; i++ )
{ ptr->bnext = FreeBond;
FreeBond = ptr++;
}
} else FreeBond = ptr->bnext;
ptr->flag = flag | SelectFlag;
ptr->srcatom = src;
ptr->dstatom = dst;
ptr->radius = 0;
ptr->col = 0;
return( ptr );
}
static void CreateHydrogenBond( srcCA, dstCA, src, dst, energy, offset )
Atom __far *srcCA, __far *dstCA;
Atom __far *src, __far *dst;
int energy, offset;
{
register HBond __far *ptr;
register int i,flag;
if( !(ptr = FreeHBond) )
{ MemSize += HBondPool*sizeof(HBond);
ptr = (HBond __far *)_fmalloc( HBondPool*sizeof(HBond) );
if( !ptr ) FatalDataError("Memory allocation failed");
RegisterAlloc( ptr );
for( i=1; i<HBondPool; i++ )
{ ptr->hnext = FreeHBond;
FreeHBond = ptr++;
}
} else FreeHBond = ptr->hnext;
if( (offset>=-128) && (offset<127) )
{ ptr->offset = (Char)offset;
} else ptr->offset = 0;
flag = ZoneBoth? src->flag&dst->flag : src->flag|dst->flag;
ptr->flag = flag & SelectFlag;
ptr->src = src;
ptr->dst = dst;
ptr->srcCA = srcCA;
ptr->dstCA = dstCA;
ptr->energy = energy;
ptr->col = 0;
*CurHBond = ptr;
ptr->hnext = (void __far*)0;
CurHBond = &ptr->hnext;
InfoHBondCount++;
}
void CreateBond( src, dst, flag )
int src, dst, flag;
{
register Chain __far *chain;
register Group __far *group;
register Atom __far *aptr;
register Atom __far *sptr;
register Atom __far *dptr;
register Bond __far *bptr;
register int done;
if( src == dst )
return;
sptr = (void __far*)0;
dptr = (void __far*)0;
done = False;
for( chain=Database->clist; chain && !done; chain=chain->cnext )
for( group=chain->glist; group && !done; group=group->gnext )
for( aptr=group->alist; aptr; aptr=aptr->anext )
{ if( aptr->serno == src )
{ sptr = aptr;
if( dptr )
{ done = True;
break;
}
} else if( aptr->serno == dst )
{ dptr = aptr;
if( sptr )
{ done = True;
break;
}
}
}
/* Both found! */
if( done )
{ if( flag )
{ /* Reset Non-bonded flags! */
sptr->flag &= ~NonBondFlag;
dptr->flag &= ~NonBondFlag;
bptr = ProcessBond( sptr, dptr, flag );
bptr->bnext = CurMolecule->blist;
CurMolecule->blist = bptr;
InfoBondCount++;
} else /* Hydrogen Bond! */
{ if( InfoHBondCount<0 )
{ CurHBond = &CurMolecule->hlist;
InfoHBondCount = 0;
}
CreateHydrogenBond( NULL, NULL, sptr, dptr, 0, 0 );
}
}
}
void CreateBondOrder( src, dst )
int src, dst;
{
register Bond __far *bptr;
register int bs,bd;
ForEachBond
{ bs = bptr->srcatom->serno;
bd = bptr->dstatom->serno;
if( ((bs==src)&&(bd==dst)) || ((bs==dst)&&(bd==src)) )
{ DrawDoubleBonds = True;
if( bptr->flag & NormBondFlag )
{ /* Convert Single to Double */
bptr->flag &= ~(NormBondFlag);
bptr->flag |= DoubBondFlag;
} else if( bptr->flag & DoubBondFlag )
{ /* Convert Double to Triple */
bptr->flag &= ~(DoubBondFlag);
bptr->flag |= TripBondFlag;
}
return;
}
}
CreateBond( src, dst, NormBondFlag );
}
static void TestBonded( sptr, dptr, flag )
Atom __far *sptr, __far *dptr;
int flag;
{
register Bond __far *bptr;
register Long dx, dy, dz;
register Long max, dist;
if( flag )
{ /* Sum of covalent radii with 0.56A tolerance */
dist = Element[sptr->elemno].covalrad +
Element[dptr->elemno].covalrad + 140;
max = dist*dist;
} else
{ /* Fast Bio-Macromolecule Bonding Calculation */
if( (sptr->flag|dptr->flag) & HydrogenFlag )
{ max = MaxHBondDist;
} else max = MaxBondDist;
}
dx = sptr->xorg-dptr->xorg; if( (dist=dx*dx)>max ) return;
dy = sptr->yorg-dptr->yorg; if( (dist+=dy*dy)>max ) return;
dz = sptr->zorg-dptr->zorg; if( (dist+=dz*dz)>max ) return;
if( dist > MinBondDist )
{ /* Reset Non-bonded flags! */
sptr->flag &= ~NonBondFlag;
dptr->flag &= ~NonBondFlag;
bptr = ProcessBond(sptr,dptr,NormBondFlag);
bptr->bnext = CurMolecule->blist;
CurMolecule->blist = bptr;
InfoBondCount++;
}
}
static void ReclaimHBonds( ptr )
HBond __far *ptr;
{
register HBond __far *temp;
if( (temp = ptr) )
{ while( temp->hnext )
temp=temp->hnext;
temp->hnext = FreeHBond;
FreeHBond = ptr;
}
}
static void ReclaimBonds( ptr )
Bond __far *ptr;
{
register Bond __far *temp;
if( (temp = ptr) )
{ while( temp->bnext )
temp=temp->bnext;
temp->bnext = FreeBond;
FreeBond = ptr;
}
}
void CreateMoleculeBonds( info, flag )
int info, flag;
{
register int i, x, y, z;
register Long tx, ty, tz;
register Long mx, my, mz;
register Long dx, dy, dz;
register int lx, ly, lz, hx, hy, hz;
register Atom __far *aptr, __far *dptr;
register Chain __far *chain;
register Group __far *group;
char buffer[40];
if( !Database )
return;
dx = (MaxX-MinX)+1;
dy = (MaxY-MinY)+1;
dz = (MaxZ-MinZ)+1;
InfoBondCount = 0;
ReclaimBonds( CurMolecule->blist );
CurMolecule->blist = (void __far*)0;
ResetVoxelData();
for( chain=Database->clist; chain; chain=chain->cnext )
{ ResetVoxelData();
for( group=chain->glist; group; group=group->gnext )
for( aptr=group->alist; aptr; aptr=aptr->anext )
{ /* Initially non-bonded! */
aptr->flag |= NonBondFlag;
mx = aptr->xorg-MinX;
my = aptr->yorg-MinY;
mz = aptr->zorg-MinZ;
tx = mx-AbsMaxBondDist;
ty = my-AbsMaxBondDist;
tz = mz-AbsMaxBondDist;
lx = (tx>0)? (int)((VOXORDER*tx)/dx) : 0;
ly = (ty>0)? (int)((VOXORDER*ty)/dy) : 0;
lz = (tz>0)? (int)((VOXORDER*tz)/dz) : 0;
tx = mx+AbsMaxBondDist;
ty = my+AbsMaxBondDist;
tz = mz+AbsMaxBondDist;
hx = (tx<dx)? (int)((VOXORDER*tx)/dx) : VOXORDER-1;
hy = (ty<dy)? (int)((VOXORDER*ty)/dy) : VOXORDER-1;
hz = (tz<dz)? (int)((VOXORDER*tz)/dz) : VOXORDER-1;
for( x=lx; x<=hx; x++ )
{ i = VOXORDER2*x + VOXORDER*ly;
for( y=ly; y<=hy; y++ )
{ for( z=lz; z<=hz; z++ )
if( (dptr = (Atom __far*)HashTable[i+z]) )
do { TestBonded(aptr,dptr,flag);
} while( (dptr = dptr->next) );
i += VOXORDER;
}
}
x = (int)((VOXORDER*mx)/dx);
y = (int)((VOXORDER*my)/dy);
z = (int)((VOXORDER*mz)/dz);
i = VOXORDER2*x + VOXORDER*y + z;
aptr->next = (Atom __far*)HashTable[i];
HashTable[i] = (void __far*)aptr;
}
VoxelsClean = False;
}
if( info )
{ if( CommandActive )
WriteChar('\n');
CommandActive=False;
sprintf(buffer,"Number of Bonds ..... %ld\n\n",(long)InfoBondCount);
WriteString(buffer);
}
}
/*===============================*/
/* Disulphide bridging functions */
/*===============================*/
#ifdef FUNCPROTO
static Atom __far *FindCysSulphur( Group __far* );
#endif
static Atom __far *FindCysSulphur( group )
Group __far *group;
{
register Atom __far *ptr;
register char *elem;
for( ptr=group->alist; ptr; ptr=ptr->anext )
{ elem = ElemDesc[ptr->refno];
if( (elem[1]=='S') && (elem[0]==' ') )
return( ptr );
}
return( (Atom __far*)0 );
}
static void TestDisulphideBridge( group1, group2, cys1 )
Group __far *group1, __far *group2;
Atom __far *cys1;
{
register HBond __far *ptr;
register Atom __far *cys2;
register int dx, dy, dz;
register Long max,dist;
if( !(cys2=FindCysSulphur(group2)) )
return;
max = (Long)750*750;
dx = (int)(cys1->xorg-cys2->xorg); if( (dist=(Long)dx*dx)>max ) return;
dy = (int)(cys1->yorg-cys2->yorg); if( (dist+=(Long)dy*dy)>max ) return;
dz = (int)(cys1->zorg-cys2->zorg); if( (dist+=(Long)dz*dz)>max ) return;
if( !(ptr = FreeHBond) )
{ MemSize += sizeof(HBond);
ptr = (HBond __far*)_fmalloc(sizeof(HBond));
if( !ptr ) FatalDataError("Memory allocation failed");
RegisterAlloc( ptr );
} else FreeHBond = ptr->hnext;
ptr->dst = cys1;
if( !(ptr->dstCA=FindGroupAtom(group1,1)) )
ptr->dstCA = cys1;
ptr->src = cys2;
if( !(ptr->srcCA=FindGroupAtom(group2,1)) )
ptr->srcCA = cys2;
ptr->offset = 0;
ptr->energy = 0;
ptr->flag = 0;
ptr->col = 0;
ptr->hnext = CurMolecule->slist;
CurMolecule->slist = ptr;
group1->flag |= CystineFlag;
group2->flag |= CystineFlag;
InfoSSBondCount++;
}
void FindDisulphideBridges()
{
register Chain __far *chn1;
register Chain __far *chn2;
register Group __far *group1;
register Group __far *group2;
register Atom __far *cys;
char buffer[40];
if( !Database ) return;
ReclaimHBonds( CurMolecule->slist );
InfoSSBondCount = 0;
for(chn1=Database->clist;chn1;chn1=chn1->cnext)
for(group1=chn1->glist;group1;group1=group1->gnext)
if( IsCysteine(group1->refno) && (cys=FindCysSulphur(group1)) )
{ for(group2=group1->gnext;group2;group2=group2->gnext)
if( IsCysteine(group2->refno) )
TestDisulphideBridge(group1,group2,cys);
for(chn2=chn1->cnext;chn2;chn2=chn2->cnext)
for(group2=chn2->glist;group2;group2=group2->gnext)
if( IsCysteine(group2->refno) )
TestDisulphideBridge(group1,group2,cys);
}
if( FileDepth == -1 )
{ if( CommandActive )
WriteChar('\n');
CommandActive=False;
sprintf(buffer,"Number of Bridges ... %d\n\n",InfoSSBondCount);
WriteString(buffer);
}
}
/*=========================================*/
/* Kabsch & Sander Structure Determination */
/*=========================================*/
#ifdef FUNCPROTO
static int CalculateBondEnergy( Group __far* );
static void CalcProteinHBonds( Chain __far* );
static void CalcNucleicHBonds( Chain __far* );
static int IsHBonded( Atom __far*, Atom __far*, HBond __far* );
static void TestLadder( Chain __far* );
#endif
/* Coupling constant for Electrostatic Energy */
/* QConst = -332 * 0.42 * 0.2 * 1000.0 [*250.0] */
#define QConst (-6972000.0)
#define MaxHDist ((Long)2250*2250)
#define MinHDist ((Long)125*125)
/* Protein Donor Atom Coordinates */
static int hxorg,hyorg,hzorg;
static int nxorg,nyorg,nzorg;
static Atom __far *best1CA;
static Atom __far *best2CA;
static Atom __far *best1;
static Atom __far *best2;
static Atom __far *optr;
static int res1,res2;
static int off1,off2;
static int CalculateBondEnergy( group )
Group __far *group;
{
register double dho,dhc;
register double dnc,dno;
register Atom __far *cptr;
register Long dx,dy,dz;
register Long dist;
register int result;
if( !(cptr=FindGroupAtom(group,2)) ) return(0);
if( !(optr=FindGroupAtom(group,3)) ) return(0);
dx = hxorg-optr->xorg;
dy = hyorg-optr->yorg;
dz = hzorg-optr->zorg;
dist = dx*dx+dy*dy+dz*dz;
if( dist < MinHDist )
return( -9900 );
dho = sqrt((double)dist);
dx = hxorg-cptr->xorg;
dy = hyorg-cptr->yorg;
dz = hzorg-cptr->zorg;
dist = dx*dx+dy*dy+dz*dz;
if( dist < MinHDist )
return( -9900 );
dhc = sqrt((double)dist);
dx = nxorg-cptr->xorg;
dy = nyorg-cptr->yorg;
dz = nzorg-cptr->zorg;
dist = dx*dx+dy*dy+dz*dz;
if( dist < MinHDist )
return( -9900 );
dnc = sqrt((double)dist);
dx = nxorg-optr->xorg;
dy = nyorg-optr->yorg;
dz = nzorg-optr->zorg;
dist = dx*dx+dy*dy+dz*dz;
if( dist < MinHDist )
return( -9900 );
dno = sqrt((double)dist);
result = (int)(QConst/dho - QConst/dhc + QConst/dnc - QConst/dno);
if( result<-9900 )
{ return( -9900 );
} else if( result>-500 )
return( 0 );
return( result );
}
static void CalcProteinHBonds( chn1 )
Chain __far *chn1;
{
register int energy, offset;
register Chain __far *chn2;
register Group __far *group1;
register Group __far *group2;
register Atom __far *ca1;
register Atom __far *ca2;
register Atom __far *pc1;
register Atom __far *po1;
register Atom __far *n1;
register int pos1,pos2;
register int dx,dy,dz;
register double dco;
register Long dist;
pos1 = 0;
pc1 = po1 = (void __far*)0;
for(group1=chn1->glist;group1;group1=group1->gnext)
{ pos1++;
if( pc1 && po1 )
{ dx = (int)(pc1->xorg - po1->xorg);
dy = (int)(pc1->yorg - po1->yorg);
dz = (int)(pc1->zorg - po1->zorg);
} else
{ pc1 = FindGroupAtom(group1,2);
po1 = FindGroupAtom(group1,3);
continue;
}
pc1 = FindGroupAtom(group1,2);
po1 = FindGroupAtom(group1,3);
if( !IsAmino(group1->refno) || IsProline(group1->refno) )
continue;
if( !(ca1=FindGroupAtom(group1,1)) ) continue;
if( !(n1=FindGroupAtom(group1,0)) ) continue;
dist = (Long)dx*dx + (Long)dy*dy + (Long)dz*dz;
dco = sqrt( (double)dist )/250.0;
nxorg = (int)n1->xorg; hxorg = nxorg + (int)(dx/dco);
nyorg = (int)n1->yorg; hyorg = nyorg + (int)(dy/dco);
nzorg = (int)n1->zorg; hzorg = nzorg + (int)(dz/dco);
res1 = res2 = 0;
/* Only Hydrogen Bond within a single chain! */
/* for(chn2=Database->clist;chn2;chn2=chn2->cnext) */
chn2 = chn1;
{ /* Only consider non-empty peptide chains! */
/* if( !chn2->glist || !IsProtein(chn2->glist->refno) ) */
/* continue; */
pos2 = 0;
for(group2=chn2->glist;group2;group2=group2->gnext)
{ pos2++;
if( (group2==group1) || (group2->gnext==group1) )
continue;
if( !IsAmino(group2->refno) )
continue;
if( !(ca2=FindGroupAtom(group2,1)) )
continue;
dx = (int)(ca1->xorg-ca2->xorg);
if( (dist=(Long)dx*dx) > MaxHDist )
continue;
dy = (int)(ca1->yorg-ca2->yorg);
if( (dist+=(Long)dy*dy) > MaxHDist )
continue;
dz = (int)(ca1->zorg-ca2->zorg);
if( (dist+=(Long)dz*dz) > MaxHDist )
continue;
if( (energy = CalculateBondEnergy(group2)) )
{ if( chn1 == chn2 )
{ offset = pos1 - pos2;
} else offset = 0;
if( energy<res1 )
{ best2CA = best1CA; best1CA = ca2;
best2 = best1; best1 = optr;
res2 = res1; res1 = energy;
off2 = off1; off1 = offset;
} else if( energy<res2 )
{ best2CA = ca2;
best2 = optr;
res2 = energy;
off2 = offset;
}
}
} /* group2 */
} /* chn2 */
if( res1 )
{ if( res2 )
CreateHydrogenBond(ca1,best2CA,n1,best2,res2,off2);
CreateHydrogenBond(ca1,best1CA,n1,best1,res1,off1);
}
}
}
static void CalcNucleicHBonds( chn1 )
Chain __far *chn1;
{
register Chain __far *chn2;
register Group __far *group1;
register Group __far *group2;
register Group __far *best;
register Atom __far *ca1;
register Atom __far *ca2;
register Atom __far *n1;
register Long max,dist;
register int dx,dy,dz;
register int refno;
for(group1=chn1->glist;group1;group1=group1->gnext)
{ if( !IsPurine(group1->refno) ) continue;
/* Find N1 of Purine Group */
if( !(n1=FindGroupAtom(group1,21)) )
continue;
/* Maximum N1-N3 distance 5A */
refno = NucleicCompl(group1->refno);
max = (Long)1250*1250;
best = (void __far*)0;
for(chn2=Database->clist;chn2;chn2=chn2->cnext)
{ /* Only consider non-empty nucleic acid chains! */
if( (chn1==chn2) || !chn2->glist ||
!IsDNA(chn2->glist->refno) )
continue;
for(group2=chn2->glist;group2;group2=group2->gnext)
if( group2->refno == (Byte)refno )
{ /* Find N3 of Pyramidine Group */
if( !(ca1=FindGroupAtom(group2,23)) )
continue;
dx = (int)(ca1->xorg - n1->xorg);
if( (dist=(Long)dx*dx) >= max )
continue;
dy = (int)(ca1->yorg - n1->yorg);
if( (dist+=(Long)dy*dy) >= max )
continue;
dz = (int)(ca1->zorg - n1->zorg);
if( (dist+=(Long)dz*dz) >= max )
continue;
best1 = ca1;
best = group2;
max = dist;
}
}
if( best )
{ /* Find the sugar phosphorous atoms */
ca1 = FindGroupAtom( group1, 7 );
ca2 = FindGroupAtom( best, 7 );
CreateHydrogenBond( ca1, ca2, n1, best1, 0, 0 );
if( IsGuanine(group1->refno) )
{ /* Guanine-Cytosine */
if( (ca1=FindGroupAtom(group1,22)) && /* G.N2 */
(ca2=FindGroupAtom(best,26)) ) /* C.O2 */
CreateHydrogenBond( (void __far*)0, (void __far*)0,
ca1, ca2, 0, 0 );
if( (ca1=FindGroupAtom(group1,28)) && /* G.O6 */
(ca2=FindGroupAtom(best,24)) ) /* C.N4 */
CreateHydrogenBond( (void __far*)0, (void __far*)0,
ca1, ca2, 0, 0 );
} else /* Adenine-Thymine */
if( (ca1=FindGroupAtom(group1,25)) && /* A.N6 */
(ca2=FindGroupAtom(best,27)) ) /* T.O4 */
CreateHydrogenBond( (void __far*)0, (void __far*)0,
ca1, ca2, 0, 0 );
}
}
}
void CalcHydrogenBonds()
{
register Chain __far *chn1;
char buffer[40];
if( !Database ) return;
ReclaimHBonds( CurMolecule->hlist );
CurMolecule->hlist = (void __far*)0;
CurHBond = &CurMolecule->hlist;
InfoHBondCount = 0;
if( MainAtomCount > 10000 )
{ if( CommandActive )
WriteChar('\n');
WriteString("Please wait... ");
CommandActive=True;
}
for(chn1=Database->clist; chn1; chn1=chn1->cnext)
if( chn1->glist )
{ if( IsProtein(chn1->glist->refno) )
{ CalcProteinHBonds(chn1);
} else if( IsDNA(chn1->glist->refno) )
CalcNucleicHBonds(chn1);
}
if( FileDepth == -1 )
{ if( CommandActive )
WriteChar('\n');
CommandActive=False;
sprintf(buffer,"Number of H-Bonds ... %d\n",InfoHBondCount);
WriteString(buffer);
}
}
static int IsHBonded( src, dst, ptr )
Atom __far *src, __far *dst;
HBond __far *ptr;
{
while( ptr && (ptr->srcCA==src) )
if( ptr->dstCA == dst )
{ return( True );
} else ptr=ptr->hnext;
return( False );
}
static void FindAlphaHelix( pitch, flag )
int pitch, flag;
{
register HBond __far *hbond;
register Chain __far *chain;
register Group __far *group;
register Group __far *first;
register Group __far *ptr;
register Atom __far *srcCA;
register Atom __far *dstCA;
register int res,dist,prev;
/* Protein chains only! */
hbond = Database->hlist;
for( chain=Database->clist; chain; chain=chain->cnext )
if( (first=chain->glist) && IsProtein(first->refno) )
{ prev = False; dist = 0;
for( group=chain->glist; hbond && group; group=group->gnext )
{ if( IsAmino(group->refno) && (srcCA=FindGroupAtom(group,1)) )
{ if( dist==pitch )
{ res = False;
dstCA=FindGroupAtom(first,1);
while( hbond && hbond->srcCA == srcCA )
{ if( hbond->dstCA==dstCA ) res=True;
hbond = hbond->hnext;
}
if( res )
{ if( prev )
{ if( !(first->struc&HelixFlag) )
InfoHelixCount++;
ptr = first;
do {
ptr->struc |= flag;
ptr = ptr->gnext;
} while( ptr != group );
} else prev = True;
} else prev = False;
} else while( hbond && hbond->srcCA==srcCA )
hbond = hbond->hnext;
} else prev = False;
if( group->struc&HelixFlag )
{ first = group; prev = False; dist = 1;
} else if( dist==pitch )
{ first = first->gnext;
} else dist++;
}
} else if( first && IsNucleo(first->refno) )
while( hbond && !IsAminoBackbone(hbond->src->refno) )
hbond = hbond->hnext;
}
static Atom __far *cprevi, __far *ccurri, __far *cnexti;
static HBond __far *hcurri, __far *hnexti;
static Group __far *curri, __far *nexti;
static void TestLadder( chain )
Chain __far *chain;
{
register Atom __far *cprevj, __far *ccurrj, __far *cnextj;
register HBond __far *hcurrj, __far *hnextj;
register Group __far *currj, __far *nextj;
register int count, result, found;
/* Already part of atleast one ladder */
found = curri->flag & SheetFlag;
nextj = nexti->gnext;
hnextj = hnexti;
while( hnextj && hnextj->srcCA==cnexti )
hnextj = hnextj->hnext;
while( True )
{ if( nextj )
if( IsProtein(chain->glist->refno) )
{ count = 1;
do {
cnextj = FindGroupAtom(nextj,1);
if( count == 3 )
{ if( IsHBonded(cnexti,ccurrj,hnexti) &&
IsHBonded(ccurrj,cprevi,hcurrj) )
{ result = ParaLadder;
} else if( IsHBonded(cnextj,ccurri,hnextj) &&
IsHBonded(ccurri,cprevj,hcurri) )
{ result = ParaLadder;
} else if( IsHBonded(cnexti,cprevj,hnexti) &&
IsHBonded(cnextj,cprevi,hnextj) )
{ result = AntiLadder;
} else if( IsHBonded(ccurrj,ccurri,hcurrj) &&
IsHBonded(ccurri,ccurrj,hcurri) )
{ result = AntiLadder;
} else result = NoLadder;
if( result )
{ curri->struc |= SheetFlag;
currj->struc |= SheetFlag;
if( found ) return;
found = True;
}
} else count++;
cprevj = ccurrj; ccurrj = cnextj;
currj = nextj; hcurrj = hnextj;
while( hnextj && hnextj->srcCA==cnextj )
hnextj = hnextj->hnext;
} while( (nextj = nextj->gnext) );
} else if( IsNucleo(chain->glist->refno) )
while( hnextj && !IsAminoBackbone(hnextj->src->refno) )
hnextj = hnextj->hnext;
if( (chain = chain->cnext) )
{ nextj = chain->glist;
} else return;
}
}
static void FindBetaSheets()
{
register Chain __far *chain;
register int ladder;
register int count;
hnexti = Database->hlist;
for( chain=Database->clist; chain; chain=chain->cnext )
if( (nexti = chain->glist) )
if( IsProtein(nexti->refno) )
{ count = 1;
ladder = False;
do {
cnexti = FindGroupAtom(nexti,1);
if( count == 3 )
{ TestLadder( chain );
if( curri->struc & SheetFlag )
{ if( !ladder )
{ InfoLadderCount++;
ladder = True;
}
} else ladder = False;
} else count++;
cprevi = ccurri; ccurri = cnexti;
curri = nexti; hcurri = hnexti;
while( hnexti && hnexti->srcCA==cnexti )
hnexti = hnexti->hnext;
} while( (nexti = nexti->gnext) );
} else if( IsNucleo(nexti->refno) )
while( hnexti && !IsAminoBackbone(hnexti->src->refno) )
hnexti = hnexti->hnext;
}
static void FindTurnStructure()
{
static Atom __far *aptr[5];
register Chain __far *chain;
register Group __far *group;
register Group __far *prev;
register Atom __far *ptr;
register Long ux,uy,uz,mu;
register Long vx,vy,vz,mv;
register int i,found,len;
register Real CosKappa;
for( chain=Database->clist; chain; chain=chain->cnext )
if( chain->glist && IsProtein(chain->glist->refno) )
{ len = 0; found = False;
for( group=chain->glist; group; group=group->gnext )
{ ptr = FindGroupAtom(group,1);
if( ptr && (ptr->flag&BreakFlag) )
{ found = False;
len = 0;
} else if( len==5 )
{ for( i=0; i<4; i++ )
aptr[i] = aptr[i+1];
len = 4;
} else if( len==2 )
prev = group;
aptr[len++] = ptr;
if( len==5 )
{ if( !(prev->struc&(HelixFlag|SheetFlag)) &&
aptr[0] && aptr[2] && aptr[4] )
{ ux = aptr[2]->xorg - aptr[0]->xorg;
uy = aptr[2]->yorg - aptr[0]->yorg;
uz = aptr[2]->zorg - aptr[0]->zorg;
vx = aptr[4]->xorg - aptr[2]->xorg;
vy = aptr[4]->yorg - aptr[2]->yorg;
vz = aptr[4]->zorg - aptr[2]->zorg;
mu = ux*ux + uy*uy + uz*uz;
mv = vx*vx + vz*vz + vy*vy;
if( mu && mv )
{ CosKappa = (Real)(ux*vx + uy*vy + uz*vz);
CosKappa /= sqrt( (Real)mu*mv );
if( CosKappa<Cos70Deg )
{ if( !found )
InfoTurnCount++;
prev->struc |= TurnFlag;
}
}
}
found = prev->struc&TurnFlag;
prev = prev->gnext;
} /* len==5 */
}
}
}
static void FindBetaTurns()
{
static Atom __far *aptr[4];
register Chain __far *chain;
register Group __far *group;
register Group __far *prev;
register Group __far *next;
register Atom __far *ptr;
register Long dx,dy,dz;
register int found,len;
register int flag;
for( chain=Database->clist; chain; chain=chain->cnext )
if( chain->glist && IsProtein(chain->glist->refno) )
{ prev = chain->glist;
len = 0; found = False;
for( next=chain->glist; next; next=next->gnext )
{ ptr = FindGroupAtom(next,1);
if( ptr && (ptr->flag&BreakFlag) )
{ found = False;
prev = next;
len = 0;
} else if( len==4 )
{ aptr[0] = aptr[1];
aptr[1] = aptr[2];
aptr[2] = aptr[3];
aptr[3] = ptr;
} else aptr[len++] = ptr;
if( len==4 )
{ flag = False;
if( aptr[0] && aptr[3] )
{ dx = aptr[3]->xorg - aptr[0]->xorg;
dy = aptr[3]->yorg - aptr[0]->yorg;
dz = aptr[3]->zorg - aptr[0]->zorg;
if( dx*dx + dy*dy + dz*dz < (Long)1750*1750 )
{ group = prev;
while( group!=next->gnext )
{ if( !(group->struc&(HelixFlag|SheetFlag)) )
{ group->struc |= TurnFlag;
flag = True;
}
group = group->gnext;
}
if( !found && flag )
InfoTurnCount++;
}
}
prev = prev->gnext;
found = flag;
} /* len==4 */
}
}
}
void DetermineStructure( flag )
int flag;
{
register Chain __far *chain;
register Group __far *group;
char buffer[40];
if( !Database )
return;
if( InfoHBondCount<0 )
CalcHydrogenBonds();
if( InfoHelixCount>=0 )
for( chain=Database->clist; chain; chain=chain->cnext )
for( group=chain->glist; group; group=group->gnext )
group->struc = 0;
InfoStrucSrc = SourceCalc;
InfoLadderCount = 0;
InfoHelixCount = 0;
InfoTurnCount = 0;
if( InfoHBondCount )
{ FindAlphaHelix(4,Helix4Flag);
FindBetaSheets();
FindAlphaHelix(3,Helix3Flag);
FindAlphaHelix(5,Helix5Flag);
if( !flag )
{ FindTurnStructure();
} else FindBetaTurns();
}
if( FileDepth == -1 )
{ if( CommandActive )
WriteChar('\n');
CommandActive=False;
sprintf(buffer,"Number of Helices ... %d\n",InfoHelixCount);
WriteString(buffer);
sprintf(buffer,"Number of Strands ... %d\n",InfoLadderCount);
WriteString(buffer);
sprintf(buffer,"Number of Turns ..... %d\n",InfoTurnCount);
WriteString(buffer);
}
}
void RenumberMolecule( start )
int start;
{
register Chain __far *chain;
register Group __far *group;
register int hinit, minit;
register int resno;
if( !Database )
return;
hinit = minit = False;
for( chain=Database->clist; chain; chain=chain->cnext )
{ resno = start;
for( group=chain->glist; group; group=group->gnext )
{ if( group->alist->flag & HeteroFlag )
{ if( hinit )
{ if( resno > MaxHetaRes )
{ MaxHetaRes = resno;
} else if( resno < MinHetaRes )
MinHetaRes = resno;
} else MinHetaRes = MaxHetaRes = resno;
hinit = True;
} else
{ if( minit )
{ if( resno > MaxMainRes )
{ MaxMainRes = resno;
} else if( resno < MinMainRes )
MinMainRes = resno;
} else MinMainRes = MaxMainRes = resno;
minit = True;
}
group->serno = resno++;
}
}
}
/*===============================*/
/* Molecule Database Maintenance */
/*===============================*/
static void ReclaimAtoms( ptr )
Atom __far *ptr;
{
register Atom __far *temp;
if( (temp = ptr) )
{ while( temp->anext )
temp=temp->anext;
temp->anext = FreeAtom;
FreeAtom = ptr;
}
}
static void ResetDatabase()
{
Database = CurMolecule = (void __far*)0;
MainGroupCount = HetaGroupCount = 0;
InfoChainCount = HetaAtomCount = 0;
MainAtomCount = InfoBondCount = 0;
SelectCount = 0;
InfoStrucSrc = SourceNone;
InfoSSBondCount = InfoHBondCount = -1;
InfoHelixCount = InfoLadderCount = -1;
InfoTurnCount = -1;
CurGroup = (void __far*)0;
CurChain = (void __far*)0;
CurAtom = (void __far*)0;
MinX = MinY = MinZ = 0;
MaxX = MaxY = MaxZ = 0;
MinMainTemp = MaxMainTemp = 0;
MinHetaTemp = MaxHetaTemp = 0;
MinMainRes = MaxMainRes = 0;
MinHetaRes = MaxHetaRes = 0;
*InfoMoleculeName = 0;
*InfoClassification = 0;
*InfoIdentCode = 0;
*InfoSpaceGroup = 0;
*InfoFileName = 0;
VoxelsClean = False;
HMinMaxFlag = False;
MMinMaxFlag = False;
HasHydrogen = False;
ElemNo = MINELEM;
ResNo = MINRES;
MaskCount = 0;
NMRModel = 0;
}
void DestroyDatabase()
{
register void __far *temp;
register Group __far *gptr;
if( Database )
{ ReclaimHBonds( Database->slist );
ReclaimHBonds( Database->hlist );
ReclaimBonds( Database->blist );
while( Database->clist )
{ ReclaimBonds(Database->clist->blist);
if( (gptr = Database->clist->glist) )
{ ReclaimAtoms(gptr->alist);
while( gptr->gnext )
{ gptr = gptr->gnext;
ReclaimAtoms(gptr->alist);
}
gptr->gnext = FreeGroup;
FreeGroup = Database->clist->glist;
}
temp = Database->clist->cnext;
Database->clist->cnext = FreeChain;
FreeChain = Database->clist;
Database->clist = temp;
}
FreeMolecule = Database;
Database = (void __far*)0;
}
ResetDatabase();
}
void PurgeDatabase()
{
#ifdef APPLEMAC
register AllocRef *ptr;
register AllocRef *tmp;
register int i;
/* Avoid Memory Leaks */
for( ptr=AllocList; ptr; ptr=tmp )
{ for( i=0; i<ptr->count; i++ )
_ffree( ptr->data[i] );
tmp = ptr->next;
_ffree( ptr );
}
#endif
}
void InitialiseDatabase()
{
FreeMolecule = (void __far*)0;
FreeHBond = (void __far*)0;
FreeChain = (void __far*)0;
FreeGroup = (void __far*)0;
FreeAtom = (void __far*)0;
FreeBond = (void __far*)0;
#ifdef APPLEMAC
AllocList = (void*)0;
#endif
ResetDatabase();
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.