ga_link.c

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

		 // applying each fragment in the list

		 for(j=0;j<num_frag;j++)
			{
			 mark=fraglib->Get_A_Fragment(frag_list[j],frag);
			 if(mark==FALSE) continue; // get fragment failed 

			 num_frag_seed=frag.Make_Seed_List(frag_seed_list);

			 for(k=0;k<num_frag_seed;k++)
			{
			 tmp_record.ligand_id=1;
			 tmp_record.lig_seed=core_seed_list[i];
			 tmp_record.fragment_id=frag_list[j];
			 tmp_record.frag_seed=frag_seed_list[k];

			 mark=Growing_Record_Check(tmp_record);
			 if(mark==FALSE) continue;
			 else frag_seed_list[k]=0;  // has been tried before 
			}

			 tmp=0;		 
			 for(k=0;k<num_frag_seed;k++)
			{
			 if(frag_seed_list[k]==0) continue;
			 else {tmp_seed_list[tmp]=frag_seed_list[k]; tmp++;}
			}

			 if(tmp==0) 
				{
				 // puts("no available seed H on the fragment");
				 continue;
				}
			 else k=(int)(drand48()*tmp);

                         tmp_record.ligand_id=1;
                         tmp_record.lig_seed=core_seed_list[i];
                         tmp_record.fragment_id=frag_list[j];
                         tmp_record.frag_seed=tmp_seed_list[k];

			 growing_record[num_record]=tmp_record;
			 if(num_record<(max_record-1)) num_record++;

			 mark=lig.Generate_New_Structure
			      (core_seed_list[i],frag,tmp_seed_list[k],result);

			 if(mark==FALSE) 
				{
				 // puts("growing fails");
				 continue;
				}
			 else if(mark>(parm->max_population-num_member))
				goto End_of_growing; // too many! stop growing

			 // add the new molecules

			 count=0;	

			 for(k=0;k<24;k++)
			{
			 if(result[k].valid==0) continue;
			 else if(Duplicate_Check(result[k])==TRUE) continue;
			 else
				{
			 	 member[num_member]=result[k];
			 	 sprintf(member[num_member].name,"No_%d_%d", 
				 	 num_generation, num_member+1);
			 	 member[num_member].id=num_member+1;
			 	 num_member++; count++;
				}
			}

			 num_new+=count;
			}
		}
	 if(num_new<num_core_seed) break; 
	}

	End_of_growing:

	delete [] core_seed_list;
	delete [] frag_seed_list;
	delete [] tmp_seed_list;
	delete [] frag_list;
	delete [] growing_record;
	
	return;
}

int Population::Growing_Record_Check(Growing_Record &record) const
{
	int i;

	for(i=0;i<num_record;i++)
	{
 	 if(record.ligand_id!=growing_record[i].ligand_id) continue; 
	 else if(record.lig_seed!=growing_record[i].lig_seed) continue; 
	 else if(record.fragment_id!=growing_record[i].fragment_id) continue;
	 else if(record.frag_seed!=growing_record[i].frag_seed) continue;
	 else return TRUE; 
	}

	return FALSE;
}

int Population::Select_A_Mother_Molecule() const
{
	int i,total,tmp,mark;
	struct Wheel
		{
		 int min;
		 int max;
		} *wheel;

	if(num_member==0) return FALSE;

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

	// initialize the wheel according to wins 

	total=0;

	for(i=0;i<num_member;i++)
		{
		 wheel[i].min=total;
		 total+=member[i].wins;
		 wheel[i].max=total;
		}

	// select the mother molecule

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

	 mark=0;

	 for(i=0;i<num_member;i++)
		{
		 if(tmp<wheel[i].min) continue;
		 else if(tmp>=wheel[i].max) continue;
		 else {mark=i+1;break;}
		}

	delete [] wheel;

	return mark;
}

int Population::Growing_On_Molecule(int mother_id, Ligand result[24])
{
	extern Parameter *parm;
	extern FragLibrary *fraglib;
	Ligand lig,frag;
	int num;
	int core_seed,frag_seed;
	int *frag_list; 

	num=0;

	frag_list=new int[fraglib->num_frag];
	if(frag_list==NULL) Memory_Allocation_Error();

	lig=member[mother_id-1];

	fraglib->Make_Fragment_List(1,frag_list);
	if(fraglib->Get_A_Fragment(frag_list[0],frag)==FALSE) goto End; 

	core_seed=lig.Select_A_Seed_Hydrogen();
	frag_seed=frag.Select_A_Seed_Hydrogen();

	if(core_seed==0||frag_seed==0) goto End;  // no available seed H 

	num=lig.Generate_New_Structure(core_seed,frag,frag_seed,result);

	End:

	delete [] frag_list;

	return num;
}

void Population::Add_A_Molecule_To_Population(Ligand &lig)
{
	extern Parameter *parm;
	int i,mark,min_wins;

	//if the population is not full yet, add the new molecule directly

	if(num_member<parm->max_population) 
	{
	 // duplicate check first

	 if(Duplicate_Check(lig)==TRUE) return;
	 else
		{
	 	 member[num_member]=lig; 
	 	 sprintf(member[num_member].name,"No_%d", num_member+1);
	 	 num_member++;
	 	 Count_Wins();
	 	 return;
		}
	}

	// otherwise, count the wins of lig

	lig.wins=0;

	for(i=0;i<num_member;i++)
	{
	 if(lig.binding_score<member[i].binding_score) lig.wins--;
	 else if(lig.binding_score>member[i].binding_score) lig.wins++;
	 else continue;
	}

	for(i=0;i<num_member;i++)
	{
	 if(lig.chemical_score<member[i].chemical_score) lig.wins--;
	 else if(lig.chemical_score>member[i].chemical_score) lig.wins++;
	 else continue;
	}

	if(lig.wins<=0) lig.wins=0;

	// find the worst one in the population

	mark=0; min_wins=lig.wins;

	for(i=0;i<num_member;i++)
		{
		 if(member[i].wins>=min_wins) continue;
		 else
			{
			 mark=i+1;
			 min_wins=member[i].wins;
			}
		}

	// finally, replace the worst one by lig

	if(mark>0) 
	{
	 // duplicate check first

	 if(Duplicate_Check(lig)==TRUE) return;
	 else
		{
	 	 strcpy(lig.name,member[mark-1].name);
	 	 member[mark-1]=lig; 
	 	 Count_Wins();
		}
	}

	return;
}

int Population::Add_A_Molecule_To_Result(Ligand &lig) const
{
	extern Parameter *parm;
	FILE *fp;
	int i;
	char filename[160];

	strcpy(filename,parm->ligands_file);

	// check whether all the fragments have been linked together

	if(lig.num_part!=1) return FALSE;

	// chemical viability check

	if(lig.Chemical_Viability_Check()==FALSE) return FALSE;

	if((fp=fopen(filename,"a"))==NULL) Openning_File_Error(filename);

	fprintf(fp,"<MOLECULE> %s\n", lig.name);
	fprintf(fp,"<FORMULA> %s\n", lig.formula);
	fprintf(fp,"<WEIGHT> %d\n", lig.weight);
	fprintf(fp,"<LOGP> %6.2f\n", lig.logp);
	fprintf(fp,"<BINDING_SCORE> %5.2f\n", lig.binding_score);
	fprintf(fp,"<CHEMICAL_SCORE> %6.1f\n", lig.chemical_score);

	fprintf(fp,"<ATOM> %d\n", lig.num_atom);

	for(i=0;i<lig.num_atom;i++)
                {
                 fprintf(fp,"%-3d  %-6s  %8.3f  %8.3f  %8.3f  %5s  ",
                         lig.atom[i].id,
                         lig.atom[i].name,
                         lig.atom[i].coor[0],
                         lig.atom[i].coor[1],
                         lig.atom[i].coor[2],
                         lig.atom[i].type);
                 fprintf(fp,"%-20s <%1d> %-2d  %6.2f\n",
                         lig.atom[i].xtype,
                         lig.atom[i].part,
                         lig.atom[i].valid,
                         lig.atom[i].score);
                }

	fprintf(fp,"<BOND> %d\n", lig.num_bond);

	for(i=0;i<lig.num_bond;i++)
                {
                 fprintf(fp,"%-3d  %-3d  %-3d  %3s <%1d> %1d\n",
                         lig.bond[i].id,
                         lig.bond[i].atom_1,
                         lig.bond[i].atom_2,
                         lig.bond[i].type,
                         lig.bond[i].part,
                         lig.bond[i].valid);
                }

        fprintf(fp,"<INDEX> %s\n", lig.index);
        fprintf(fp,"<END>\n");

	fclose(fp);

	return TRUE;
}

int Population::Duplicate_Check(Ligand &lig) const
{
        int i,mark;

        mark=FALSE;

        for(i=0;i<num_member;i++)
                {
		 if(member[i].valid==0) continue;
                 else if(Molecule_Duplicate_Check(lig,member[i])==TRUE)
                        {
                         mark=TRUE; break;
                        }
                 else continue;
                }

        return mark;
}