ga_link.c

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

# include "link.h"

Population::Population()
{
	extern Parameter *parm;

	num_member=num_record=0;

	member=new Ligand[parm->max_population];
	if(member==NULL) Memory_Allocation_Error();
}

Population::~Population()
{
	delete [] member;
}

void Population::Read_From_Lig(char *filename)
{
	extern Parameter *parm;
	FILE *fp;
	int i,j,mark;
	char buf[160],head[80];

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

	num_member=0;

	for(i=0;i<parm->max_population;i++)
		{
		 mark=FALSE;	// indicator of the end of the file

		 do 
			{
			 if(fgets(buf,160,fp)==NULL) {mark=TRUE;break;} 
			 else sscanf(buf,"%s",head);
			}
        	 while(strcmp(head,"<MOLECULE>"));

		 if(mark==TRUE) break;

		 sscanf(buf,"%*s%s",member[i].name);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<FORMULA>"));

		 sscanf(buf,"%*s%s",member[i].formula);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<WEIGHT>"));

		 sscanf(buf,"%*s%d",&member[i].weight);
			
		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<LOGP>"));

		 sscanf(buf,"%*s%f",&member[i].logp);
		 	
		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<BINDING_SCORE>"));

		 sscanf(buf,"%*s%f",&member[i].binding_score);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<CHEMICAL_SCORE>"));

                 sscanf(buf,"%*s%f",&member[i].chemical_score);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<ATOM>"));

		 sscanf(buf,"%*s%d",&member[i].num_atom);

		 delete [] member[i].atom;

		 member[i].atom=new Atom[member[i].num_atom];

		 for(j=0;j<member[i].num_atom;j++)
		{
		 fgets(buf,160,fp);
		 sscanf(buf,"%d%s%f%f%f%s%s%*s%d%f", 
			&member[i].atom[j].id, 
			 member[i].atom[j].name,
			&member[i].atom[j].coor[0], 
			&member[i].atom[j].coor[1], 
			&member[i].atom[j].coor[2],
			 member[i].atom[j].type, 
			 member[i].atom[j].xtype, 
			&member[i].atom[j].valid, 
			&member[i].atom[j].score);
		}

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<BOND>"));

		 sscanf(buf,"%*s%d",&member[i].num_bond);

		 delete [] member[i].bond;

		 member[i].bond=new Bond[member[i].num_bond];

		 for(j=0;j<member[i].num_bond;j++)
		{
		 fgets(buf,160,fp);
		 sscanf(buf,"%d%d%d%s%*s%d", 
			&member[i].bond[j].id, 
			&member[i].bond[j].atom_1,
			&member[i].bond[j].atom_2,
			 member[i].bond[j].type,  
			&member[i].bond[j].valid);
		}

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<INDEX>"));

		 sscanf(buf,"%*s%s",member[i].index);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<END>"));

		 member[i].valid=1; num_member++;
		}

	fclose(fp);

	if(num_member==0) Lig_Format_Error(filename);

	printf("%d molecules are read in.\n", num_member);

	for(i=0;i<num_member;i++)
                {
                 member[i].Detect_Connections();
                 member[i].Assign_Atom_Parameters();
                }

	return;
}

void Population::Write_Out_Lig(char *filename) const
{
	FILE *fp;	
	int i,j;

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

	for(i=0;i<num_member;i++)
		{
		 fprintf(fp,"<MOLECULE> %s\n", member[i].name);
		 fprintf(fp,"<FORMULA> %s\n", member[i].formula);
		 fprintf(fp,"<WEIGHT> %d\n", member[i].weight);
		 fprintf(fp,"<LOGP> %6.2f\n", member[i].logp);
		 fprintf(fp,"<BINDING_SCORE> %5.2f\n", member[i].binding_score);
		 fprintf(fp,"<CHEMICAL_SCORE> %6.1f\n", 
				member[i].chemical_score);

		 fprintf(fp,"<ATOM> %d\n", member[i].num_atom);

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

		 fprintf(fp,"<BOND> %d\n", member[i].num_bond);

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

		 fprintf(fp,"<INDEX> %s\n", member[i].index);

		 fprintf(fp,"<END>\n");
		}

	fclose(fp);

	return;
}

void Population::Count_Wins()
{
	int i,j,wins;

	if(num_member==0) return;

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

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

	 if(wins>0) member[i].wins=wins;
	 else member[i].wins=1;
	}

	return;
}

void Population::Statistics()
{
	int i;

	if(num_member==0) {printf("\nWarning: empty population.\n"); return;}

	max_weight=0;
        min_weight=1000;
        ave_weight=0;

        max_binding_score=-100.00;
        min_binding_score= 100.00;
        ave_binding_score=   0.00;

        max_chemical_score=-1000.00;
        min_chemical_score= 1000.00;
        ave_chemical_score=    0.00;

        max_logp=-10.00;
        min_logp= 10.00;
        ave_logp=  0.00;

	for(i=0;i<num_member;i++)
                {
                 if(member[i].weight>max_weight) max_weight=member[i].weight;
                 if(member[i].weight<min_weight) min_weight=member[i].weight;
                 ave_weight+=(member[i].weight);

                 if(member[i].binding_score>max_binding_score) max_binding_score=member[i].binding_score;
                 if(member[i].binding_score<min_binding_score) min_binding_score=member[i].binding_score;
                 ave_binding_score+=(member[i].binding_score);

                 if(member[i].chemical_score>max_chemical_score) max_chemical_score=member[i].chemical_score;
                 if(member[i].chemical_score<min_chemical_score) min_chemical_score=member[i].chemical_score;
                 ave_chemical_score+=(member[i].chemical_score);

                 if(member[i].logp>max_logp) max_logp=member[i].logp;
                 if(member[i].logp<min_logp) min_logp=member[i].logp;
                 ave_logp+=(member[i].logp);
                }

	ave_weight/=(num_member);
        ave_binding_score/=(num_member);
        ave_chemical_score/=(num_member);
        ave_logp/=(num_member);

	printf("\n****************************************************************\n");

	printf("Current population: %d\n", num_member);

	printf("Molecular weight range: %d --- %d,  average = %d;\n",  
                min_weight, max_weight, ave_weight);
        printf("LogP range: %6.2f --- %6.2f,  average = %6.2f;\n",
                min_logp, max_logp, ave_logp);
        printf("Binding score range: %6.2f --- %6.2f,  average = %6.2f;\n",
                min_binding_score, max_binding_score, ave_binding_score);
        // printf("Chemical score range: %7.2f --- %7.2f,  average = %7.2f;\n",
        //        min_chemical_score, max_chemical_score, ave_chemical_score);

	printf("****************************************************************\n\n");

	return;
}

void Population::Initialize(Ligand core, int num_generation)
{
	extern Parameter *parm;
	extern FragLibrary *fraglib;
	int i,j,k,tmp,mark,count;
	int max_record,num_core_seed,num_frag_seed,num_frag,num_new;
	Ligand lig,frag,result[24];
	int *core_seed_list,*frag_seed_list,*frag_list,*tmp_seed_list;
	Growing_Record tmp_record;

	core_seed_list=new int[core.num_atom];
	if(core_seed_list==NULL) Memory_Allocation_Error();

	frag_seed_list=new int[500];
	if(frag_seed_list==NULL) Memory_Allocation_Error();

	tmp_seed_list=new int[500];
	if(tmp_seed_list==NULL) Memory_Allocation_Error();

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

	// first, get the seed hydrogens on the core structure

	num_core_seed=core.Make_Seed_List(core_seed_list);

	// growing on each seed hydrogen on the core structure

	max_record=parm->max_population*fraglib->num_frag;

        growing_record=new Growing_Record[max_record];
        if(growing_record==NULL) Memory_Allocation_Error();

	num_member=num_record=0;

	lig=core;

	while(1)
	{
	 num_new=0;

	 // first, estimate the number of fragments which will be
	 // grown on each growing site and make the fragment list.

	 num_frag=(parm->max_population-num_member)/num_core_seed;

	 for(i=0;i<num_core_seed;i++)
		{
	 	 num_frag=fraglib->Make_Fragment_List(num_frag,frag_list);