grow.c

来自「药物开发中的基于结构的从头设计代码」· C语言 代码 · 共 1,238 行 · 第 1/3 页

        bond[num_bond].atom_1=result.id;
        bond[num_bond].atom_2=result_1.id;
        bond[num_bond].valid=1;
        strcpy(bond[num_bond].type,"1");

        bond[num_bond+1].id=num_bond+2;
        bond[num_bond+1].atom_1=result.id;
        bond[num_bond+1].atom_2=result_2.id;
        bond[num_bond+1].valid=1;
        strcpy(bond[num_bond+1].type,"1");

	num_bond+=2;

/*
        printf("Bridging happens between %d and %d\n", 
		root_1.center.id, root_2.center.id);
*/

        return TRUE;
}

int Ligand::Joining(int id1, int id2)
{
        extern ForceField *ff;
        Atom seed_1,seed_2;
        Atom result_1,result_2;
        Group root_1, root_2;
        int i;
        float tmp1,tmp2,tmp3;
        Bond new_bond;

        // find the two seed hydrogens and their roots

        seed_1=atom[id1-1]; seed_2=atom[id2-1];

        if(Find_Root_Of_A_Hydrogen(id1,root_1)==FALSE) return FALSE;
        if(Find_Root_Of_A_Hydrogen(id2,root_2)==FALSE) return FALSE;

	if(Connection_1_2_Check(root_1.center.id,root_2.center.id)==TRUE) 
		return FALSE;

	strcpy(new_bond.type,New_Bond_Viability_Check(root_1,root_2));

	if(!strcmp(new_bond.type,"no")) return FALSE;	
        else
                {
		 new_bond.id=num_bond+1;
                 new_bond.atom_1=root_1.center.id;
                 new_bond.atom_2=root_2.center.id;
                 new_bond.length=ff->Get_Bond_Length(root_1.center.type,
                                                     root_2.center.type,
                                                     new_bond.type);
		 new_bond.valid=1;
                }

        tmp1=Distance(seed_1.coor,root_1.center.coor);
	tmp2=Distance(seed_2.coor,root_2.center.coor);
        tmp3=new_bond.length;

        for(i=0;i<3;i++)
        {
         result_1.coor[i]=root_1.center.coor[i];
         result_1.coor[i]+=((seed_1.coor[i]-root_1.center.coor[i])*tmp3/tmp1);
         result_2.coor[i]=root_2.center.coor[i];
         result_2.coor[i]+=((seed_2.coor[i]-root_2.center.coor[i])*tmp3/tmp2);
        }

	if(Atom_Overlap_Check(result_1,root_2.center)==FALSE) return FALSE;
	if(Atom_Overlap_Check(result_2,root_1.center)==FALSE) return FALSE;

        // the new bond may need to change the atom type of root_1 or root_2

        atom[root_1.center.id-1]=root_1.center;
        atom[root_2.center.id-1]=root_2.center;

        // delete the two hydrogens

        Delete_An_Atom(id1); Delete_An_Atom(id2);

	// add the new bond

	bond[num_bond]=new_bond;

        num_bond++;

/*
        printf("Linking happens between %d and %d\n", 
		root_1.center.id, root_2.center.id);
*/

        return TRUE;
}

int Ligand::Link_Structure(Linking_Record &record)
{
	int i,j,mark;
	float d,d0;
	int old_num_bond,num_candidate;
	struct Candidate
		{
		 int valid;
		 int id1;
		 int id2;
		} candidate[500];

	old_num_bond=num_bond;

	record.num_link=0; record.rmsd=0.000;

	// check whether there is any atom in overlap 
		
	for(i=0;i<num_atom-1;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].type[0]=='H') continue;

		 for(j=i+1;j<num_atom;j++)
		{
		 if(atom[j].valid==0) continue;
		 else if(atom[j].type[0]=='H') continue;
		 else if(Atom_Overlap_Check(atom[i],atom[j])==FALSE) continue;
		 else return FALSE;
		}
		}

	// apply bridging and joining algorithm

	while(1)
	{
	 // first, make the candidate list of H-H bumps

	 Detect_Connections();

	 num_candidate=0;

	 for(i=0;i<500;i++) 
		{
		 candidate[i].valid=0; 
		 candidate[i].id1=0;
		 candidate[i].id2=0;
		}

	 for(i=0;i<num_atom-1;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].type[0]!='H') continue;

		 for(j=i+1;j<num_atom;j++)
		{
		 if(atom[j].valid==0) continue;
		 else if(atom[j].type[0]!='H') continue;

		 d0=atom[i].r+atom[j].r;
		 d=Distance(atom[i].coor,atom[j].coor);

		 if(d>(d0-VDW_BUMP_RANGE)) continue;
		 else if(Connection_1_2_Check(atom[i].id,atom[j].id)) continue;
		 else if(Connection_1_3_Check(atom[i].id,atom[j].id)) continue;
			{
			 candidate[num_candidate].valid=1;
			 candidate[num_candidate].id1=atom[i].id;
			 candidate[num_candidate].id2=atom[j].id;
			 num_candidate++;
			}
		}
		}

	 if(num_candidate==0) break;

	 // try to deal with each candidate

	 mark=FALSE;

	 for(i=0;i<num_candidate;i++)
		{
		 if(Bridging(candidate[i].id1,candidate[i].id2)==TRUE) 
			{
			 record.id1[record.num_link]=candidate[i].id1;
			 record.id2[record.num_link]=candidate[i].id2;
			 d=Distance(atom[candidate[i].id1-1].coor,
				    atom[candidate[i].id2-1].coor);
			 record.rmsd+=(d*d);
			 record.num_link++;
			 mark=TRUE;
			 break;
			}
		 else if(Joining(candidate[i].id1,candidate[i].id2)==TRUE)
			{
                         record.id1[record.num_link]=candidate[i].id1;
                         record.id2[record.num_link]=candidate[i].id2;
                         d=Distance(atom[candidate[i].id1-1].coor,
                                    atom[candidate[i].id2-1].coor);
                         record.rmsd+=(d*d);
                         record.num_link++;
			 mark=TRUE;
			 break;
			}
		 else continue;
		}

	if(mark==TRUE&&num_candidate==1) break;		// finish linking
	else if(mark==TRUE&&num_candidate!=1) continue;	// continue linking
	else if(mark==FALSE&&num_candidate!=0) return FALSE; // bad structure
	else break; // no more candidate, finish linking
	}

	if(num_bond==old_num_bond) return FALSE;	// no link made 

	record.rmsd=sqrt((double)(record.rmsd/record.num_link));

	this->valid=Arrange_IDs();

	if(this->valid==FALSE) return FALSE;
	else return TRUE;
}

int Ligand::Link_Structure()
{
	int i,j,mark;
	float d,d0;
	int old_num_bond,num_candidate;
	struct Candidate
		{
		 int valid;
		 int id1;
		 int id2;
		} candidate[500];

	old_num_bond=num_bond;

	// check whether there is any atom in overlap 
		
	for(i=0;i<num_atom-1;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].type[0]=='H') continue;

		 for(j=i+1;j<num_atom;j++)
		{
		 if(atom[j].valid==0) continue;
		 else if(atom[j].type[0]=='H') continue;
		 else if(Atom_Overlap_Check(atom[i],atom[j])==FALSE) continue;
		 else return FALSE;
		}
		}

	// apply bridging and joining algorithm

	while(1)
	{
	 // first, make the candidate list of H-H bumps

	 Detect_Connections();

	 num_candidate=0;

	 for(i=0;i<500;i++) 
		{
		 candidate[i].valid=0; 
		 candidate[i].id1=0;
		 candidate[i].id2=0;
		}

	 for(i=0;i<num_atom-1;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].type[0]!='H') continue;

		 for(j=i+1;j<num_atom;j++)
		{
		 if(atom[j].valid==0) continue;
		 else if(atom[j].type[0]!='H') continue;

		 d0=atom[i].r+atom[j].r;
		 d=Distance(atom[i].coor,atom[j].coor);

		 if(d>(d0-VDW_BUMP_RANGE)) continue;
		 else if(Connection_1_2_Check(atom[i].id,atom[j].id)) continue;
		 else if(Connection_1_3_Check(atom[i].id,atom[j].id)) continue;
			{
			 candidate[num_candidate].valid=1;
			 candidate[num_candidate].id1=atom[i].id;
			 candidate[num_candidate].id2=atom[j].id;
			 num_candidate++;
			}
		}
		}

	 if(num_candidate==0) break;	// no more possible linking

	 // try to deal with each candidate

	 mark=FALSE;

	 for(i=0;i<num_candidate;i++)
		{
		 if(Bridging(candidate[i].id1,candidate[i].id2)==TRUE) 
			{
			 mark=TRUE;
			 break;
			}
		 else if(Joining(candidate[i].id1,candidate[i].id2)==TRUE)
			{
			 mark=TRUE;
			 break;
			}
		 else continue;
		}

	if(mark==TRUE&&num_candidate==1) break;		// finish linking
	else if(mark==TRUE&&num_candidate!=1) continue;	// continue linking
	else if(mark==FALSE&&num_candidate!=0) return FALSE; // bad structure
	else break; // no more candidate, finish linking
	}

	if(num_bond==old_num_bond) return FALSE;	// no link made 

	this->valid=Arrange_IDs();

        if(this->valid==FALSE) return FALSE;
        else return TRUE;
}

void Ligand::Seed_Structure_Check() const
{
        extern Pocket *pok;
        int i,mark,*seed_list;
        char grid;

	if(num_part!=1)
	{
	 printf("\nThere are %d fragments in the seed structure\n",num_part);
	 printf("You are suggested to use LigBuilder/Link for instead.\n");
	 exit(1);
	}

	// check whether the seed is in the box

	mark=0;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue;
		 else if(atom[i].coor[0]>pok->max_x) mark++;
		 else if(atom[i].coor[0]<pok->min_x) mark++;
		 else if(atom[i].coor[1]>pok->max_y) mark++;
		 else if(atom[i].coor[1]<pok->min_y) mark++;
		 else if(atom[i].coor[2]>pok->max_z) mark++;
		 else if(atom[i].coor[2]<pok->min_z) mark++;
		 else continue;
		}

	if(mark)
		{
		 puts("\nError: part of the seed structure is out of the box.");
		 puts("Are you sure it has been docked to the right place?");
		 exit(1);
		}

	// then check whether any part of the seed bumps to the protein

	mark=0;

        for(i=0;i<num_atom;i++)
                {
                 grid=pok->Get_Grid_Property(atom[i].coor);
                 if(grid=='E')
                        {
                         if(!strcmp(atom[i].type,"H")) continue;
                         else if(!strcmp(atom[i].type,"H.spc"))
                                {strcpy(atom[i].type,"H"); continue;}
                         else
                                {
                                 printf("Atom %d %s bumps to the protein!\n",
                                        atom[i].id, atom[i].type);
                                 mark++;
                                 continue;
                                }
                        }
                 else continue;
                }

        if(mark) 
                {
                 puts("\nError: some atoms in the seed bump to the protein.");
		 puts("Please modify the seed structure to release the bumps.");
                 exit(1);
                }

	// notice here the seeds on 'S' grid are allowed to survive 

	// check whether there is any growing site on the seed

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

        mark=Make_Seed_List(seed_list);

        if(mark==0)
                {
                 puts("\nError: No available growing site found on the seed.");
                 puts("Program has no way to develop new molecules.");
                 puts("Please assign proper growing sites on the seed.");
                 exit(1);
                }

        delete [] seed_list;

	return;
}