link.c

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

# include "link.h"

int Ligand::Generate_New_Structure(int seed_1, Ligand &frag, int seed_2, Ligand result[24]) const
// growing frag/seed_2 with ligand/seed_1, results stored in result[]
{ 
	extern Parameter *parm;
	extern FragLibrary *forbidlib;
	int i,mark;
	int old_ring, new_ring;
	Ligand core;

	core=*this;

	for(i=0;i<24;i++) result[i].Clean_Members();

	// growing & linking 

	if(drand48()<parm->grow_ratio)
		{
		 mark=Grow_A_Fragment(core,frag,seed_1,seed_2,result);
		 if(mark==FALSE) return FALSE;
		 // else printf("%d molecules generated\n", mark);
		}
	else return FALSE;

	// mutation

	if(drand48()<parm->mutate_ratio)
		{
		 for(i=0;i<24;i++) 
			{
		 	 if(result[i].valid==0) continue;
		 	 else 
				{
				 result[i].Assign_Atom_Parameters();
				 result[i].Rational_Mutate();
				}
			}
		}

	// count the newly formed rings 

	old_ring=new_ring=0;
	old_ring+=(this->num_bond-this->num_atom+this->num_part);
	old_ring+=(frag.num_bond-frag.num_atom+frag.num_part);

	for(i=0;i<24;i++)
		{
		 if(result[i].valid==0) continue;
		 result[i].Assign_Atom_Parameters();
		 new_ring=result[i].num_bond-result[i].num_atom+result[i].num_part;
		 result[i].num_fused_ring=new_ring-old_ring;
		}

	// check whether bump with the receptor 

	for(i=0;i<24;i++)
		{
		 if(result[i].valid==0) continue;
		 else if(result[i].Outer_Collision_Check()==FALSE) 
			{
			 // puts("Outer collision!");
			 continue; 
			}
		 else result[i].valid=0;
		}

	// check whether the newly generated molecules are reasonable

	for(i=0;i<24;i++)
		{
		 if(result[i].valid==0) continue;
		 else if(result[i].Inner_Collision_Check()==FALSE) 
			{
			 // puts("Inner collision!");
			 continue;
			}
		 else result[i].valid=0;
		}

	// scoring

	for(i=0;i<24;i++)
		{
		 if(result[i].valid==0) continue;
		 else 
			{
			 result[i].Get_Molecular_Properties();
			 result[i].Calculate_Binding_Score();		 
			 result[i].Calculate_Chemical_Score();
			}
		}

	// forbidden substructure checking and
	// counting the survived children molecules

	mark=0;

	for(i=0;i<24;i++)
		{
		 if(result[i].valid==0) continue;
		 else if(!strcmp(parm->apply_forbidden_check,"YES"))
			{
			 if(forbidlib->Map_A_Molecule(result[i])==FALSE) mark++;
			 else result[i].valid=0; 
			}
		 else mark++;
		}

/*
        result[0].Write_Out_Mol2("test000.mol2");
        result[1].Write_Out_Mol2("test015.mol2");
        result[2].Write_Out_Mol2("test030.mol2");
        result[3].Write_Out_Mol2("test045.mol2");
        result[4].Write_Out_Mol2("test060.mol2");
        result[5].Write_Out_Mol2("test075.mol2");
        result[6].Write_Out_Mol2("test090.mol2");
        result[7].Write_Out_Mol2("test105.mol2");
        result[8].Write_Out_Mol2("test120.mol2");
        result[9].Write_Out_Mol2("test135.mol2");
        result[10].Write_Out_Mol2("test150.mol2");
        result[11].Write_Out_Mol2("test165.mol2");
        result[12].Write_Out_Mol2("test180.mol2");
        result[13].Write_Out_Mol2("test195.mol2");
        result[14].Write_Out_Mol2("test210.mol2");
        result[15].Write_Out_Mol2("test225.mol2");
        result[16].Write_Out_Mol2("test240.mol2");
        result[17].Write_Out_Mol2("test255.mol2");
        result[18].Write_Out_Mol2("test270.mol2");
        result[19].Write_Out_Mol2("test285.mol2");
        result[20].Write_Out_Mol2("test300.mol2");
        result[21].Write_Out_Mol2("test315.mol2");
        result[22].Write_Out_Mol2("test330.mol2");
        result[23].Write_Out_Mol2("test345.mol2");
*/

	return mark;
}

int Ligand::Grow_A_Fragment(Ligand &core, Ligand &frag,
			    int seed_1, int seed_2, Ligand result[24]) const
/* 
   'seed_1' is the ID of the seed hydrogen on the core; 
   'seed_2' is the ID of the seed hydrogen on the fragment;
   'result[24]' will return the new molecules.
*/
{
	extern Parameter *parm;
	extern ForceField *ff;
	Group root_1,root_2;
	Bond new_bond;
	Torsion new_torsion;
	int i,j,k,l,mark,mark1,mark2;
	int id1,id2,id3,id4,id5,id6;
	float tmp1,tmp2,tmp3,angle;
	float move[3],v1[3],v2[3],anchor[3],axis[3];

/*
   set up atom "grades" for the ligand and the fragment. this builds a "tree" 
   of the molecular structure and it will help to select the seed hydrogens 
   during the later building-up process.
*/

	mark=core.atom[seed_1-1].valid;

	if(mark>=2)
		{
		 for(i=0;i<frag.num_atom;i++) frag.atom[i].valid=mark-1;
		}
	else if(mark==1)
		{
		 for(i=0;i<core.num_atom;i++) core.atom[i].valid++;
		}
	else 
		{
		 return FALSE;
		}

	// find the root groups for the two seed hydrogens

	if(core.Find_Root_Of_A_Hydrogen(seed_1,root_1)==FALSE) 
		{
		 // puts("Error: cannot figure out the seed H on the core");
		 return FALSE; 
		}
	if(frag.Find_Root_Of_A_Hydrogen(seed_2,root_2)==FALSE) 
		{
		 // puts("Error: cannot figure out the seed H on the fragment");
		 return FALSE; 
		}

	id1=root_1.neib[0].id; id2=root_1.center.id; 
	id3=root_2.center.id; id4=root_2.neib[0].id; 
	id5=seed_1; id6=seed_2;

/*
 		id1                         id4
	  	 \                          /
  	 core  	 id2 - id5 (H)  (H) id6 - id3    frag

        please notice: for improper torsion, id1 may equal to id5 
	and id4 may equal to id6
*/

	// now check the chemical viability of the new bond

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

	if(!strcmp(new_bond.type,"no")) 
		{
		 return FALSE;	// Invalid new bond 
		}
	else
		{
		 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);

		 // notice that here atom type may have been changed by 
		 // New_Bond_Viability_Check() above

		 core.atom[id2-1]=root_1.center; 
		 frag.atom[id3-1]=root_2.center; 
		}

	// now place the fragment on the anchor point

	tmp1=new_bond.length;
        tmp2=Distance(core.atom[id5-1].coor,core.atom[id2-1].coor);

	for(i=0;i<3;i++)
                {
                 anchor[i]=core.atom[id2-1].coor[i];
                 tmp3=core.atom[id5-1].coor[i]-core.atom[id2-1].coor[i];
                 anchor[i]+=(tmp3*tmp1/tmp2);
                 move[i]=anchor[i]-frag.atom[id3-1].coor[i];
                }

	frag.Translate(move);

	// now rotate the fragment to overlap the two seed bonds

	for(i=0;i<3;i++)
                {
                 v1[i]=frag.atom[id6-1].coor[i]-frag.atom[id3-1].coor[i];
                 v2[i]=core.atom[id2-1].coor[i]-frag.atom[id3-1].coor[i];
                }

        angle=Angle_Of_Two_Vectors(v1,v2);
        Cross_Multiply(v1,v2,axis);
        frag.Rotate(angle,axis,anchor);

	// now start to build the new molecule

	Ligand tmp_frag[24];
        Ligand tmp_result(core.num_atom+frag.num_atom+100,
                          core.num_bond+frag.num_bond+100);

	if(root_1.num_neib<=1||root_2.num_neib<=1) goto Improper_Torsion;
	else goto Proper_Torsion; 

Proper_Torsion:

	float e,tors_e,vdw_e;
	struct 
		{
		 float angle;
		 float e;
		 int valid;
		 int grow;
		 int link;
		} profile[24];

	// now rotate the fragment along the new bond to overlap
        // atom id1, id2, id3, id4 on the same plane

        angle=Torsion_Angle(core.atom[id1-1].coor,core.atom[id2-1].coor,
			    frag.atom[id3-1].coor,frag.atom[id4-1].coor);

        for(i=0;i<3;i++) 
		axis[i]=frag.atom[id3-1].coor[i]-core.atom[id2-1].coor[i];

        frag.Rotate(-angle,axis,anchor);

	// now rotate the fragment thoroughly along the new bond
	// record the energy profile and the ring_formation information 

	new_torsion.atom_1=core.atom[id1-1];
	new_torsion.atom_2=core.atom[id2-1];
	new_torsion.atom_3=frag.atom[id3-1];
	new_torsion.atom_4=frag.atom[id4-1];
	strcpy(new_torsion.type,new_bond.type);
	ff->Assign_Torsion_Parameters(new_torsion);

	for(i=0;i<3;i++)
		{
		 anchor[i]=frag.atom[id3-1].coor[i];
		 axis[i]=frag.atom[id3-1].coor[i]-core.atom[id2-1].coor[i];
		}

	for(k=0;k<24;k++)
		{
		 angle=k*15.00;

		 profile[k].angle=angle;
		 profile[k].e=0.000;
		 profile[k].valid=FALSE;
		 profile[k].grow=FALSE;
		 profile[k].link=FALSE;

		 tmp_frag[k]=frag;
		 tmp_frag[k].Rotate(angle,axis,anchor);

		 // calculate the VDW energy between ligand and fragment

		 e=tors_e=vdw_e=0.000;

		 for(i=0;i<core.num_atom;i++)
                        {
                         if(core.atom[i].id==id5) continue;
                         else if(core.atom[i].id==id2) continue;

			// check if atom i is the neighbor of atom id2 

			 mark1=FALSE;

			 for(l=0;l<root_1.num_neib;l++) 
			{
			 if(core.atom[i].id==root_1.neib[l].id) 
				{mark1=TRUE;break;}
			 else continue;
			}

                         for(j=0;j<tmp_frag[k].num_atom;j++)
                        {
                         if(tmp_frag[k].atom[j].id==id6) continue;
                         else if(tmp_frag[k].atom[j].id==id3) continue;

			// check if atom j is the neighbor of atom id3 

			 mark2=FALSE;

			 for(l=0;l<root_2.num_neib;l++)
				{
				 if(tmp_frag[k].atom[j].id==root_2.neib[l].id)
					{mark2=TRUE;break;}
				 else continue;
				}

			 if(mark1==TRUE&&mark2==TRUE) // 1-4 interaction
				{
			 	 vdw_e+=(ff->Cal_VDW_Energy(core.atom[i],tmp_frag[k].atom[j],1));
				}				
			 else
				{
                         	 vdw_e+=(ff->Cal_VDW_Energy(core.atom[i],tmp_frag[k].atom[j]));
				}
			}
                        }

		// calculate the torsion energy of the new bond

		 new_torsion.angle=(int)angle;
                 new_torsion.e=ff->Cal_Torsion_Energy(new_torsion);
                 tors_e=new_torsion.e;
                 e=tors_e+vdw_e; profile[k].e=e;
		}

	// find the energy minima and the possible candidate for linking

	if((profile[0].e<profile[23].e)&&
	   (profile[0].e<profile[1].e)&&
	   (profile[0].e<MINIMA_ENERGY_CUTOFF)) profile[0].grow=TRUE;

	for(i=1;i<=22;i++)
	{
	 if((profile[i].e<profile[i-1].e)&&
	    (profile[i].e<profile[i+1].e)&&
	    (profile[i].e<MINIMA_ENERGY_CUTOFF)) profile[i].grow=TRUE;
	 else continue; 
	}

	if((profile[23].e<profile[22].e)&&
	   (profile[23].e<profile[0].e)&&
	   (profile[23].e<MINIMA_ENERGY_CUTOFF)) profile[23].grow=TRUE;

	if(drand48()<parm->link_ratio)
		{
		 for(i=0;i<24;i++)
			{
		 	 if(profile[i].e<LINK_ENERGY_CUTOFF) continue;
		 	 else profile[i].link=TRUE;
			}
		}

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

	// eliminate the duplicates due to symmetric fragments

	for(i=0;i<23;i++)
                {
		 if(profile[i].valid==0) continue;
                 for(j=i+1;j<24;j++)
                        {
			 if(profile[j].valid==0) continue;
                         else if(Conformation_Duplicate_Check
                                 (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;