ligand.c

药物开发中的基于结构的从头设计代码

		 d=Distance(atom[i].coor,atom[j].coor);	
		 d0=atom[i].r+atom[j].r;

		 if(d>(d0-VDW_BUMP_RANGE)) continue;	// no bump;
		 else if(Connection_1_2_Check(atom[i].id, atom[j].id)==TRUE) 
			continue; 
		 else if(Connection_1_3_Check(atom[i].id, atom[j].id)==TRUE) 
			continue;
		 else if(Connection_1_4_Check(atom[i].id, atom[j].id)==TRUE) 
			continue;
		 else return TRUE; 
		}

	return FALSE;
}

int Ligand::Connection_1_2_Check(int id1, int id2) const
{
	int i;

	if(id1==id2) return FALSE;

	for(i=0;i<atom[id2-1].num_neib;i++)
		{
		 if(id1==atom[id2-1].neib[i]) return TRUE; 
		 else continue;
		}
	
	return FALSE;
}

int Ligand::Connection_1_3_Check(int id1, int id2) const
{
	int i,j;

	if(id1==id2) return FALSE;
	if(Connection_1_2_Check(id1,id2)==TRUE) return FALSE;

	id1--; id2--;

	for(i=0;i<atom[id1].num_neib;i++)
	for(j=0;j<atom[id2].num_neib;j++)
		{
		 if(atom[id1].neib[i]==atom[id2].neib[j]) return TRUE;
		 else continue;
		}

	return FALSE;
}

int Ligand::Connection_1_4_Check(int id1, int id2) const
{
	int i,j;

	if(id1==id2) return FALSE;
	if(Connection_1_2_Check(id1,id2)==TRUE) return FALSE;
	if(Connection_1_3_Check(id1,id2)==TRUE) return FALSE;

	id1--; id2--;

	for(i=0;i<atom[id1].num_neib;i++)
	for(j=0;j<atom[id2].num_neib;j++)
		{
		 if(Connection_1_2_Check(atom[id1].neib[i],
		    atom[id2].neib[j])==TRUE) return TRUE;
		 else continue;
		}

	return FALSE;
}

int Ligand::Connection_1_5_Check(int id1, int id2) const
{
        int i,j;

        if(id1==id2) return FALSE;
        if(Connection_1_2_Check(id1,id2)==TRUE) return FALSE;
        if(Connection_1_3_Check(id1,id2)==TRUE) return FALSE;
        if(Connection_1_4_Check(id1,id2)==TRUE) return FALSE;

        id1--; id2--;

        for(i=0;i<atom[id1].num_neib;i++)
        for(j=0;j<atom[id2].num_neib;j++)
                {
                 if(Connection_1_3_Check(atom[id1].neib[i],
                    atom[id2].neib[j])==TRUE) return TRUE;
                 else continue;
                }

        return FALSE;
}

int Ligand::Get_Connection_Table(int ctable[]) const
{
	int i,j,k,tmp1,tmp2;
	int size;

	// first, initialize the connection table

	size=num_atom;

	for(i=0;i<size;i++)
	for(j=0;j<size;j++)
		{
		 if(i!=j) ctable[i*size+j]=0;
                 else if(!strcmp(atom[i].type,"H")) ctable[i*size+j]=1;
                 else if(!strcmp(atom[i].type,"H.spc")) ctable[i*size+j]=1;
                 else if(!strcmp(atom[i].type,"C.1")) ctable[i*size+j]=2;
                 else if(!strcmp(atom[i].type,"C.2")) ctable[i*size+j]=3;
                 else if(!strcmp(atom[i].type,"C.3")) ctable[i*size+j]=4;
                 else if(!strcmp(atom[i].type,"C.ar")) ctable[i*size+j]=5;
                 else if(!strcmp(atom[i].type,"C.cat")) ctable[i*size+j]=6;
                 else if(!strcmp(atom[i].type,"N.1")) ctable[i*size+j]=7;
                 else if(!strcmp(atom[i].type,"N.2")) ctable[i*size+j]=8;
                 else if(!strcmp(atom[i].type,"N.3")) ctable[i*size+j]=9;
                 else if(!strcmp(atom[i].type,"N.4")) ctable[i*size+j]=9;
                 else if(!strcmp(atom[i].type,"N.am")) ctable[i*size+j]=10;
                 else if(!strcmp(atom[i].type,"N.ar")) ctable[i*size+j]=11;
                 else if(!strcmp(atom[i].type,"N.pl3")) ctable[i*size+j]=9;
                 else if(!strcmp(atom[i].type,"O.2")) ctable[i*size+j]=13;
                 else if(!strcmp(atom[i].type,"O.3")) ctable[i*size+j]=14;
                 else if(!strcmp(atom[i].type,"O.co2")) ctable[i*size+j]=13;
                 else if(!strcmp(atom[i].type,"P.3")) ctable[i*size+j]=15;
                 else if(!strcmp(atom[i].type,"S.2")) ctable[i*size+j]=16;
                 else if(!strcmp(atom[i].type,"S.3")) ctable[i*size+j]=17;
                 else if(!strcmp(atom[i].type,"S.o")) ctable[i*size+j]=18;
                 else if(!strcmp(atom[i].type,"S.o2")) ctable[i*size+j]=19;
                 else if(!strcmp(atom[i].type,"F")) ctable[i*size+j]=20;
                 else if(!strcmp(atom[i].type,"Cl")) ctable[i*size+j]=21;
                 else if(!strcmp(atom[i].type,"Br")) ctable[i*size+j]=22;
                 else if(!strcmp(atom[i].type,"I")) ctable[i*size+j]=23;
                 else ctable[i*size+j]=0;	// Unknown atom type
		}

	// fill in the connection table

	for(i=0;i<size-1;i++)
	for(j=i+1;j<size;j++)
		{
		 if(atom[i].type[0]=='H'||atom[j].type[0]=='H') continue;
		 else if(Connection_1_2_Check(atom[i].id,atom[j].id)==TRUE)
			{
			 ctable[i*size+j]=2;
			 ctable[j*size+i]=2;
			}
		 else if(Connection_1_3_Check(atom[i].id,atom[j].id)==TRUE)
			{
                         ctable[i*size+j]=3;
                         ctable[j*size+i]=3;
                        }
                 else if(Connection_1_4_Check(atom[i].id,atom[j].id)==TRUE)
                        {
			 ctable[i*size+j]=4;
                         ctable[j*size+i]=4;
                        }
		 else continue;
		}

	// rank the connection table

	for(i=0;i<size-1;i++)
	for(j=i+1;j<size;j++)
		{
		 if(ctable[i*size+i]>ctable[j*size+j]) continue;
		 else if(ctable[i*size+i]<ctable[j*size+j])
			{
			 Rearrange_Matrix(size,ctable,i+1,j+1);
			 continue;
			}
		 else if(ctable[i*size+i]==ctable[j*size+j])  // atom type same
			{
			 // check who has more connections

			 tmp1=0;
			 for(k=0;k<size;k++)
				{
				 if(k==i) continue;
				 else if(ctable[i*size+k]==2) tmp1+=10000;
				 else if(ctable[i*size+k]==3) tmp1+=100;
				 else if(ctable[i*size+k]==4) tmp1++;
				 else continue;
				}

			 tmp2=0;
			 for(k=0;k<size;k++)
                                {
                                 if(k==j) continue;
                                 else if(ctable[j*size+k]==2) tmp2+=10000;
				 else if(ctable[j*size+k]==3) tmp2+=100;
				 else if(ctable[j*size+k]==4) tmp2++;
				 else continue;
                                }

			 if(tmp1>tmp2) continue;
			 else if(tmp1<tmp2)
				{
				 Rearrange_Matrix(size,ctable,i+1,j+1);
				 continue;
				}
			 else continue;
			} 
		 else continue;
		}

	// finally, find out the number of non-H atoms

	for(i=0;i<size;i++)
		{
		 if(ctable[i*size+i]>1) continue;
		 else {size=i; break;}
		}	

	return size;
}

void Ligand::Detect_Rings()
{
	int i,j,mark,count;

	// first of all, check whether there is any ring in the molecule

	if((num_bond-num_atom+num_part)==0)
		{
		 for(i=0;i<num_atom;i++) atom[i].ring=0;
		 for(i=0;i<num_bond;i++) bond[i].ring=0;
		 return;
		}

	// detect terminal atoms

	for(i=0;i<num_atom;i++) 	
		{
		 if(atom[i].num_nonh==1) atom[i].ring=0;
		 else atom[i].ring=-1;	// ambiguous at present
		}

	// collapse the structure to eliminate branches

	do		
		{
		 mark=0;
		 for(i=0;i<num_atom;i++)
			{
			 if(atom[i].ring!=-1) continue;

			 count=0;	// count possible ring neighbors	
		 	 for(j=0;j<atom[i].num_neib;j++)
			{
			 if(atom[atom[i].neib[j]-1].ring==0) continue;
			 else count++;
			}

			 if(count<=1) {atom[i].ring=0; mark++;}
			}
		}
	while(mark);

	// detect branching bonds

	mark=0;

	for(i=0;i<num_bond;i++) 
		{
		 if(atom[bond[i].atom_1-1].ring==0) bond[i].ring=0;
		 else if(atom[bond[i].atom_2-1].ring==0) bond[i].ring=0;
		 else {bond[i].ring=-1; mark++;}	// ambiguous
		}

	if(mark==0) return; 	// finished, no ring detected

	// check all the ambiguous bonds
	
	for(i=0;i<num_bond;i++)	
		{
		 if(bond[i].ring!=-1) continue;
		 else bond[i].ring=Judge_If_An_Bond_In_Ring(bond[i].id);
		}

	// find all the atoms in rings

	for(i=0;i<num_bond;i++) 
		{
		 if(bond[i].ring!=1) continue;
		 else
			{
			 atom[bond[i].atom_1-1].ring=1;
			 atom[bond[i].atom_2-1].ring=1;
			}
		}

	for(i=0;i<num_atom;i++)	// define all the other atoms
		{
		 if(atom[i].ring!=-1) continue;
		 else atom[i].ring=0;
		}
	
	return;
}

void Ligand::Reset_Choices(int id,int choice[]) const
// id is ID of the bond
{
	int i,tmp;

	for(i=0;i<10;i++) choice[i]=0;

	if(bond[id-1].ring!=0) 
		{
		 for(i=0;i<bond[id-1].num_neib;i++)
			{
			 tmp=bond[id-1].neib[i];
			 if(bond[tmp-1].ring==0) choice[i]=0;
			 else choice[i]=tmp; 
			}
		}

	return;
}

void Ligand::Check_Choices(int head, int choice[]) const
// head: ID of the head atom of this bond
{
	int i;

	for(i=0;i<10;i++)
		{
		 if(choice[i]==0) continue;
		 else if(choice[i]>num_bond) choice[i]=0;
		 else if(bond[choice[i]-1].ring==0) choice[i]=0;
		 else if(bond[choice[i]-1].atom_1==head) choice[i]=0;
		 else if(bond[choice[i]-1].atom_2==head) choice[i]=0;
		 else continue;
		}

	return;
}

void Ligand::Clean_Choices(int wrong, int choice[]) const
// wrong: the ID of the unwanted bond
{
        int i;

        for(i=0;i<10;i++)
                {
                 if(choice[i]==0) continue;
                 else if(choice[i]>num_bond) choice[i]=0;
                 else if(choice[i]==wrong) choice[i]=0;
                 else continue;
                }

        return;
}

int Ligand::Judge_If_An_Bond_In_Ring(int bond_id) const
// this is a bond-based algorithm to detect rings in a molecule.
{
	int i,j,tmp;
	int p,next;
	int atom_chain[500], bond_chain[500];
	int choice[500][10];

	for(i=0;i<500;i++)
		{
		 atom_chain[i]=bond_chain[i]=0;
		 for(j=0;j<10;j++) choice[i][j]=0;
		}

	for(i=0;i<num_bond;i++) Reset_Choices(bond[i].id,choice[i]);

	// push the bond into the searching path

	p=0;	// pointer of the searching path

	bond_chain[p]=bond_id; 
	atom_chain[p]=bond[bond_id-1].atom_1;	// head of this bond

	for(;;)
	{
		// find the next bond to enlong the searching path
		do
		{
		 tmp=bond_chain[p]-1;

		 Check_Choices(atom_chain[p], choice[tmp]);
		
		 if(p>0) Clean_Choices(bond_chain[p-1], choice[tmp]);

		 next=FALSE;

		 for(i=0;i<10;i++)	// select the next bond
			{
			 if(choice[tmp][i]==0) continue;
			 else {next=choice[tmp][i]; break;}
			}

		 if(next==FALSE)	// no choice is available
			{
			 if(p==0) return FALSE;	// this bond is not in ring.
			 else		// pop out the bond of the path
				{
				 Reset_Choices(bond_chain[p],choice[tmp]);
				 bond_chain[p]=0;
				 atom_chain[p]=0;
				 p--;
				 continue;
				}
			}

                 for(i=0;i<=p;i++)
                        {
                         if(next!=bond_chain[i]) continue;
                         else           // this is a wrong choice
                                {
                                 Clean_Choices(next,choice[tmp]);
                                 next=FALSE; break;
                                }
                        }

		 if(next!=FALSE) 
			{
			 if(Two_Bonds_Connection_Check(bond[bond_chain[0]-1],
			    bond[next-1])==atom_chain[0])    // find a ring
				{
				 // puts("Ring detected:");
				 for(i=0;i<=p;i++) 
					{
					 bond[bond_chain[i]-1].ring=1;
					/*
					 printf("Atom %d - %d\n",
					 bond[bond_chain[i]-1].atom_1,
					 bond[bond_chain[i]-1].atom_2);
					*/
					}
				 bond[next-1].ring=1;
				/*
				 printf("Atom %d - %d\n\n",
					bond[next-1].atom_1,
					bond[next-1].atom_2);
			        */
				 return TRUE;
				}
			}
			 
		} while(next==FALSE);

		Clean_Choices(next,choice[tmp]);
		tmp=Two_Bonds_Connection_Check(bond[tmp],bond[next-1]);
		p++;
		bond_chain[p]=next;
		atom_chain[p]=tmp;
	}
}