ligand.c

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

	delete [] atom; 

	atom=new Atom[num_atom];
	if(atom==NULL) Memory_Allocation_Error();

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

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

	sscanf(buf,"%*s%d",&num_bond);

	delete [] bond; 

	bond=new Bond[num_bond];
	if(bond==NULL) Memory_Allocation_Error();

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

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

	sscanf(buf,"%*s%s",index);

	fclose(fp);

	valid=1;	// read the molecule okay

	return;
}

void Ligand::Write_Out_Lig(char *filename) const
{
	FILE *fp;
	int i;

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

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

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

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

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

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

	fclose(fp);

	return;
}

void Ligand::Assign_Atom_Parameters()
{
	extern ForceField *ff;
	int i;

	Detect_Rings();

	Assign_Atom_Xtype();

	for(i=0;i<num_atom;i++) ff->Assign_Atom_Parameters(atom[i]);

	for(i=0;i<num_atom;i++) ff->Assign_Atom_Xparameters(atom[i]);

	return;
}

void Ligand::Get_Molecular_Properties()
{
	this->weight=Get_Weight();
	
	this->Get_Formula();

	this->Get_Structure_Index();

	this->logp=Get_LogP();

	this->num_donor=Get_Num_Donor_Atom();

	this->num_acceptor=Get_Num_Acceptor_Atom();
	
	return;
}

int Ligand::Get_Weight() const
{
        int i,tmp_weight;

        tmp_weight=0;

        for(i=0;i<num_atom;i++) tmp_weight+=atom[i].weight;

        return tmp_weight;
}

void Ligand::Get_Formula()
{
	int i;
        int num[11]={0,0,0,0,0,0,0,0,0,0,0}; // C,H,N,O,P,S,F,Cl,Br,I,Un
	char element[11][3];
	char tmp[10],tmp_formula[40];

	strcpy(element[0],"C");
	strcpy(element[1],"H");
	strcpy(element[2],"N");
	strcpy(element[3],"O");
	strcpy(element[4],"P");
	strcpy(element[5],"S");
	strcpy(element[6],"F");
	strcpy(element[7],"Cl");
	strcpy(element[8],"Br");
	strcpy(element[9],"I");
	strcpy(element[10],"Un");

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) num[10]++;
		 else if(!strcmp(atom[i].type,"C.3")) num[0]++;
		 else if(!strcmp(atom[i].type,"C.2")) num[0]++;
		 else if(!strcmp(atom[i].type,"C.1")) num[0]++;
		 else if(!strcmp(atom[i].type,"C.cat")) num[0]++;
                 else if(!strcmp(atom[i].type,"C.ar")) num[0]++;
                 else if(!strcmp(atom[i].type,"H")) num[1]++;
                 else if(!strcmp(atom[i].type,"H.spc")) num[1]++;
		 else if(!strcmp(atom[i].type,"N.4")) num[2]++;
                 else if(!strcmp(atom[i].type,"N.3")) num[2]++;
                 else if(!strcmp(atom[i].type,"N.2")) num[2]++;
		 else if(!strcmp(atom[i].type,"N.1")) num[2]++;
                 else if(!strcmp(atom[i].type,"N.ar")) num[2]++;
                 else if(!strcmp(atom[i].type,"N.pl3")) num[2]++;
		 else if(!strcmp(atom[i].type,"N.am")) num[2]++;
                 else if(!strcmp(atom[i].type,"O.3")) num[3]++;
		 else if(!strcmp(atom[i].type,"O.2")) num[3]++;
                 else if(!strcmp(atom[i].type,"O.co2")) num[3]++;
                 else if(!strcmp(atom[i].type,"P.3")) num[4]++;
		 else if(!strcmp(atom[i].type,"S.3")) num[5]++;
                 else if(!strcmp(atom[i].type,"S.2")) num[5]++;
                 else if(!strcmp(atom[i].type,"S.o")) num[5]++;
		 else if(!strcmp(atom[i].type,"S.o2")) num[5]++;
                 else if(!strcmp(atom[i].type,"F")) num[6]++;
                 else if(!strcmp(atom[i].type,"Cl")) num[7]++;
                 else if(!strcmp(atom[i].type,"Br")) num[8]++;
		 else if(!strcmp(atom[i].type,"I")) num[9]++;
                 else num[10]++;	// Unknown
		}

	strcpy(tmp_formula,"");

	for(i=0;i<11;i++)
		{
		 if(num[i]==0) continue;
		 else
			{
			 strcat(tmp_formula,element[i]);
			 if(num[i]==1) continue;
			 else 
				{
				 sprintf(tmp,"%d", num[i]);
				 strcat(tmp_formula,tmp);
				}
			}
		}

	strcpy(this->formula,tmp_formula);

        return;
}

void Ligand::Get_Structure_Index()
{
	// The structure index is a 60-bit string, recording the number
	// of 25 types of atoms and 5 types of bonds in hexdecimal number. 
	// Each atom type or bond type is assigned two bits.
	// The string will look like: 1A03000E......

	int i,tmp[30],high_tmp,low_tmp;

	strcpy(this->index,"");

	for(i=0;i<30;i++) tmp[i]=0;

        for(i=0;i<num_atom;i++)
                {
                 if(atom[i].valid==0) continue;
                 else if(!strcmp(atom[i].type,"C.3")) tmp[0]++;
		 else if(!strcmp(atom[i].type,"C.2")) tmp[1]++;
		 else if(!strcmp(atom[i].type,"C.cat")) tmp[2]++;
		 else if(!strcmp(atom[i].type,"C.ar")) tmp[3]++;
		 else if(!strcmp(atom[i].type,"C.1")) tmp[4]++;
                 else if(!strcmp(atom[i].type,"O.3"))  tmp[5]++;
                 else if(!strcmp(atom[i].type,"O.2")) tmp[6]++;
                 else if(!strcmp(atom[i].type,"O.co2")) tmp[7]++;
                 else if(!strcmp(atom[i].type,"N.4")) tmp[8]++;
                 else if(!strcmp(atom[i].type,"N.3")) tmp[9]++;
                 else if(!strcmp(atom[i].type,"N.am")) tmp[10]++;
                 else if(!strcmp(atom[i].type,"N.pl3")) tmp[11]++;
                 else if(!strcmp(atom[i].type,"N.2")) tmp[12]++;
                 else if(!strcmp(atom[i].type,"N.ar")) tmp[13]++;
                 else if(!strcmp(atom[i].type,"N.1")) tmp[14]++;
                 else if(!strcmp(atom[i].type,"S.3")) tmp[15]++;
                 else if(!strcmp(atom[i].type,"S.2")) tmp[16]++;
                 else if(!strcmp(atom[i].type,"S.o")) tmp[17]++;
                 else if(!strcmp(atom[i].type,"S.o2")) tmp[18]++;
                 else if(!strcmp(atom[i].type,"P.3")) tmp[19]++;
                 else if(!strcmp(atom[i].type,"F")) tmp[20]++;
                 else if(!strcmp(atom[i].type,"Cl")) tmp[21]++;
                 else if(!strcmp(atom[i].type,"Br")) tmp[22]++;
                 else if(!strcmp(atom[i].type,"I")) tmp[23]++;
                 else tmp[24]++;        // H plus Unknown
                }

	for(i=0;i<num_bond;i++)
		{
		 if(bond[i].valid==0) continue;
		 else if(atom[bond[i].atom_1-1].type[0]=='H') continue;
		 else if(atom[bond[i].atom_2-1].type[0]=='H') continue;
		 else if(!strcmp(bond[i].type,"1")) tmp[25]++;
		 else if(!strcmp(bond[i].type,"2")) tmp[26]++;
		 else if(!strcmp(bond[i].type,"3")) tmp[27]++;
		 else if(!strcmp(bond[i].type,"ar")) tmp[28]++;
		 else if(!strcmp(bond[i].type,"am")) tmp[29]++;
		 else continue;
		}

	for(i=0;i<30;i++)
		{
		 high_tmp=tmp[i]/16; low_tmp=tmp[i]%16;

		 if(high_tmp==0) this->index[i*2]='0';
		 else if(high_tmp==1) this->index[i*2]='1';
		 else if(high_tmp==2) this->index[i*2]='2';
		 else if(high_tmp==3) this->index[i*2]='3';
                 else if(high_tmp==4) this->index[i*2]='4';
		 else if(high_tmp==5) this->index[i*2]='5';
                 else if(high_tmp==6) this->index[i*2]='6';
		 else if(high_tmp==7) this->index[i*2]='7';
                 else if(high_tmp==8) this->index[i*2]='8';
		 else if(high_tmp==9) this->index[i*2]='9';
                 else if(high_tmp==10) this->index[i*2]='A';
		 else if(high_tmp==11) this->index[i*2]='B';
                 else if(high_tmp==12) this->index[i*2]='C';
		 else if(high_tmp==13) this->index[i*2]='D';
                 else if(high_tmp==14) this->index[i*2]='E';
		 else this->index[i*2]='F';

		 if(low_tmp==0) this->index[i*2+1]='0';
                 else if(low_tmp==1) this->index[i*2+1]='1';
                 else if(low_tmp==2) this->index[i*2+1]='2';
                 else if(low_tmp==3) this->index[i*2+1]='3';
                 else if(low_tmp==4) this->index[i*2+1]='4';
                 else if(low_tmp==5) this->index[i*2+1]='5';
                 else if(low_tmp==6) this->index[i*2+1]='6';
                 else if(low_tmp==7) this->index[i*2+1]='7';
                 else if(low_tmp==8) this->index[i*2+1]='8';
                 else if(low_tmp==9) this->index[i*2+1]='9';
                 else if(low_tmp==10) this->index[i*2+1]='A';
                 else if(low_tmp==11) this->index[i*2+1]='B';
                 else if(low_tmp==12) this->index[i*2+1]='C';
                 else if(low_tmp==13) this->index[i*2+1]='D';
                 else if(low_tmp==14) this->index[i*2+1]='E';
                 else this->index[i*2+1]='F';
		}

	this->index[60]='\0';

	return;
}

int Ligand::Get_Num_Donor_Atom() const
{
	int i,num=0;

        for(i=0;i<num_atom;i++)
                {
                 if(atom[i].valid==0) continue;
                 else if(!strcmp(atom[i].hb,"D")) num++;
		 else if(!strcmp(atom[i].hb,"DA")) num++;
                 else continue;
                }

	return num;
}

int Ligand::Get_Num_Acceptor_Atom() const
{
	int i,num=0;

        for(i=0;i<num_atom;i++)
                {
                 if(atom[i].valid==0) continue;
                 else if(!strcmp(atom[i].hb,"A")) num++;
		 else if(!strcmp(atom[i].hb,"DA")) num++;
                 else continue;
                }

	return num;
}

char *Ligand::Get_Atom_Hybridizing_Type(Atom atm) const
{
         if(atm.valid==0) return "none";
         else if(!strcmp(atm.type,"C.3"))   return "sp3";
         else if(!strcmp(atm.type,"C.2"))   return "sp2";
         else if(!strcmp(atm.type,"C.1"))   return "sp";
         else if(!strcmp(atm.type,"C.cat")) return "sp2";
         else if(!strcmp(atm.type,"C.ar"))  return "sp2";
         else if(!strcmp(atm.type,"H"))     return "none";
         else if(!strcmp(atm.type,"H.spc")) return "none";
         else if(!strcmp(atm.type,"N.4"))   return "sp3";
         else if(!strcmp(atm.type,"N.3"))   return "sp3";
         else if(!strcmp(atm.type,"N.2"))   return "sp2";
         else if(!strcmp(atm.type,"N.1"))   return "sp";
         else if(!strcmp(atm.type,"N.ar"))  return "sp2";
         else if(!strcmp(atm.type,"N.pl3")) return "sp2";
         else if(!strcmp(atm.type,"N.am"))  return "sp2";
         else if(!strcmp(atm.type,"O.3"))   return "sp3";
         else if(!strcmp(atm.type,"O.2"))   return "sp2";
         else if(!strcmp(atm.type,"O.co2")) return "sp2";
         else if(!strcmp(atm.type,"P.3"))   return "sp3";
         else if(!strcmp(atm.type,"S.3"))   return "sp3";
         else if(!strcmp(atm.type,"S.2"))   return "sp2";
         else if(!strcmp(atm.type,"S.o"))   return "sp3";
         else if(!strcmp(atm.type,"S.o2"))  return "sp3";
         else if(!strcmp(atm.type,"F"))     return "none";
         else if(!strcmp(atm.type,"Cl"))    return "none";
         else if(!strcmp(atm.type,"Br"))    return "none";
         else if(!strcmp(atm.type,"I"))     return "none";
         else return "none";
}

void Ligand::Detect_Connections()
{
	int i,j,k,mark,tmp1,tmp2,old_part,new_part,id1,id2;

	// first, detect how many fragments are in the molecule

	// initialize the molecule as the assembly of isolated atoms

	for(i=0;i<num_atom;i++) atom[i].part=atom[i].id;

	// second, merge the atoms into fragments 

	for(i=0;i<num_bond;i++)
		{
		 tmp1=bond[i].atom_1; tmp2=bond[i].atom_2;
		 if(atom[tmp1-1].part>=atom[tmp2-1].part) 
			{
			 old_part=atom[tmp1-1].part; 
			 new_part=atom[tmp2-1].part;
			}
		 else  
			{
			 old_part=atom[tmp2-1].part; 
			 new_part=atom[tmp1-1].part;
			}

		 for(j=0;j<num_atom;j++)
			{
			 if(atom[j].part!=old_part) continue;
			 else atom[j].part=new_part;
			} 
		}

	// third, count the number of fragments 
	// also re-arrange the ids to make them continuous

	int *part_list;

	part_list=new int[num_atom];
	if(part_list==NULL) Memory_Allocation_Error();

	for(i=0;i<num_atom;i++) part_list[i]=0;

	num_part=0;

	for(i=0;i<num_atom;i++)
		{
		 mark=FALSE;

		 for(j=0;j<num_part;j++)
			{
			 if(atom[i].part!=part_list[j]) continue;
			 else {mark=TRUE; break;}
			}

		 if(mark==TRUE) {atom[i].part=j+1;}	// not a new component
		 else					// a new component
			{
			 part_list[num_part]=atom[i].part;
			 num_part++;
			 atom[i].part=num_part;
			}
		}

	delete [] part_list;

	// also define all the bonds

	for(i=0;i<num_bond;i++) bond[i].part=atom[bond[i].atom_1-1].part;