ga_grow.c
来自「药物开发中的基于结构的从头设计代码」· C语言 代码 · 共 879 行 · 第 1/2 页
C
879 行
# include "grow.h"
Population::Population()
{
extern Parameter *parm;
num_member=num_parent=num_record=0;
member=new Ligand[parm->max_population];
if(member==NULL) Memory_Allocation_Error();
parent=new Ligand[parm->num_parent];
if(parent==NULL) Memory_Allocation_Error();
}
Population::~Population()
{
delete [] member;
delete [] parent;
}
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::Record_Results(char *filename) const
{
FILE *fp;
int i,j;
if(num_member==0) return;
if((fp=fopen(filename,"a"))==NULL) Openning_File_Error(filename);
for(i=0;i<num_member;i++)
{
if(member[i].valid==0) continue;
else if(member[i].Chemical_Viability_Check()==FALSE) continue;
else
{
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::Select_Elites(float ratio,int num_generation)
{
int i,j,num,count,mark,tmp;
Ligand *pool;
struct Index
{
int id;
int wins;
};
Index *index,temp;
if(num_member==0) return;
// first, rank all the members according to their wins and make index
index=NULL;
index=new Index[num_member];
if(index==NULL) Memory_Allocation_Error();
for(i=0;i<num_member;i++)
{
index[i].id=i+1;
index[i].wins=member[i].wins;
}
if(num_member>1)
{
for(i=0;i<num_member-1;i++)
for(j=i+1;j<num_member;j++)
{
if(index[i].wins>=index[j].wins) continue;
else
{
temp=index[i];
index[i]=index[j];
index[j]=temp;
}
}
}
// then, select the top members
num=(int)(num_member*ratio);
pool=new Ligand[num];
if(pool==NULL) Memory_Allocation_Error();
count=0;
for(i=0;i<num_member;i++)
{
if(count>=num) break; // pool is full
tmp=index[i].id-1;
if(member[tmp].valid==0) continue;
else
{
mark=FALSE;
for(j=0;j<count;j++)
{
if(Molecule_Duplicate_Check(pool[j],member[tmp])==TRUE)
{mark=TRUE; break;}
else continue;
}
if(mark==TRUE) continue; // duplicate found
else
{
pool[count]=member[tmp];
count++;
}
}
}
// rename these elites
for(i=0;i<count;i++) sprintf(pool[i].name,"No_%d_%d",num_generation,i+1);
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?