in contrib/champ/champ.c [489:1018]
void ChampCountRings(CChamp *I,int index)
{
ListPat *pat;
ListBond *bd,*bd0,*bd1,*bd2,*bd3;
ListAtom *at,*at0,*at1,*at2,*at3,*at4;
int a,a0,a1,ai,bi,ai0,ai1,ai2,a2,bd_a0,bd_a1,ai3,ai4;
int i2,i0,i1,i3;
int bi0,bi1,bi2,bi3;
int ni1;
int n_atom = 0;
int ring_size;
pat = I->Pat + index;
ai = pat->atom;
while(ai) { /* count number of atoms */
n_atom++;
at = I->Atom + ai;
ai = at->link;
}
if(n_atom) {
ListAtom **atom = NULL;
int *atom_idx;
int *bonds = NULL;
int *neighbors = NULL;
int lin_mask_size;
int **lin_mask_vla = NULL;
int *lin_mask;
int lin_index;
int n_expand;
/* create a packed array of the atom pointers -- these will be our
atom indices for this routine */
atom = mac_malloc_array(ListAtom*,n_atom);
atom_idx = mac_malloc_array(int,n_atom);
a = 0;
ai = pat->atom;
while(ai) { /* load array and provide back-reference in mark_targ */
at = I->Atom + ai;
at->cycle = 0; /* initialize */
at->class = 0;
atom[a] = at;
atom_idx[a] = ai;
at->mark_targ = a; /* used for back-link */
at->mark_read = 0; /* used to mark exclusion */
at->mark_tmpl = 0; /* used to mark as dirty */
at->first_targ = 0; /* used for lineage list */
at->first_tmpl = 0; /* user for lineage mask index */
at->first_base = 0; /* used for storing branches */
ai = at->link;
a++;
}
/* build an efficient neighbor/bond list for all atoms */
bonds = os_calloc(n_atom,sizeof(int)); /* list-starts of global bonds indices */
neighbors = os_calloc(n_atom,sizeof(int)); /* list-starts of local atom indices */
bi = pat->bond;
while(bi) {
bd = I->Bond + bi;
bd->cycle = 0; /* initialize */
bd->class = 0;
ai0 = bd->atom[0]; /* global index */
at0 = I->Atom+ai0;
a0 = at0->mark_targ; /* local */
ai1 = bd->atom[1];
at1 = I->Atom+ai1;
a1 = at1->mark_targ;
if((bd->order)&(cH_Double|cH_Triple)) { /* record bond & atoms with Pi bonds */
at0->class |= cH_Pi;
at1->class |= cH_Pi;
bd->class |= cH_Pi;
}
bonds[a0] = ListElemPushInt(&I->Int,bonds[a0],bi); /* enter this bond */
bonds[a1] = ListElemPushInt(&I->Int,bonds[a1],bi);
neighbors[a0] = ListElemPushInt(&I->Int,neighbors[a0],a1); /* cross-enter neighbors */
#ifdef RINGDEBUG
printf(" ring: neighbor of %d is %d (%d)\n",a0,a1,neighbors[a0]);
#endif
neighbors[a1] = ListElemPushInt(&I->Int,neighbors[a1],a0);
#ifdef RINGDEBUG
printf(" ring: neighbor of %d is %d (%d)\n",a1,a0,neighbors[a1]);
#endif
bi = bd->link;
}
/* set up data structures we'll need for ring finding */
lin_mask_size = sizeof(int)*n_atom;
lin_mask_vla = (int**)VLAMalloc(100,lin_mask_size,5,1); /* auto-zero */
/* okay, repeat the ring finding process for each atom in the molecule
(optimize later to exclude dead-ends) */
for(a=0;a<n_atom;a++) {
int expand;
int clean_up;
int sz;
int next;
lin_index = 1;
#ifdef RINGDEBUG
printf("==========\n");
printf(" ring: starting with atom %d\n",a);
ChampMemoryDump(I);
#endif
expand = 0;
clean_up = 0;
sz = 1;
at = atom[a];
clean_up = ListElemPushInt(&I->Int,clean_up,a); /* mark this atom as dirty */
at->mark_tmpl = true;
at->mark_read = true; /* exclude this atom */
expand = ListElemPushInt(&I->Int,expand,a); /* store this atom for the expansion cycle */
vla_check(lin_mask_vla,int*,lin_index); /* assign blank lineage for this atom */
at->first_tmpl = lin_index;
lin_index++;
while(sz<(RING_SEARCH_CUTOFF+1)) { /* while the potential loop is small...*/
next = 0; /* start list of atoms for growth pass */
/* find odd cycles and note growth opportunies */
while(expand) { /* iterate over each atom to expand */
expand = ListElemPopInt(I->Int,expand,&a0);
#ifdef RINGDEBUG
printf(" ring: expand sz %d a0 %d\n",sz,a0);
#endif
ni1 = neighbors[a0];
while(ni1) { /* iterate over each neighbor of that atom */
#ifdef RINGDEBUG
printf(" ring: ni1 %d\n",ni1);
#endif
ni1 = ListElemGetInt(I->Int,ni1,&a1);
#ifdef RINGDEBUG
printf(" ring: neighbor of %d is %d\n",a0,a1);
printf(" ring: neighbor %d mark_read %d\n",a1,atom[a1]->mark_read);
#endif
if(!atom[a1]->mark_read) { /* if that neighbor hasn't yet been covered... */
if(!atom[a1]->first_base) {
next = ListElemPushInt(&I->Int,next,a1); /* store it for growth pass */
}
atom[a1]->first_base = /* and record this entry to it */
ListElemPushInt(&I->Int,atom[a1]->first_base,a0);
if(!atom[a1]->mark_tmpl) {/* mark for later clean_up */
clean_up = ListElemPushInt(&I->Int,clean_up,a1);
atom[a1]->mark_tmpl = true;
}
} else if(sz==atom[a1]->mark_read) { /* ...or if the neighbor has an equal exclusion level... */
#ifdef RINGDEBUG
printf(" ring: checking lineage between %d %d\n",a0,a1);
#endif
if(unrelated_lineage(I,atom[a0]->first_targ, /* unrelated? */
(int*)(((char*)lin_mask_vla)+(atom[a1]->first_tmpl*lin_mask_size))))
{
ring_size = sz*2-1;
#ifdef RINGDEBUG
printf(" ring: #### found odd cycle %d\n",ring_size);
#endif
/* then we have a cycle of size sz*2-1 (odd) */
at0 = atom[a0];
at1 = atom[a1];
merge_lineages(I,
&at0->first_targ,
(int*)(((char*)lin_mask_vla)+(at0->first_tmpl*lin_mask_size)),
&at1->first_targ,
(int*)(((char*)lin_mask_vla)+(at1->first_tmpl*lin_mask_size)));
/* set the atom bits appropriately */
at0->cycle|=num_to_ring[ring_size];
at1->cycle|=num_to_ring[ring_size];
i0 = bonds[a0];
ai0 = atom_idx[a0];
ai1 = atom_idx[a1];
while(i0) {
i0 = ListElemGetInt(I->Int,i0,&bi);
bd = I->Bond+bi;
bd_a0=bd->atom[0];
bd_a1=bd->atom[1];
if(((bd_a0==ai0)&&(bd_a1==ai1))||
((bd_a1==ai0)&&(bd_a0==ai1)))
bd->cycle|=num_to_ring[ring_size];
}
}
}
}
}
/* find even cycles and grow pass */
n_expand = 0;
while(next) { /* iterate over potential branches */
next = ListElemPopInt(I->Int,next,&a1);
#ifdef RINGDEBUG
printf(" ring: next %d\n",a1);
#endif
at1 = atom[a1];
/* see if the base atoms form distinct cycles */
i0 = atom[a1]->first_base;
while(i0) {
a0 = I->Int[i0].value;
i2 = I->Int[i0].link;
while(i2) {
a2 = I->Int[i2].value;
if(unrelated_lineage(I,atom[a0]->first_targ,
(int*)(((char*)lin_mask_vla)+(atom[a2]->first_tmpl*lin_mask_size))))
{
ring_size = sz*2;
at0 = atom[a0];
at2 = atom[a2];
#ifdef RINGDEBUG
printf(" ring: #### found even cycle %d %d %d\n",ring_size,i0,i2);
#endif
at0->cycle|=num_to_ring[ring_size];
at1->cycle|=num_to_ring[ring_size];
at2->cycle|=num_to_ring[ring_size];
i1 = bonds[a1];
ai0 = atom_idx[a0];
ai1 = atom_idx[a1];
ai2 = atom_idx[a2];
while(i1) {
i1 = ListElemGetInt(I->Int,i1,&bi);
bd = I->Bond+bi;
bd_a0=bd->atom[0];
bd_a1=bd->atom[1];
if(((bd_a0==ai0)&&(bd_a1==ai1))||
((bd_a1==ai0)&&(bd_a0==ai1))||
((bd_a0==ai2)&&(bd_a1==ai1))||
((bd_a1==ai2)&&(bd_a0==ai1)))
bd->cycle|=num_to_ring[ring_size];
}
/* we have a cycle of size sz*2 */
/* set atom and bond bits appropriately */
}
i2 = I->Int[i2].link;
}
i0 = I->Int[i0].link;
}
/* allocate a new lineage for this branch atom */
vla_check(lin_mask_vla,int*,lin_index);
at1->first_tmpl = lin_index;
lin_mask = (int*)(((char*)lin_mask_vla)+(lin_mask_size*at1->first_tmpl));
lin_index++;
/* extend this particular lineage with each linking atom's lineage*/
i0 = atom[a1]->first_base;
a0 = I->Int[i0].value;
if(a0!=a) {
lin_mask[a0] = 1; /* add base atom */
at1->first_targ = ListElemPushInt(&I->Int,at1->first_targ,a0);
}
while(i0) { /* and lineage of base atom */
a0 = I->Int[i0].value;
at0 = atom[a0];
at1->first_targ = combine_lineage(I,at0->first_targ,at1->first_targ,lin_mask);
i0 = I->Int[i0].link;
}
expand = ListElemPushInt(&I->Int,expand,a1);
n_expand++;
at1->mark_read = sz+1;
if(!at1->mark_tmpl) {/* mark for later clean_up */
clean_up = ListElemPushInt(&I->Int,clean_up,a1);
at1->mark_tmpl = 1;
}
/* clean up */
ListElemFreeChain(I->Int,atom[a1]->first_base); /* free base atom list */
atom[a1]->first_base = 0;
/* copy the lineage to the new atom */
}
sz++;
if(n_expand<2) {/* on a uniconnected branch, so no point in continuing */
break;
}
}
/* clean-up expansion list */
ListElemFreeChain(I->Int,expand);
expand = 0;
while(clean_up) {
clean_up=ListElemPopInt(I->Int,clean_up,&a0);
at = atom[a0];
/* reset lineage list and mask */
lin_mask = (int*)(((char*)lin_mask_vla)+ (lin_mask_size*at->first_tmpl));
while(at->first_targ) {
at->first_targ = ListElemPopInt(I->Int, at->first_targ,&a1);
lin_mask[a1] = 0;
}
#ifdef RINGDEBUG
for(i0=0;i0<n_atom;i0++) {
if(lin_mask[i0])
printf("lineage mask unclean !\n");
}
#endif
at->mark_read = 0; /* used to mark exclusion */
at->mark_tmpl = 0; /* used to mark as dirty */
at->first_targ = 0; /* used for lineage list */
at->first_tmpl = 0; /* user for lineage mask index */
at->first_base = 0; /* used for storing branches */
}
}
#ifdef RINGDEBUG
for(a=0;a<n_atom;a++) {
i1 = bonds[a];
while(i1) { /* bonds of atom 1 */
i1 = ListElemGetInt(I->Int,i1,&bi1);
bd1 = I->Bond + bi1;
ai1 = bd1->atom[0];
ai2 = bd1->atom[1];
at1 = I->Atom + ai1;
at2 = I->Atom + ai2;
printf("%d %d %x\n",at1->mark_targ+1,at2->mark_targ+1,bd1->cycle);
}
}
#endif
/* now determine which atoms/bonds are aromatic */
/* the following substructures are considered "aromatic" if they occur in rings of 5 or 6
* *=*-*=*
* *=*-*-*=*
* NOTE this is not a chemically accurate definition of aromaticity, just one
* that is easy to program, and which covers most common occurences -- an 80-90% soln.
*
* I'll figure out something better later on...
*/
ai0 = pat->atom;
while(ai0) { /* cycle through each atom */
at0 = I->Atom + ai0;
if(at0->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 0 is cyclic */
i0 = bonds[at0->mark_targ];
while(i0) { /* bonds of atom 0 */
i0 = ListElemGetInt(I->Int,i0,&bi0);
bd0 = I->Bond + bi0;
if((bd0->order==cH_Double)&&(bd0->cycle&(cH_Ring4|cH_Ring5|cH_Ring6))) {
/* bond 0 is double */
ai1 = bd0->atom[0];
if(ai0==ai1) ai1 = bd0->atom[1];
at1 = I->Atom + ai1;
if(at1->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 1 is cyclic */
i1 = bonds[at1->mark_targ];
while(i1) { /* bonds of atom 1 */
i1 = ListElemGetInt(I->Int,i1,&bi1);
bd1 = I->Bond + bi1;
if((bd1->order==cH_Single)&&(bd1->cycle&(cH_Ring4|cH_Ring5|cH_Ring6))) {
/* bond 1 is single */
ai2 = bd1->atom[0];
if(ai1==ai2) ai2 = bd1->atom[1];
at2 = I->Atom + ai2;
if(at2->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 2 is cyclic */
i2 = bonds[at2->mark_targ];
while(i2) {
i2 = ListElemGetInt(I->Int,i2,&bi2);
bd2 = I->Bond + bi2;
if((bd2->order==cH_Double)&&(bd2->cycle&(cH_Ring4|cH_Ring5|cH_Ring6))) {
/* bond 2 is double */
ai3 = bd2->atom[0];
if(ai2==ai3) ai3 = bd2->atom[1];
at3 = I->Atom + ai3;
if(at3->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 3 is cyclic, therefore system is aromatic */
at3->class|=cH_Aromatic|cH_Pi;
at2->class|=cH_Aromatic|cH_Pi;
at1->class|=cH_Aromatic|cH_Pi;
at0->class|=cH_Aromatic|cH_Pi;
bd2->class|=cH_Aromatic|cH_Pi;
bd1->class|=cH_Aromatic|cH_Pi;
bd0->class|=cH_Aromatic|cH_Pi;
}
}
}
i2 = bonds[at2->mark_targ];
while(i2) {
i2 = ListElemGetInt(I->Int,i2,&bi2);
bd2 = I->Bond + bi2;
if((bd2->order==cH_Single)&&(bd2->cycle&(cH_Ring4|cH_Ring5|cH_Ring6))) {
/* bond 2 is single */
ai3 = bd2->atom[0];
if(ai2==ai3) ai3 = bd2->atom[1];
if(ai3!=ai1) {
/* avoid backtracking... */
at3 = I->Atom + ai3;
if(at3->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 3 is cyclic */
i3 = bonds[at3->mark_targ];
while(i3) {
i3 = ListElemGetInt(I->Int,i3,&bi3);
bd3 = I->Bond + bi3;
if((bd3->order==cH_Double)&&(bd3->cycle&(cH_Ring4|cH_Ring5|cH_Ring6))) {
/* bond 3 is double */
ai4 = bd3->atom[0];
if(ai3==ai4) ai4 = bd3->atom[1];
at4 = I->Atom + ai4;
if(at4->cycle&(cH_Ring4|cH_Ring5|cH_Ring6)) {
/* atom 4 is cyclic, therefore system is aromatic */
at4->class|=cH_Aromatic|cH_Pi;
at3->class|=cH_Aromatic|cH_Pi;
at2->class|=cH_Aromatic|cH_Pi;
at1->class|=cH_Aromatic|cH_Pi;
at0->class|=cH_Aromatic|cH_Pi;
bd3->class|=cH_Aromatic|cH_Pi;
bd2->class|=cH_Aromatic|cH_Pi;
bd1->class|=cH_Aromatic|cH_Pi;
bd0->class|=cH_Aromatic|cH_Pi;
}
}
}
}
}
}
}
}
}
}
}
}
}
}
ai0 = at0->link;
}
for(a=0;a<n_atom;a++) { /* default conclusions */
at = atom[a];
if(!at->cycle)
at->cycle=cH_Acyclic;
if(!(at->class&(cH_Aromatic|cH_Aliphatic)))
at->class|=cH_Aliphatic;
}
bi = pat->bond;
while(bi) {
bd = I->Bond + bi;
if(!bd->cycle)
bd->cycle=cH_Acyclic;
if(!(bd->class&(cH_Aromatic|cH_Aliphatic)))
bd->class=cH_Aliphatic;
bi = bd->link;
}
/* now, clean up our efficient bond and neighbor lists */
for(a=0;a<n_atom;a++) {
ListElemFreeChain(I->Int,neighbors[a]);
ListElemFreeChain(I->Int,bonds[a]);
}
mac_free(neighbors);
mac_free(atom);
mac_free(atom_idx);
mac_free(bonds);
/* and free the other data structures */
vla_free(lin_mask_vla);
}
#ifdef RINGDEBUG
ChampPatDump(I,index);
#endif
}