⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 ligand.c

📁 药物开发中的基于结构的从头设计代码
💻 C
📖 第 1 页 / 共 4 页
字号:
	// then assign atomic weight to atoms
	// it is necessary to rank a certain atom's neighbors

        for(i=0;i<num_atom;i++)
                {
                 if(!strcmp(atom[i].type,"Cl")) atom[i].weight=35;
                 else if(!strcmp(atom[i].type,"Br")) atom[i].weight=80;
		 else if(atom[i].type[0]=='C') atom[i].weight=12;
		 else if(atom[i].type[0]=='H') atom[i].weight=1;
		 else if(atom[i].type[0]=='N') atom[i].weight=14;
		 else if(atom[i].type[0]=='O') atom[i].weight=16;
		 else if(atom[i].type[0]=='P') atom[i].weight=31;
		 else if(atom[i].type[0]=='S') atom[i].weight=32;
		 else if(atom[i].type[0]=='F') atom[i].weight=19;
                 else if(atom[i].type[0]=='I') atom[i].weight=127;
		 else atom[i].weight=0;
                }

	// now detect the neighboring atoms for each atom

	for(i=0;i<num_atom;i++) atom[i].num_neib=atom[i].num_nonh=0;

	for(i=0;i<num_bond;i++)
		{
		 if(bond[i].valid==0) continue;

		 id1=bond[i].atom_1; id2=bond[i].atom_2;
		 tmp1=atom[id1-1].num_neib; tmp2=atom[id2-1].num_neib;

		 atom[id1-1].neib[tmp1]=id2; atom[id1-1].num_neib++;
		 atom[id2-1].neib[tmp2]=id1; atom[id2-1].num_neib++;
		}

	// now arrange the neighboring atoms in terms of atomic weights

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].num_neib<=1) continue;

		 for(j=0;j<atom[i].num_neib-1;j++)
		 for(k=j+1;k<atom[i].num_neib;k++)
			{
			 tmp1=atom[i].neib[j]; tmp2=atom[i].neib[k];
			 if(atom[tmp1-1].weight>=atom[tmp2-1].weight) continue;
			 else
				{
				 mark=atom[i].neib[j];
				 atom[i].neib[j]=atom[i].neib[k];
				 atom[i].neib[k]=mark;
				}
			}
		}

	for(i=0;i<num_atom;i++)
		{
		 for(j=0;j<atom[i].num_neib;j++)
			{
			 if(atom[atom[i].neib[j]-1].type[0]=='H') continue;
		 	 else atom[i].num_nonh++;
			}

		}

	// now detect the neighboring bonds

	for(i=0;i<num_bond;i++) bond[i].num_neib=0;

	for(i=0;i<num_bond-1;i++)
	for(j=i+1;j<num_bond;j++)
		{
		 tmp1=bond[i].num_neib; tmp2=bond[j].num_neib;

		 if(Two_Bonds_Connection_Check(bond[i],bond[j]))
			{
			 bond[i].neib[tmp1]=bond[j].id; bond[i].num_neib++;
			 bond[j].neib[tmp2]=bond[i].id; bond[j].num_neib++;
			}

		 else continue;
		}

	return;
}

void Ligand::Translate(float move[3])
{
	int i;

	for(i=0;i<num_atom;i++)
		{
		 atom[i].coor[0]+=move[0];
		 atom[i].coor[1]+=move[1];
		 atom[i].coor[2]+=move[2];
		}

	return;
}

void Ligand::Rotate(float angle,float axis[3],float origin[3])
{
	int i,j,k;
	float old_vector[3],new_vector[3];
	float mat[3][3];	
	float a,b;

	Unify_Vector(axis);

	a=sin((double)(angle*3.1416/180.0));
        b=cos((double)(angle*3.1416/180.0));

	mat[0][0]=axis[0]*axis[0]+(1-axis[0]*axis[0])*b;
        mat[1][0]=axis[0]*axis[1]*(1-b)-axis[2]*a;
        mat[2][0]=axis[0]*axis[2]*(1-b)+axis[1]*a;
        mat[0][1]=axis[0]*axis[1]*(1-b)+axis[2]*a;
        mat[1][1]=axis[1]*axis[1]+(1-axis[1]*axis[1])*b;
        mat[2][1]=axis[1]*axis[2]*(1-b)-axis[0]*a;
        mat[0][2]=axis[0]*axis[2]*(1-b)-axis[1]*a;
        mat[1][2]=axis[1]*axis[2]*(1-b)+axis[0]*a;
        mat[2][2]=axis[2]*axis[2]+(1-axis[2]*axis[2])*b;

	for(i=0;i<num_atom;i++)
		{
		 for(j=0;j<3;j++) 
			{
			 old_vector[j]=atom[i].coor[j]-origin[j];
			 new_vector[j]=0.000;
			}

		 for(j=0;j<3;j++)
			{
			 for(k=0;k<3;k++) 
				new_vector[j]+=(old_vector[k]*mat[k][j]);
			}

		 for(j=0;j<3;j++)
			{
			 atom[i].coor[j]=new_vector[j]+origin[j];
			}
		}

	return;
}

Group Ligand::Find_A_Group(int id) const
// id is the ID of the central atom
{
        int i,j,num;
        Group group;

        // define the center

        group.center=atom[id-1];

        // find the center's neighbours

        num=group.center.num_neib;

        for(i=0;i<num;i++)
                {
                 group.neib[i]=atom[group.center.neib[i]-1];

                 for(j=0;j<num_bond;j++)
                        {
                         if((bond[j].atom_1==group.center.neib[i])&&
                            (bond[j].atom_2==group.center.id))
                                {
                                 group.bond[i]=bond[j];
                                 break;
                                }
                         else if((bond[j].atom_2==group.center.neib[i])&&
                                 (bond[j].atom_1==group.center.id))
                                {
                                 group.bond[i]=bond[j];
                                 break;
                                }
                         else continue;
                        }
                }

        // count the necessary parameters

        group.num_neib=num;

        group.num_h=0; group.num_nonh=0;

        for(i=0;i<num;i++)
                {
                 if(group.neib[i].type[0]=='H') group.num_h++;
                 else group.num_nonh++;
                }

        group.num_hetero=0;

        for(i=0;i<group.num_nonh;i++)
                {
                 if(!strcmp(group.bond[i].type,"2")) continue;
                 else if(!strcmp(group.bond[i].type,"3")) continue;
                 else if(!strcmp(group.bond[i].type,"ar")) continue;
                 else if(group.neib[i].type[0]=='N') group.num_hetero++;
                 else if(group.neib[i].type[0]=='O') group.num_hetero++;
                 else continue;
                }

        group.num_pi=0;

        for(i=0;i<group.num_nonh;i++)
                {
                 if(!strcmp(group.bond[i].type,"2")) continue;
                 else if(!strcmp(group.bond[i].type,"3")) continue;
                 else if(!strcmp(group.neib[i].type,"C.ar")) group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"C.2")) group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"C.1"))group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"C.cat")) group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"N.2")) group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"N.1")) group.num_pi++;
                 else if(!strcmp(group.neib[i].type,"N.ar")) group.num_pi++;
                 else continue;
                }

        group.num_nar=group.num_car=0;

        for(i=0;i<group.num_nonh;i++)
                {
                 if(strcmp(group.bond[i].type,"ar")) continue;
                 else if(!strcmp(group.neib[i].type,"N.ar")) group.num_nar++;
                 else if(!strcmp(group.neib[i].type,"C.ar")) group.num_car++;
                 else if(group.neib[i].type[0]=='C') group.num_car++;
                 else group.num_nar++;
                }

        group.db_type=0;

        for(i=0;i<group.num_nonh;i++)
                {
                 if(!strcmp(group.center.type,"O.co2")||
                    !strcmp(group.center.type,"O.2")||
                    !strcmp(group.center.type,"S.2"))
                        {
                         if(group.neib[i].type[0]=='C') group.db_type=1;
                         else if(group.neib[i].type[0]=='N') group.db_type=2;
                         else if(group.neib[i].type[0]=='O') group.db_type=3;
                         else if(group.neib[i].type[0]=='S') group.db_type=4;
                         else if(group.neib[i].type[0]=='P') group.db_type=5;
                         else continue;
                        }
                 else if(!strcmp(group.bond[i].type,"2")||
                         !strcmp(group.neib[i].type,"O.co2")||
                         !strcmp(group.neib[i].type,"O.2")||
                         !strcmp(group.neib[i].type,"S.2"))
                        {
                         if(group.neib[i].type[0]=='C') group.db_type=1;
                         else if(group.neib[i].type[0]=='N') group.db_type=2;
                         else if(group.neib[i].type[0]=='O') group.db_type=3;
                         else if(group.neib[i].type[0]=='S') group.db_type=4;
                         else if(group.neib[i].type[0]=='P') group.db_type=5;
                         else continue;
                        }
                 else continue;
                }

        group.valid=1; return group;
}

int Ligand::Two_Bonds_Connection_Check(Bond bond1, Bond bond2) const
// this function returns the ID of the joint atom
{
        int id;

        id=bond1.atom_1;

        if(id==bond2.atom_1) return id;
        else if(id==bond2.atom_2) return id;

        id=bond1.atom_2;

        if(id==bond2.atom_1) return id;
        else if(id==bond2.atom_2) return id;

        return 0;       // two bonds are not connected
}

void Ligand::Delete_An_Atom(int atom_id)
// atom_id: the id of the hydrogen to be deleted
{
	int i;

	// first, set the atom to invalid

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].id!=atom_id) continue;
		 else {atom[i].valid=0; break;}
		}

	// second, set all the related bonds to invalid 

	for(i=0;i<num_bond;i++)
		{
		 if((bond[i].atom_1==atom_id)||(bond[i].atom_2==atom_id))
			{bond[i].valid=0;}
		 else continue;
		}

	return;

	// NOTICE: this function leaves the molecule in a dangerous way.
	// It must be followed by Arrange_ID() to re-arrange the structure.
}

int Ligand::Arrange_IDs()
// delete all the invalid atoms and bonds and re-arrange the ids 
// if return TRUE, operation okay
// if return FALSE, something wrong with the molecule
{
	int i,j,tmp,num_a,num_b;

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

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

	struct Record {int id_old; int id_new;} *record;
	record=new Record[num_atom];
	if(record==NULL) Memory_Allocation_Error();

	// find the left valid atoms

	num_a=0;	// number of left valid atoms

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue;
		 else
			{
			 record[num_a].id_old=atom[i].id;
			 record[num_a].id_new=num_a+1;
			 tmp_atom[num_a]=atom[i];
			 tmp_atom[num_a].id=num_a+1;
			 num_a++;
			}
		}

	// find the left valid bonds

	num_b=0;	// number of left valid bonds

	for(i=0;i<num_bond;i++)
	{
	 if(bond[i].valid==0) continue;
	 else
		{
		 tmp_bond[num_b]=bond[i];
		 tmp_bond[num_b].id=num_b+1;

		 tmp=bond[i].atom_1;
		 for(j=0;j<num_atom;j++)
			{
			 if(tmp!=record[j].id_old) continue;
			 else
			{
			 tmp_bond[num_b].atom_1=atom[record[j].id_new-1].id;
			 break;
			}
			}

		 tmp=bond[i].atom_2;
		 for(j=0;j<num_atom;j++)
                        {
                         if(tmp!=record[j].id_old) continue;
                         else
                        {
                         tmp_bond[num_b].atom_2=atom[record[j].id_new-1].id;
                         break;
                        }
                        }

		 num_b++;
		}
	}

	// check whether the molecule seems to be okay

	if(num_a<=0||num_b<=0) return FALSE;

	// generate the new list

	num_atom=num_a;	
	for(i=0;i<num_atom;i++) atom[i]=tmp_atom[i];

	num_bond=num_b;	
	for(i=0;i<num_bond;i++) bond[i]=tmp_bond[i];

	delete [] tmp_atom;
	delete [] tmp_bond;
	delete [] record;
	
	this->Detect_Connections();

	return TRUE;
}

int Ligand::Outer_Collision_Check() 
{
	extern Pocket *pok;
	int i,mark;
	char grid;

	mark=FALSE;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue; 
		 else if(!strcmp(atom[i].type,"H")) continue;

		 grid=pok->Get_Grid_Property(atom[i].coor);

		 if(grid=='E')
			{
			 if(!strcmp(atom[i].type,"H.spc"))
				{
				 strcpy(atom[i].type,"H"); 
				 continue;
				}
			 else 
				{
				 mark++; 
				 continue;
				}
			}
		 else if(grid=='S')
			{
			 if(!strcmp(atom[i].type,"H.spc"))
                                {
				 strcpy(atom[i].type,"H"); 
				 continue;
				}
			 else continue;
			}
		 else continue;
		}

	return mark;
}

int Ligand::Inner_Collision_Check() const
{
	int i,j;
	float d,d0;

	for(i=0;i<num_atom-1;i++)
	for(j=i+1;j<num_atom;j++)
		{

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -