search.c

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

# include "link.h"

int *lig_ctable;		// connection table for lig
int *frag_ctable;		// connection table for frag
int N;				// lig.num_atom;
int M;				// frag.num_atom;
int n;				// lig.nonh_num_atom;
int m;				// frag.nonh_num_atom;


Candidate::Candidate()
{
	int i;

	for(i=0;i<500;i++) ID[i]=0;

	num=0;
}

Candidate::~Candidate()
{
}

void Candidate::Reset(int target)
{
	extern int *lig_ctable;
	extern int *frag_ctable;
	extern int N,M,n,m;	
	int i,j,k;
	int tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;

	target--;

	num=0;

	for(i=0;i<500;i++) ID[i]=0;
		
	tmp1=tmp2=tmp3=0;

	for(i=0;i<m;i++)
		{
		 if(i==target) continue;
		 else if(frag_ctable[target*M+i]==2) tmp1++;
		 else if(frag_ctable[target*M+i]==3) tmp2++;
		 else if(frag_ctable[target*M+i]==4) tmp3++;
		 else continue;
		}

	for(j=0;j<n;j++)
		{
		 if(lig_ctable[j*N+j]!=frag_ctable[target*M+target]) continue;

		 tmp4=tmp5=tmp6=0;
		
		 for(k=0;k<n;k++)
			{
			 if(k==j) continue;
			 else if(lig_ctable[j*N+k]==2) tmp4++;
			 else if(lig_ctable[j*N+k]==3) tmp5++;
			 else if(lig_ctable[j*N+k]==4) tmp6++;
			 else continue;
			}

		 if(tmp4<tmp1) continue;
		 else if(tmp5<tmp2) continue;
		 else if(tmp6<tmp3) continue;
		 else	
			{
			 ID[num]=j+1;
			 num++;
			}
		}

	return;
}

void Candidate::Clean(int wrong)
{
	int i;

	for(i=0;i<num;i++)
		{
		 if(ID[i]!=wrong) continue;
		 else ID[i]=0;
		}

	return;
}

int Candidate::Select()
{
	int i,mark=0;

	for(i=0;i<num;i++)
		{
		 if(ID[i]==0) continue;
		 else {mark=ID[i];break;}
		}

	return mark;
}

int Substructure_Search(const Ligand &lig, const Ligand &frag)
{
	extern int *lig_ctable;
	extern int *frag_ctable;
	extern int N,M,n,m;
	int i,mark;

	if(First_Level_Check(lig,frag)==FALSE) return FALSE;

	// generate the connection table for lig and frag

	N=lig.num_atom; M=frag.num_atom;

	lig_ctable=new int[N*N];
	if(lig_ctable==NULL) Memory_Allocation_Error();

	frag_ctable=new int[M*M];
	if(frag_ctable==NULL) Memory_Allocation_Error();

	n=lig.Get_Connection_Table(lig_ctable);
	m=frag.Get_Connection_Table(frag_ctable);

	// Show_Matrix(n,N,lig_ctable);
	// Show_Matrix(m,M,frag_ctable);

	// now start the searching chain

	Candidate *choice;

	choice=new Candidate[m];
	if(choice==NULL) Memory_Allocation_Error();

	int *chain;
	int p,match;
	int id1,id2;  		// latest exchanged two atoms 
	
	chain=new int[m];
	if(chain==NULL) Memory_Allocation_Error();

	for(i=0;i<m;i++) chain[i]=0;

	p=0; id1=id2=0; 

	while(1)
	{
	 choice[p].Reset(p+1);

	 if(p>0)
		{
		 for(i=0;i<=(p-1);i++) choice[p].Clean(i+1);
		}
	
	 while(1) 
		{
		 mark=choice[p].Select();

		 if(mark==FALSE)
			{
			 if(p==0) {match=FALSE; goto End;}
			 else
			{	 
			 p--;
			 Rearrange_Matrix(N,lig_ctable,chain[p],p+1);
			 continue; 
			}
			}

		 id1=p+1; id2=mark;

		 Rearrange_Matrix(N,lig_ctable,id1,id2);

		 match=Compare_Two_Matrix(p+1,N,lig_ctable,M,frag_ctable);

		 if(match==TRUE) break;
		 else
			{
			 choice[p].Clean(mark);
			 Rearrange_Matrix(N,lig_ctable,id1,id2);
			}
    		} 

	 if((p+1)==m) {match=TRUE; goto End;}
	 else
		{
		 chain[p]=mark;
		 choice[p].Clean(mark);
		 p++;
		}
	}

	End:

	// Show_Matrix(n,N,lig_ctable);
	// Show_Matrix(m,M,frag_ctable);
	
	delete [] lig_ctable; 
	delete [] frag_ctable;
	delete [] choice;
	delete [] chain;

	return match;
}

int First_Level_Check(const Ligand &lig, const Ligand &frag)
{
	int i,j,mark;
	Ligand tmp_lig,tmp_frag;
	int id1,id2,id3,id4;

	// check whether all the atoms in the substructure are fulfilled

	tmp_lig=lig; tmp_frag=frag;

	for(i=0;i<tmp_frag.num_atom;i++)
		{
		 if(tmp_frag.atom[i].type[0]=='H') continue;

		 mark=FALSE;

		 for(j=0;j<tmp_lig.num_atom;j++)
			{
			 if(tmp_lig.atom[j].valid==0) continue;
			 else if(tmp_lig.atom[j].type[0]=='H') continue;
			 else if(strcmp(tmp_frag.atom[i].type,
				 	tmp_lig.atom[j].type)) continue;
			 else {tmp_lig.atom[j].valid=0; mark=TRUE; break;}
			}

		 if(mark==FALSE) return FALSE;
		 else continue;
		}

	// check whether all the bonds in the substructure are fulfilled

	tmp_lig=lig; tmp_frag=frag;

	for(i=0;i<tmp_frag.num_bond;i++)
		{
		 id1=tmp_frag.bond[i].atom_1;
		 id2=tmp_frag.bond[i].atom_2;

		 if(tmp_frag.atom[id1-1].type[0]=='H') continue;
		 else if(tmp_frag.atom[id2-1].type[0]=='H') continue;

		 mark=FALSE;

		 for(j=0;j<tmp_lig.num_bond;j++)
			{
			 if(tmp_lig.bond[j].valid==0) continue;
			 else if(strcmp(tmp_frag.bond[i].type,
					tmp_lig.bond[j].type)) continue;

			 id3=tmp_lig.bond[j].atom_1;
			 id4=tmp_lig.bond[j].atom_2;

			 if(tmp_lig.atom[id3-1].type[0]=='H') continue;
			 else if(tmp_lig.atom[id4-1].type[0]=='H') continue;
			 else if(!strcmp(tmp_lig.atom[id3-1].type,
					 tmp_frag.atom[id1-1].type)&&
				 !strcmp(tmp_lig.atom[id4-1].type,
					tmp_frag.atom[id2-1].type)) mark=TRUE;
			 else if(!strcmp(tmp_lig.atom[id3-1].type,
                                         tmp_frag.atom[id2-1].type)&&
                                 !strcmp(tmp_lig.atom[id4-1].type,
                                        tmp_frag.atom[id1-1].type)) mark=TRUE;
			 else continue;

			 if(mark==TRUE) 
				{
				 tmp_lig.bond[j].valid=0;
				 mark=TRUE; break;
				}
			}

		 if(mark==FALSE) return FALSE;
		 else continue;
		}

	return TRUE;
}

void Rearrange_Matrix(int size, int matrix[], int id1, int id2)
{
	int i,tmp;

	if(id1<=0||id2<=0) return;

	id1--; id2--;

	// first, exchange line id1 and id2

	for(i=0;i<size;i++)
		{
		 tmp=matrix[id1*size+i];
		 matrix[id1*size+i]=matrix[id2*size+i];
		 matrix[id2*size+i]=tmp;
		}

	// second, exchange column id1 and id2

	for(i=0;i<size;i++)
		{
		 tmp=matrix[i*size+id1];
                 matrix[i*size+id1]=matrix[i*size+id2];
                 matrix[i*size+id2]=tmp;
                }

	return;
}

int Compare_Two_Matrix(int subsize, 
		       int size1, int m1[], int size2, int m2[])
{
	int i,j;

	for(i=0;i<subsize;i++)
	for(j=0;j<subsize;j++)
		{
		 if(m1[i*size1+j]==m2[i*size2+j]) continue;
		 else return FALSE;
		}

	return TRUE;
}

void Show_Matrix(int subsize, int size, int matrix[])
{
	int i,j;

	printf("\n");

	for(i=0;i<subsize;i++)
		{
		 for(j=0;j<subsize;j++)
			{
			 printf("%2d ", matrix[i*size+j]);
			}
		 printf("\n");
		}

	return;
}