grow.c

来自「药物开发中的基于结构的从头设计代码」· C语言 代码 · 共 1,238 行 · 第 1/3 页

                                 (tmp_frag[i],tmp_frag[j])==FALSE) continue;
                         else 
				{
				 profile[j].valid=0;
				 profile[j].grow=0;
				 profile[j].link=0;
				}
                        }
                }

	// now make new molecules

	for(i=0;i<24;i++)
		{
		 if(profile[i].valid==FALSE) continue;
		 else
			{
		 	 mark=Join_Two_Molecules(core,tmp_frag[i],
					 	 id5,id2,id6,id3,new_bond.type,
					 	 tmp_result);
		 	 if(mark==FALSE) return FALSE;
		 	 else result[i]=tmp_result;
			}
		}

	// try to make linkings

	 Linking_Record record[24];

	 for(i=0;i<24;i++)
		{
		 if(profile[i].link==FALSE) continue;

		 tmp_result<=result[i];	// notice this operator!

		 if(tmp_result.Link_Structure(record[i])==FALSE) 
			{
			 profile[i].link=FALSE;
			 profile[i].valid=FALSE;
			}
		 else
			{
			 result[i]=tmp_result;
			}
		}

	// find the best conformer among linked structures 

	 for(i=0;i<23;i++)
		{
		 if(profile[i].link==FALSE) continue;
		 for(j=i+1;j<24;j++)
			{
			 if(profile[j].link==FALSE) continue;
			 if(Record_Sameness_Check(record[i],record[j])==TRUE)
				{
				 if(record[i].rmsd>record[j].rmsd) 
					{
					 profile[i].link=FALSE;
					 profile[i].valid=FALSE;
					}
				 else 
					{
					 profile[j].link=FALSE;
					 profile[j].valid=FALSE;
					}
				}
			 else continue;
			}
		}

	// finally, count the generated molecules

	mark=0;

	for(i=0;i<24;i++)
		{
		 if(profile[i].valid) {result[i].valid=TRUE; mark++;}
		 else result[i].valid=FALSE;
		}

/*
	for(i=0;i<24;i++)
		{
		 printf("Angle: %5.1f\tEnergy: %7.2f\tGrow: %d\tLink: %d\t",
                        profile[i].angle, profile[i].e, 
			profile[i].grow, profile[i].link);
                 printf("Result: %d\n", result[i].valid);
		}
*/

	return mark;	// return the number of molecules generated


Improper_Torsion:

	mark=Join_Two_Molecules(core,frag,
                                id5,id2,id6,id3,new_bond.type,
                                tmp_result);

        if(mark==FALSE) return FALSE;
        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;