forcefield.c

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

# include "link.h"

ForceField::ForceField(char *dirname)
{
	char filename[160];

	num_atomtype=0;	atom=NULL;
	num_bondtype=0; bond=NULL;
	num_torstype=0; torsion=NULL;
	num_xatomtype=0; xatom=NULL;

	strcpy(filename,dirname); strcat(filename,"ATOM_DEF");
	Read_ATOM_DEF(filename);

	strcpy(filename,dirname); strcat(filename,"BOND_DEF");
	Read_BOND_DEF(filename);

	strcpy(filename,dirname); strcat(filename,"TORSION_DEF");
	Read_TORSION_DEF(filename);

	strcpy(filename,dirname); strcat(filename,"XATOM_DEF");
	Read_XATOM_DEF(filename);
}

ForceField::~ForceField()
{
	delete [] atom;
	delete [] bond;
	delete [] torsion;
	delete [] xatom;
}

void ForceField::Read_ATOM_DEF(char *filename)
{
	FILE *fp;
	int tmp;
	char buf[160];

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

	// first, determine num_atomtype

	tmp=0;

	for(;;)
		{
		 if(fgets(buf,160,fp)==NULL) break;
		 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else tmp++;
		}

	rewind(fp);

	num_atomtype=tmp;

	atom=new ATOM_DEF[num_atomtype];
	if(atom==NULL) Memory_Allocation_Error();

	// second, read atom parameters

	tmp=0;

	for(;;)
		{
		 if(fgets(buf,160,fp)==NULL) break;
		 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else
			{
			 sscanf(buf,"%*d%s%f%f%d", 
				 atom[tmp].type, 
				&atom[tmp].r,
				&atom[tmp].eps, 
				&atom[tmp].weight);
			 tmp++;
			}
		}

	fclose(fp);

	if(tmp!=num_atomtype) Reading_File_Error(filename);

	return;
}

void ForceField::Read_BOND_DEF(char *filename)
{
	FILE *fp;
	int tmp;
	char buf[160];

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

	// first, determine num_bondtype

        tmp=0;

	for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else tmp++;
		}

	rewind(fp);

	num_bondtype=tmp;

	bond=new BOND_DEF[num_bondtype];
	if(bond==NULL) Memory_Allocation_Error();

	// second, read bond parameters

	tmp=0;

        for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else 
			{
			 sscanf(buf,"%s%s%s%f", 
				 bond[tmp].atom_1, 
				 bond[tmp].atom_2,
				 bond[tmp].type, 
				&bond[tmp].length);
			 tmp++;
			}
		}

	fclose(fp);

	if(num_bondtype!=tmp) Reading_File_Error(filename);

	return;
}

void ForceField::Read_TORSION_DEF(char *filename)
{
	FILE *fp;
	int tmp,mark;
	char buf[160];

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

	// first, determine num_torstype
 
	tmp=0;

        for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
                 else tmp++;
                }

	rewind(fp);

	num_torstype=tmp;

	torsion=new TORS_DEF[num_torstype];
	if(torsion==NULL) Memory_Allocation_Error();

	// second, read torsion parameters

        tmp=0;

        for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
                 else
                        {
                         sscanf(buf,"%s%s%s%s%s%f%d",
                                 torsion[tmp].atom_1, 
				 torsion[tmp].atom_2,
				 torsion[tmp].atom_3, 
				 torsion[tmp].atom_4,
                                 torsion[tmp].type, 
				&torsion[tmp].V, 
				&mark);

			 if(mark>0) {torsion[tmp].n=mark; torsion[tmp].S=1;}
			 else {torsion[tmp].n=-mark; torsion[tmp].S=-1;}
                         tmp++;
                        }
                }

        fclose(fp);

        if(num_torstype!=tmp) Reading_File_Error(filename);

	return;
}

void ForceField::Read_XATOM_DEF(char *filename)
{
	FILE *fp;
	int tmp;
	char buf[160];

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

	// first, determine num_xatomtype

        tmp=0;

        for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else tmp++;
		}

	rewind(fp);

	num_xatomtype=tmp;

	xatom=new XATOM_DEF[num_xatomtype];
	if(xatom==NULL) Memory_Allocation_Error();

	// second, read xatom parameters

        tmp=0;

        for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else if(buf[0]=='#') continue;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else
			{
			 sscanf(buf,"%s%f%d%s%f", 
				 xatom[tmp].type, 
				&xatom[tmp].R,
				&xatom[tmp].weight, 
				 xatom[tmp].hb,
				&xatom[tmp].logp);
			 tmp++;
			}
		}

	fclose(fp);

	if(num_xatomtype!=tmp) Reading_File_Error(filename);

	return;
}

void ForceField::Show_Content() const
{
	int i;

        printf("Number of atom types is %d\n",num_atomtype);

        for(i=0;i<num_atomtype;i++)
                {
                 printf("%10s %5.3f %5.3f %3d\n",
                        atom[i].type,
                        atom[i].r,
                        atom[i].eps,
                        atom[i].weight);
                }

        printf("Number of bond types is %d\n", num_bondtype);

        for(i=0;i<num_bondtype;i++)
                {
                 printf("%10s %10s %5s %5.3f\n",
                        bond[i].atom_1,
                        bond[i].atom_2,
                        bond[i].type,
                        bond[i].length);
                }

        printf("Number of torsion angle types is %d\n", num_bondtype);

        for(i=0;i<num_torstype;i++)
                {
                 printf("%10s %10s %10s %10s %5s %5.3f %2d %2d\n",
                        torsion[i].atom_1,
                        torsion[i].atom_2,
                        torsion[i].atom_3,
                        torsion[i].atom_4,
                        torsion[i].type,
                        torsion[i].V,
                        torsion[i].n,
                        torsion[i].S);
                }


        printf("Number of X atom types is %d\n", num_xatomtype);

        for(i=0;i<num_xatomtype;i++)
                {
                 printf("%10s %5.3f %d %2s %6.3f\n",
                        xatom[i].type,
                        xatom[i].R,
                        xatom[i].weight,
                        xatom[i].hb,
                        xatom[i].logp);
                }

	return;
}

void ForceField::Assign_Atom_Parameters(Atom &atm) const
{
	int i,mark;

	mark=0;

	for(i=0;i<num_atomtype;i++)
		{
		 if(strcmp(atm.type,atom[i].type)) continue;
		 else
			{
			 atm.r=atom[i].r; 
			 atm.eps=atom[i].eps; 
			 atm.weight=atom[i].weight;
			 mark=1; break;
			}
		}

	if(mark==0)
		{
		 printf("Warning: atom %s has no parameter\n",atm.type);
		 puts("This atom will be overlooked.");
		 atm.valid=0;
		}

	return;
}

void ForceField::Assign_Atom_Xparameters(Atom &atm) const
{
	int i,mark;

	mark=0;

	for(i=0;i<num_xatomtype;i++)
		{
		 if(strcmp(atm.xtype,xatom[i].type)) continue;
		 else
			{
			 atm.R=xatom[i].R; 
			 strcpy(atm.hb,xatom[i].hb);
			 atm.logp=xatom[i].logp;
			 mark=1; break;
			}
		}

	if(mark==0)
		{
		 printf("Warning: atom %s has no parameter\n",atm.xtype);
		 puts("This atom will be overlooked.");
		 atm.valid=0;
		}

	return;
}

float ForceField::Get_Bond_Length(char *a1, char *a2, char *bond_type) const
{
        int i,mark;
	float length;
	char atom_1[20], atom_2[20];

	strcpy(atom_1,a1); strcpy(atom_2,a2);

        if(!strcmp(atom_1,"H.spc")) strcpy(atom_1,"H");
        if(!strcmp(atom_2,"H.spc")) strcpy(atom_2,"H");

        mark=0;

        for(i=0;i<num_bondtype;i++)
                {
                 if(strcmp(bond_type,bond[i].type)) continue;
                 else if(!strcmp(atom_1,bond[i].atom_1)&&!strcmp(atom_2,bond[i].atom_2))
                        {length=bond[i].length; mark=1; break;}
                 else if(!strcmp(atom_2,bond[i].atom_1)&&!strcmp(atom_1,bond[i].atom_2))
                        {length=bond[i].length; mark=1; break;}
                 else continue;
                }

        if(mark==0)
                {
                 printf("Warning: no suitable bond parameter for ");
                 printf("%s - %s - %s\n", atom_1, bond_type, atom_2);
		 printf("Default parameter is assigned.\n");
                 length=1.54;       // default value for C.3-C.3
                }

        return length;
}

void ForceField::Assign_Torsion_Parameters(Torsion &tors) const
{
	int i,mark,tmp;

	if(!strcmp(tors.atom_1.type,"H.spc")) strcpy(tors.atom_1.type,"H");
	if(!strcmp(tors.atom_2.type,"H.spc")) strcpy(tors.atom_2.type,"H");
	if(!strcmp(tors.atom_3.type,"H.spc")) strcpy(tors.atom_3.type,"H");
	if(!strcmp(tors.atom_4.type,"H.spc")) strcpy(tors.atom_4.type,"H");
	
	if(!strcmp(tors.atom_1.type,"N.4")) strcpy(tors.atom_1.type,"N.3");
	if(!strcmp(tors.atom_2.type,"N.4")) strcpy(tors.atom_2.type,"N.3");
	if(!strcmp(tors.atom_3.type,"N.4")) strcpy(tors.atom_3.type,"N.3");
	if(!strcmp(tors.atom_4.type,"N.4")) strcpy(tors.atom_4.type,"N.3");

	if(!strcmp(tors.atom_1.type,"O.co2")) strcpy(tors.atom_1.type,"O.2");
	if(!strcmp(tors.atom_2.type,"O.co2")) strcpy(tors.atom_2.type,"O.2");
	if(!strcmp(tors.atom_3.type,"O.co2")) strcpy(tors.atom_3.type,"O.2");
	if(!strcmp(tors.atom_4.type,"O.co2")) strcpy(tors.atom_4.type,"O.2");

	mark=0;

	for(i=0;i<num_torstype;i++)
		{
		 if(strcmp(tors.type,torsion[i].type)) continue;
		 else if(!strcmp(tors.atom_2.type,torsion[i].atom_2)&&
			 !strcmp(tors.atom_3.type,torsion[i].atom_3))
			{
			 // find the most suitable torsion parameter
			 tmp=2;
			 if(!strcmp(tors.atom_1.type,torsion[i].atom_1)) tmp++;
			 if(!strcmp(tors.atom_4.type,torsion[i].atom_4)) tmp++;
			 if(tmp>mark)
				{
				 mark=tmp;
				 tors.V=torsion[i].V;
				 tors.n=torsion[i].n;
				 tors.S=torsion[i].S;
				 if(mark==4) break;
				}
			}
		 else if(!strcmp(tors.atom_3.type,torsion[i].atom_2)&&
                         !strcmp(tors.atom_2.type,torsion[i].atom_3))
                        {
			 // find the most suitable torsion parameter
                         tmp=2;
                         if(!strcmp(tors.atom_4.type,torsion[i].atom_1)) tmp++;
                         if(!strcmp(tors.atom_1.type,torsion[i].atom_4)) tmp++;
                         if(tmp>mark)
                                {
                                 mark=tmp;
                                 tors.V=torsion[i].V;
                                 tors.n=torsion[i].n;
                                 tors.S=torsion[i].S;
				 if(mark==4) break;
                                }

                        }
		 else continue;
		}

	if(mark==0)
		{
                 printf("Warning: no suitable torsion parameter for ");
		 printf("%s  %s - %s - %s %s\n", 
			 tors.atom_1.type, tors.atom_2.type, tors.type,
			 tors.atom_3.type, tors.atom_4.type);
		 printf("Default value is assigned.\n");
               	 tors.V=0.5;   // default value for C.3-C.3-C.3-C.3
	 	 tors.n=3;
	 	 tors.S=1;
		}

        // printf("Torsion V: %5.2f n: %d  S %d\n", tors.V, tors.n, tors.S);

        return;
}

float ForceField::Cal_Torsion_Energy(Torsion tors) const
{
	float tmp1,tmp2;

	// TORSION_E = 1/2*V*[1+S*cos(n*angle)] 

	tmp1=0.50*tors.V;
	tmp2=tors.angle/180.0*3.1416*tors.n;

	tors.e=tmp1*(1+tors.S*cos((double)tmp2));
	
	return tors.e;
}

float ForceField::Cal_VDW_Energy(Atom a1, Atom a2, int mark_1_4) const
{
	float d0,d,eps,e;
	float tmp1,tmp2,tmp3,tmp4;

	// VDW_E = SQRT(e1*e2)*[(d0/d)^12-2.0*(d0/d)^6] 

	e=0.000;

	d=Distance(a1.coor,a2.coor);
	if(d>(VDW_DIST_CUTOFF)) return 0.000;

	d0=a1.r+a2.r;

	eps=sqrt((double)(a1.eps*a2.eps));

	tmp1=d0/d; tmp2=tmp1*tmp1*tmp1;
	tmp3=tmp2*tmp2; tmp4=tmp3*tmp3;

	e=eps*(tmp4-2.0*tmp3);

	if(mark_1_4) e*=(VDW_1_4_FACTOR);

	if(e>1000.0) e=1000.0;	

	return e;
}