link.c

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

        else 
		{	
		 result[0]=tmp_result;
		 result[0].valid=TRUE;
		}

	return 1;
}

int Ligand::Join_Two_Molecules(Ligand &lig_1, Ligand &lig_2,
		       	       int seed_1, int root_1, int seed_2, int root_2,
		               char new_bond_type[3], Ligand &result) const
{
	int i,id1,id2,tmp;

	result.Clean_Members();

	// first, delete the two hydrogens 

	lig_1.Delete_An_Atom(seed_1);
	lig_2.Delete_An_Atom(seed_2);

	// second, make the new molecule by adding lig_1 and lig_2

	for(i=0;i<lig_1.num_atom;i++)
		{
		 result.atom[i]=lig_1.atom[i];
		}

	for(i=0;i<lig_1.num_bond;i++)
		{
		 result.bond[i]=lig_1.bond[i];
		}

	for(i=0;i<lig_2.num_atom;i++)
		{
		 tmp=lig_1.num_atom+i;
		 result.atom[tmp]=lig_2.atom[i];
		 result.atom[tmp].id+=lig_1.num_atom;
		}

	for(i=0;i<lig_2.num_bond;i++)
		{
		 tmp=lig_1.num_bond+i;
		 result.bond[tmp]=lig_2.bond[i];
		 result.bond[tmp].id+=lig_1.num_bond;
		 result.bond[tmp].atom_1+=lig_1.num_atom;
		 result.bond[tmp].atom_2+=lig_1.num_atom;
		}

	result.num_atom=lig_1.num_atom+lig_2.num_atom;
	result.num_bond=lig_1.num_bond+lig_2.num_bond;
	strcpy(result.name,"GENERATED");

	// finally, add the new bond

	id1=root_1; id2=root_2+lig_1.num_atom;

	tmp=result.num_bond;
	result.bond[tmp].id=tmp+1;
	result.bond[tmp].atom_1=result.atom[id1-1].id;
	result.bond[tmp].atom_2=result.atom[id2-1].id;
	strcpy(result.bond[tmp].type,new_bond_type);
	result.bond[tmp].valid=1;

	result.num_bond++;

	result.valid=result.Arrange_IDs();	

	if(result.valid==FALSE) return FALSE;
	else return TRUE;
}

int Ligand::Find_Root_Of_A_Hydrogen(int id, Group &root) const
// id is the ID of the hydrogen
{
	int id_root;
	
	if(atom[id-1].type[0]!='H') return FALSE;  // not a hydrogen!

	// first, find the root atom

	id_root=atom[id-1].neib[0];

	if(id_root==0) return FALSE;

	// second, find the neighbors of the root atom

	root=Find_A_Group(id_root);

	if(root.valid==0) return FALSE;
	else return TRUE;
}

int Ligand::Select_A_Seed_Hydrogen() const
{
        int i,tmp;
	int num,total,final,*seed_list;
	struct Wheel
		{
		 int min;
		 int max;
		} *wheel;

        // first, count all the possible seed hydrogens

	seed_list=new int[num_atom];
	if(seed_list==NULL) Memory_Allocation_Error();

	num=Make_Seed_List(seed_list);
	if(num==0) {delete [] seed_list; return FALSE;}

	// then, set the wheel acoording to their grades

	wheel=new Wheel[num];
	if(wheel==NULL) Memory_Allocation_Error();

	total=0;

	for(i=0;i<num;i++)
		{
		 wheel[i].min=total;
		 total+=Pow(2,atom[seed_list[i]-1].valid);
		 wheel[i].max=total;
		}

        // then, randomly select one from the wheel 

	tmp=(int)(drand48()*total);

	final=0;

	for(i=0;i<num;i++)
		{
		 if(tmp<wheel[i].min) continue;
		 else if(tmp>=wheel[i].max) continue;
		 else {final=seed_list[i]; break;}
		}

	delete [] wheel; delete [] seed_list;

	return final;
}

int Ligand::Make_Seed_List(int seed_id[]) const
{
	int i,num;
	
	for(i=0;i<num_atom;i++) seed_id[i]=0;

	num=0;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(strcmp(atom[i].type,"H.spc")) continue;
		 else {seed_id[num]=atom[i].id; num++;}
		}

	return num;
}

int Ligand::Record_Sameness_Check(Linking_Record record_1,
				  Linking_Record record_2) const
{
	int i,mark;

	if(record_1.num_link!=record_2.num_link) return FALSE;

	mark=TRUE;

	for(i=0;i<record_1.num_link;i++)
		{
		 if(record_1.id1[i]!=record_2.id1[i]) 
			{
			 mark=FALSE; break;
			}
		 else if(record_1.id2[i]!=record_2.id2[i])
			{
			 mark=FALSE; break;
			}
		 else continue;
		}

	return mark;
}

int Ligand::Atom_Overlap_Check(Atom atom_1, Atom atom_2) const
{
	if(Distance(atom_1.coor,atom_2.coor)<ATOM_OVERLAP_RANGE) return TRUE;
	else return FALSE;
}

int Ligand::Bridging(int id1, int id2)
{
	extern ForceField *ff;
	Atom seed_1,seed_2;
        Atom result_1,result_2,result;
        Group root_1, root_2;
        int i;
        float tmp1,tmp2,tmp3,tmp4,d;
	float v1[3],v2[3],v3[3],v4[3],angle;

        // find the two seed hydrogens and their roots

        seed_1=atom[id1-1]; seed_2=atom[id2-1];

	if(Find_Root_Of_A_Hydrogen(id1,root_1)==FALSE) return FALSE;
	if(Find_Root_Of_A_Hydrogen(id2,root_2)==FALSE) return FALSE;

	if(Atom_Overlap_Check(root_1.center,root_2.center)==TRUE) return FALSE;

	// check the angle of the two vectors

	for(i=0;i<3;i++)
		{
		 v1[i]=root_1.center.coor[i]-seed_1.coor[i];
		 v2[i]=root_2.center.coor[i]-seed_2.coor[i];
		}

	angle=Angle_Of_Two_Vectors(v1,v2);

	if(angle>(109.4+ANGLE_OVERLAP_RANGE)) return FALSE;
	else if(angle<(109.4-ANGLE_OVERLAP_RANGE)) return FALSE;

	// calculate the position of hypothetic atoms

        tmp1=Distance(seed_1.coor,root_1.center.coor);
	tmp2=Distance(seed_2.coor,root_2.center.coor);
        tmp3=ff->Get_Bond_Length(root_1.center.type,"C.3","1");
	tmp4=ff->Get_Bond_Length(root_2.center.type,"C.3","1");

        for(i=0;i<3;i++)
        {
         result_1.coor[i]=root_1.center.coor[i];
         result_1.coor[i]+=((seed_1.coor[i]-root_1.center.coor[i])*tmp3/tmp1);
	 result_2.coor[i]=root_2.center.coor[i];
	 result_2.coor[i]+=((seed_2.coor[i]-root_2.center.coor[i])*tmp4/tmp2);
        }

	if(Atom_Overlap_Check(result_1,result_2)==FALSE) return FALSE;

        // first, delete the two hydrogens

        Delete_An_Atom(id1); Delete_An_Atom(id2);

        // add the new bridging carbon, stored in result

	result.id=num_atom+1;
	strcpy(result.name,"C");
        strcpy(result.type,"C.3");
        result.weight=12;
        result.valid=seed_2.valid;
	for(i=0;i<3;i++) 
	 	result.coor[i]=(result_1.coor[i]+result_2.coor[i])/2.0;

	atom[num_atom]=result;

        num_atom++;

        bond[num_bond].id=num_bond+1;
        bond[num_bond].atom_1=root_1.center.id;
        bond[num_bond].atom_2=result.id;
        bond[num_bond].valid=1;
        strcpy(bond[num_bond].type,"1");

        bond[num_bond+1].id=num_bond+2;
        bond[num_bond+1].atom_1=root_2.center.id;
        bond[num_bond+1].atom_2=result.id;
        bond[num_bond+1].valid=1;
        strcpy(bond[num_bond+1].type,"1");

        num_bond+=2;

        // add the two new hydrogens, stored in result_1 and result_2

        for(i=0;i<=2;i++)
                {
                 v1[i]=root_1.center.coor[i]-result.coor[i];
                 v2[i]=root_2.center.coor[i]-result.coor[i];
                }

        Unify_Vector(v1); Unify_Vector(v2);

        Get_Sp3_Coordinates(v1,v2,v3,v4);

        d=ff->Get_Bond_Length("C.3","H","1");

        for(i=0;i<=2;i++)
                {
                 result_1.coor[i]=result.coor[i]+v3[i]*d;
                 result_2.coor[i]=result.coor[i]+v4[i]*d;
                }

        result_1.id=num_atom+1; 
        strcpy(result_1.name,"H"); 
	strcpy(result_1.type,"H.spc");
        result_1.weight=1;
	result_1.valid=seed_2.valid;

	result_2.id=num_atom+2;
	strcpy(result_2.name,"H");
	strcpy(result_2.type,"H.spc");
	result_2.weight=1;
	result_2.valid=seed_2.valid;

        atom[num_atom]=result_1;
	atom[num_atom+1]=result_2;
	
	num_atom+=2;

        bond[num_bond].id=num_bond+1;
        bond[num_bond].atom_1=result.id;
        bond[num_bond].atom_2=result_1.id;
        bond[num_bond].valid=1;
        strcpy(bond[num_bond].type,"1");

        bond[num_bond+1].id=num_bond+2;
        bond[num_bond+1].atom_1=result.id;
        bond[num_bond+1].atom_2=result_2.id;
        bond[num_bond+1].valid=1;
        strcpy(bond[num_bond+1].type,"1");

	num_bond+=2;

/*
        printf("Bridging happens between %d and %d\n", 
		root_1.center.id, root_2.center.id);
*/

        return TRUE;
}

int Ligand::Joining(int id1, int id2)
{
        extern ForceField *ff;
        Atom seed_1,seed_2;
        Atom result_1,result_2;
        Group root_1, root_2;
        int i;
        float tmp1,tmp2,tmp3;
        Bond new_bond;

        // find the two seed hydrogens and their roots

        seed_1=atom[id1-1]; seed_2=atom[id2-1];

        if(Find_Root_Of_A_Hydrogen(id1,root_1)==FALSE) return FALSE;
        if(Find_Root_Of_A_Hydrogen(id2,root_2)==FALSE) return FALSE;

	if(Connection_1_2_Check(root_1.center.id,root_2.center.id)==TRUE) 
		return FALSE;

	strcpy(new_bond.type,New_Bond_Viability_Check(root_1,root_2));

	if(!strcmp(new_bond.type,"no")) return FALSE;	
        else
                {
		 new_bond.id=num_bond+1;
                 new_bond.atom_1=root_1.center.id;
                 new_bond.atom_2=root_2.center.id;
                 new_bond.length=ff->Get_Bond_Length(root_1.center.type,
                                                     root_2.center.type,
                                                     new_bond.type);
		 new_bond.valid=1;
                }

        tmp1=Distance(seed_1.coor,root_1.center.coor);
	tmp2=Distance(seed_2.coor,root_2.center.coor);
        tmp3=new_bond.length;

        for(i=0;i<3;i++)
        {
         result_1.coor[i]=root_1.center.coor[i];
         result_1.coor[i]+=((seed_1.coor[i]-root_1.center.coor[i])*tmp3/tmp1);
         result_2.coor[i]=root_2.center.coor[i];
         result_2.coor[i]+=((seed_2.coor[i]-root_2.center.coor[i])*tmp3/tmp2);
        }

	if(Atom_Overlap_Check(result_1,root_2.center)==FALSE) return FALSE;
	if(Atom_Overlap_Check(result_2,root_1.center)==FALSE) return FALSE;

        // the new bond may need to change the atom type of root_1 or root_2

        atom[root_1.center.id-1]=root_1.center;
        atom[root_2.center.id-1]=root_2.center;

        // delete the two hydrogens

        Delete_An_Atom(id1); Delete_An_Atom(id2);

	// add the new bond

	bond[num_bond]=new_bond;

        num_bond++;

/*
        printf("Linking happens between %d and %d\n", 
		root_1.center.id, root_2.center.id);
*/

        return TRUE;
}

int Ligand::Fusing(int id1, int id2)
{
        Group root_1, root_2;
	Bond new_bond;
        int i,mark;
	int delete_mode_1, delete_mode_2;
	char type_1[5],type_2[5];

        struct Bonding_Orientation
                {
                 int valid;
                 int id1;
                 int id2;
                } bonding_1[4], bonding_2[4];

	for(i=0;i<4;i++)
		{
		 bonding_1[i].valid=bonding_2[i].valid=0;
		 bonding_1[i].id1=bonding_2[i].id1=0;
		 bonding_1[i].id2=bonding_2[i].id2=0;
		}

	root_1=Find_A_Group(id1); root_2=Find_A_Group(id2);

	// check whether their hybridization types match or not
	// in fact, only sp3-sp3 matching is allowed

	strcpy(type_1,Get_Atom_Hybridizing_Type(root_1.center));
	strcpy(type_2,Get_Atom_Hybridizing_Type(root_2.center));

	if(!strcmp(type_1,"none")||!strcmp(type_2,"none")) return FALSE;
	else if(!strcmp(type_1,"sp2")||!strcmp(type_2,"sp2")) return FALSE;
	else if(!strcmp(type_1,"sp")||!strcmp(type_2,"sp")) return FALSE;
	else if(strcmp(type_1,type_2)) return FALSE;

	// check whether the bonding orientations are match or not

	if((root_1.num_nonh+root_2.num_nonh)>4) return FALSE;

	delete_mode_1=TRUE;

	for(i=0;i<root_1.num_nonh;i++)
		{
		 bonding_1[i].id1=root_1.neib[i].id;

		 mark=Bonding_Match_Check(root_1.neib[i],root_2);

		 if(mark==0) 
			{
			 delete_mode_1=FALSE;
			 break;
			}
		 else 
			{
			 bonding_1[i].id2=mark; 
			 bonding_1[i].valid=1;
			}
		}

	delete_mode_2=TRUE;

        for(i=0;i<root_2.num_nonh;i++)
                {
                 bonding_2[i].id1=root_2.neib[i].id;

                 mark=Bonding_Match_Check(root_2.neib[i],root_1);

                 if(mark==0)
			{
			 delete_mode_2=FALSE;
			 break;
			}
                 else 
                        {
                         bonding_2[i].id2=mark;
                         bonding_2[i].valid=1;
                        }
                }

	// determine which atom should be deleted

	if((delete_mode_1==FALSE)&&(delete_mode_2==FALSE)) return FALSE;
	else if((delete_mode_1==TRUE)&&(delete_mode_2==TRUE)) 
		{
		 if(root_1.num_nonh>root_2.num_nonh) delete_mode_1=FALSE;
		 else if(root_1.num_nonh<root_2.num_nonh) delete_mode_2=FALSE;
		 else if(root_1.center.weight>=root_2.center.weight) delete_mode_1=FALSE;
		 else delete_mode_2=FALSE;
		}

	if((delete_mode_1==TRUE)&&(delete_mode_2==FALSE))
	{
	 // delete root_1

	 for(i=0;i<root_1.num_neib;i++)
		{
		 if(root_1.neib[i].type[0]!='H') Delete_An_Atom(bonding_1[i].id2);
		 else Delete_An_Atom(root_1.neib[i].id);
		}

	 Delete_An_Atom(root_1.center.id);

	 for(i=0;i<root_1.num_nonh;i++)
		{
		 new_bond.id=num_bond+1;
		 new_bond.atom_1=root_1.neib[i].id;
		 new_bond.atom_2=root_2.center.id;
		 strcpy(new_bond.type,"1");