# include "link.h"

void Ligand::Calculate_Binding_Score()
{
	extern Pocket *pok;
	int i;
	float vb,mb,shb,mhb,whb,swh,mwh,wwh,hm,rt;
	char grid;

	Calculate_HB_Root();

	Count_Rotors();

	binding_score=2.254;	// according to SCORE

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) {atom[i].score=0.000; continue;}
		 else if(!strcmp(atom[i].type,"H")) {atom[i].score=0.000; continue;}
		 else if(!strcmp(atom[i].type,"H.spc")) {atom[i].score=0.000; continue;}
		 else
			{
			 grid=pok->Get_Grid_Property(atom[i].coor);
			 if(grid=='E') continue;

			 vb=Get_An_Atom_VDW_Score(atom[i]);
			 mb=Get_An_Atom_MB_Score(atom[i]);
			 Get_An_Atom_HB_Score(atom[i],shb,mhb,whb,swh,mwh,wwh);
			 hm=Get_An_Atom_HM_Score(atom[i]);
			 rt=Get_An_Atom_RT_Score(atom[i]);

			 atom[i].score=0.000;
			 atom[i].score+=(SCORE_VB)*vb;
			 atom[i].score+=(SCORE_MB)*mb;
			 atom[i].score+=(SCORE_SHB)*shb;
			 atom[i].score+=(SCORE_MHB)*mhb;
			 atom[i].score+=(SCORE_WHB)*whb;
			 atom[i].score+=(SCORE_SWH)*swh;
			 atom[i].score+=(SCORE_MWH)*mwh;
			 atom[i].score+=(SCORE_WWH)*wwh;
			 atom[i].score+=(SCORE_HM)*hm;
			 atom[i].score+=(SCORE_RT)*rt;

			 binding_score+=atom[i].score;
			}
		}

	return;
}

void Ligand::Calculate_Chemical_Score()
{
	extern Parameter *parm;
	extern Pocket *pok;
	int i;
	char grid;

	chemical_score=0.00;

	// first, check whether the ligand has got out of the pocket 

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

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

		 if(grid=='E'||grid=='S') chemical_score-=10.0;
		 else continue;
		}

	// second, check whether there are too many chiral carbons :-)

	for(i=0;i<num_atom;i++)
        	{
         	 if(atom[i].valid==0) continue;
         	 else if(strcmp(atom[i].type,"C.3")) continue;
	 	 else if(atom[i].num_nonh<3) continue;
	 	 else chemical_score-=20.0;
        	}

	// third, check the number of components in the molecule

	chemical_score-=((num_part-1)*50.0);

	// fourth, penalize the fused rings

	chemical_score-=(num_fused_ring*30.0);

	// finally, Lipinski rules
	// notice, they are applied only when molecules are large enough.

	if(!strcmp(parm->apply_chemical_rules,"YES"))
	{
	 if(weight>parm->max_weight) chemical_score-=1000.0;
	 if(num_donor>parm->max_donor) chemical_score-=50.0;
	 if(num_acceptor>parm->max_acceptor) chemical_score-=50.0;
	 if(logp>parm->max_logp) chemical_score-=20.0;
	 if(binding_score>parm->max_pkd) chemical_score-=20.0;
	}
	 
	return;
}

void Ligand::Count_Rotors()
// if a single bond is common, bond.valid=1;
// if a single bond is a rotor, bond.valid=2;
{
	int i,mark;
	int id1,id2;
	char type1[10],type2[10];

	// first, eliminate all other bonds than not_in_ring single bonds

	for(i=0;i<num_bond;i++)
		{
		 if(bond[i].ring==1) continue;
		 else if(strcmp(bond[i].type,"1")) continue;
		 else bond[i].valid=2;
		}

	// second, eliminate all the R-H bonds

	for(i=0;i<num_bond;i++)
		{
		 if(bond[i].valid!=2) continue;
		 else
		{
		 id1=bond[i].atom_1; id2=bond[i].atom_2;
		 if(atom[id1-1].type[0]=='H'||
		    atom[id2-1].type[0]=='H') bond[i].valid=1;
		}
		}

	// third, eliminate terminal rotors and sp2-sp2 rotors

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

		 id1=bond[i].atom_1; id2=bond[i].atom_2;

		 // detect sp2-sp2 rotors

		 mark=0; 
		 strcpy(type1,atom[id1-1].type); 
		 strcpy(type2,atom[id2-1].type);

		 if(!strcmp(type1,"C.2")) mark++;
		 else if(!strcmp(type1,"C.ar")) mark++;
		 else if(!strcmp(type1,"C.1")) mark++;
		 else if(!strcmp(type1,"C.cat")) mark++;
		 else if(!strcmp(type1,"N.2")) mark++;

		 if(!strcmp(type2,"C.2")) mark++;
                 else if(!strcmp(type2,"C.ar")) mark++;
                 else if(!strcmp(type2,"C.1")) mark++;
                 else if(!strcmp(type2,"C.cat")) mark++;
                 else if(!strcmp(type2,"N.2")) mark++;

		 if(mark==2) {bond[i].valid=1; continue;}

		 // detect terminal rotors

		 if(atom[id1-1].num_nonh==1||atom[id2-1].num_nonh==1)
			bond[i].valid=1;
		}

	// here the false rotors in  -PO3, -CF3, -CMe3, -NMe3 are not checked 

	return;
}

int Ligand::Hydrogen_Bond_Check(Atom atom_1, Atom atom_2) const
{
	float d,d0,theta1,theta2;
	float v1[3],v2[3];
	float donor[3],d_root[3],acceptor[3],a_root[3];
	float dis_cut,angle_cut;
	int i,mark=0;
	Atom tmp_atom;

	if((!strcmp(atom_1.hb,"N"))||(!strcmp(atom_2.hb,"N"))) return FALSE;
	else if((!strcmp(atom_1.hb,"M"))||
		(!strcmp(atom_2.hb,"M"))) return FALSE;
	else if((!strcmp(atom_1.hb,"D"))&&
		(!strcmp(atom_2.hb,"D"))) return FALSE;
	else if((!strcmp(atom_1.hb,"A"))&&
		(!strcmp(atom_2.hb,"A"))) return FALSE;

	d0=atom_1.R+atom_2.R; dis_cut=d0+0.40; 
	d=Distance(atom_1.coor,atom_2.coor); 

	if(d>dis_cut) return FALSE;	
	else if(d>4.00) return FALSE;

	if(strcmp(atom_1.hb,atom_2.hb)>0)
		{
		 tmp_atom=atom_1;
		 atom_1=atom_2;
		 atom_2=tmp_atom;
		}
	
	if(!strcmp(atom_1.hb,"A")&&!strcmp(atom_2.hb,"D"))
		{
		 acceptor[0]=atom_1.coor[0]; a_root[0]=atom_1.root[0];
		 acceptor[1]=atom_1.coor[1]; a_root[1]=atom_1.root[1];
		 acceptor[2]=atom_1.coor[2]; a_root[2]=atom_1.root[2];
		 donor[0]=atom_2.coor[0]; d_root[0]=atom_2.root[0];
		 donor[1]=atom_2.coor[1]; d_root[1]=atom_2.root[1];
		 donor[2]=atom_2.coor[2]; d_root[2]=atom_2.root[2];

		 if(!strcmp(atom_2.type,"O.w")) mark=1; // D is water
		 else if(!strcmp(atom_1.type,"O.w")) mark=2; // A is water
		 else mark=0; // none is water, normal h-bond
		}
	else if(!strcmp(atom_1.hb,"A")&&!strcmp(atom_2.hb,"DA"))
		{
		 acceptor[0]=atom_1.coor[0]; a_root[0]=atom_1.root[0];
		 acceptor[1]=atom_1.coor[1]; a_root[1]=atom_1.root[1];
		 acceptor[2]=atom_1.coor[2]; a_root[2]=atom_1.root[2];
		 donor[0]=atom_2.coor[0]; d_root[0]=atom_2.root[0];
		 donor[1]=atom_2.coor[1]; d_root[1]=atom_2.root[1];
		 donor[2]=atom_2.coor[2]; d_root[2]=atom_2.root[2];

		 if(!strcmp(atom_2.type,"O.w")) mark=1; // D is water
		 else if(!strcmp(atom_1.type,"O.w")) mark=2; // A is water
		 else mark=0; // none is water, normal h-bond
		}
	else if(!strcmp(atom_1.hb,"D")&&!strcmp(atom_2.hb,"DA"))
		{
		 donor[0]=atom_1.coor[0]; d_root[0]=atom_1.root[0];
		 donor[1]=atom_1.coor[1]; d_root[1]=atom_1.root[1];
		 donor[2]=atom_1.coor[2]; d_root[2]=atom_1.root[2];
		 acceptor[0]=atom_2.coor[0]; a_root[0]=atom_2.root[0];
		 acceptor[1]=atom_2.coor[1]; a_root[1]=atom_2.root[1];
		 acceptor[2]=atom_2.coor[2]; a_root[2]=atom_2.root[2];

		 if(!strcmp(atom_1.type,"O.w")) mark=1; // D is water
		 else if(!strcmp(atom_2.type,"O.w")) mark=2; // A is water
		 else mark=0; // none is water, normal h-bond	
		} 		
	else if(!strcmp(atom_1.hb,"DA")&&!strcmp(atom_2.hb,"DA"))
		{
		 donor[0]=atom_1.coor[0]; d_root[0]=atom_1.root[0];
		 donor[1]=atom_1.coor[1]; d_root[1]=atom_1.root[1];
		 donor[2]=atom_1.coor[2]; d_root[2]=atom_1.root[2];
		 acceptor[0]=atom_2.coor[0]; a_root[0]=atom_2.root[0];
		 acceptor[1]=atom_2.coor[1]; a_root[1]=atom_2.root[1];
		 acceptor[2]=atom_2.coor[2]; a_root[2]=atom_2.root[2];

		 if(!strcmp(atom_1.type,"O.w")) mark=1; // D is water
		 else if(!strcmp(atom_2.type,"O.w")) mark=2; // A is water
		 else mark=0; // none is water, normal h-bond	
		}
	else return FALSE;

	if(mark==1) {theta1=180.0;}
	else
		{
		 for(i=0;i<3;i++)
			{
		 	 v1[i]=donor[i]-d_root[i];
		 	 v2[i]=donor[i]-acceptor[i];
			}

		 theta1=Angle_Of_Two_Vectors(v1,v2);  // angle among DR-D-A
		}
	

	if(mark==2) {theta2=180.0;}
	else
		{
		 for(i=0;i<3;i++)
			{
			 v1[i]=acceptor[i]-donor[i];
			 v2[i]=acceptor[i]-a_root[i];
			}
	
		 theta2=Angle_Of_Two_Vectors(v1,v2);  // angle among D-A-AR
		}

	angle_cut=(HB_ANGLE_CUTOFF);		

	if(mark==0)
		{
		 if(theta1<angle_cut) return FALSE;
		 else if(theta2<angle_cut) return FALSE;
		 else return 1;		// ordinary h-bond
		}
	else if(mark==1)
		{
		 if(theta2<angle_cut) return FALSE;
		 else return 2;		// water-involved h-bond
		}
	else if(mark==2)
		{
		 if(theta1<angle_cut) return FALSE;
		 else return 2;		// water-involved h-bond
		}
	else return FALSE;
}

float Ligand::Get_An_Atom_MB_Score(Atom atm) const
{
	extern Pocket *pok;
	int i;
	float d,tmp;

	float mb=0.000;
	
	if(strcmp(atm.hb,"A")&&strcmp(atm.hb,"DA")) return mb;

	for(i=0;i<pok->num_atom;i++)
		{
		 if(strcmp(pok->atom[i].hb,"M")) continue;
		 else
			{
			 d=Distance(atm.coor,pok->atom[i].coor);
			 if(d<2.20) tmp=1.00;
		 	 else if(d<2.70) tmp=5.40-2*d;
		 	 else tmp=0.00;
		 	 mb+=tmp; 
			}
		}

	return mb;
}

float Ligand::Get_An_Atom_VDW_Score(Atom atm) const
{
	extern Pocket *pok;
	int i;
	float d,d0;

	float vb=0.000;

	if(!strcmp(atm.type,"H")||!strcmp(atm.type,"H.spc")) return vb;

	for(i=0;i<pok->num_atom;i++)
		{
		 if(!strcmp(pok->atom[i].type,"H")) continue;
		 else
			{
			 d=Distance(atm.coor,pok->atom[i].coor);
			 d0=atm.R+pok->atom[i].R;
			 if(d<(d0-0.60)) vb+=1.00;
			}
		}

	return vb;
}

void Ligand::Get_An_Atom_HB_Score(Atom atm, 
				  float &shb, float &mhb, float &whb,
				  float &swh, float &mwh, float &wwh) const
{
	extern Pocket *pok;
	int i,j,mark,num_candidate;
	float d,d0;
	struct 
		{
		 int valid;
		 int type;
		 float length;
		 float overlap;
		} candidate[10],tmp;

	for(i=0;i<10;i++)
		{
		 candidate[i].valid=0;
		 candidate[i].type=0;
		 candidate[i].length=0.000;
		 candidate[i].overlap=0.000;
		}

	shb=mhb=whb=swh=mwh=wwh=0.000;

	if(!strcmp(atm.hb,"N")) return;  

	// first, get all the possible hydrogen bonds

	num_candidate=0;

	for(i=0;i<pok->num_atom;i++)
		{
		 if(!strcmp(pok->atom[i].hb,"N")) continue;

		 mark=Hydrogen_Bond_Check(atm,pok->atom[i]);

		 if(mark==FALSE) continue;

		 d=Distance(atm.coor,pok->atom[i].coor);
		 d0=atm.R+pok->atom[i].R;

		 candidate[num_candidate].valid=1;
		 candidate[num_candidate].type=mark;
		 candidate[num_candidate].length=d;
		 candidate[num_candidate].overlap=d-d0;
		 num_candidate++;
		}

	// second, rank all the candidate h-bonds in terms of length

	for(i=0;i<num_candidate-1;i++)
	for(j=i+1;j<num_candidate;j++)
		{
		 if(candidate[i].length<=candidate[j].length) continue;
		 else
			{
			 tmp=candidate[i];
			 candidate[i]=candidate[j];
			 candidate[j]=tmp;
			}
		}
	
	// only keep the top three h-bonds

	for(i=3;i<10;i++) candidate[i].valid=0;
		
	// finally, count all the h-bonds 

	for(i=0;i<num_candidate;i++)
		{
		 if(candidate[i].valid==0) continue;
		 else if(candidate[i].type==1)
			{
			 if(candidate[i].overlap<-0.60) shb+=1.00;
			 else if(candidate[i].overlap<-0.30) mhb+=1.00;
			 else whb+=1.00;
			}
		 else if(candidate[i].type==2)
			{
			 if(candidate[i].overlap<-0.60) swh+=1.00;
			 else if(candidate[i].overlap<-0.20) mwh+=1.00;
			 else wwh+=1.00;
			}
		 else continue;
		}
	
	return;
}

float Ligand::Get_An_Atom_HM_Score(Atom atm) const
{
	extern Pocket *pok;
	int i,neib;
	float d,tmp;

	float hm=0.000;

	if(atm.logp<0.20) return hm;

	tmp=0.000; neib=0;

	for(i=0;i<pok->num_atom;i++)
		{
		 if(!strcmp(pok->atom[i].type,"H")) continue;
		 else
			{
			 d=Distance(atm.coor,pok->atom[i].coor);
			 if(d<5.00) {tmp+=pok->atom[i].logp;neib++;}
			}
		}

	if(neib>=2&&tmp>-1.00) hm=atm.logp;

	return hm;
}

float Ligand::Get_An_Atom_RT_Score(Atom atm) const
{
	int i;
	float rt=0.000;

	for(i=0;i<num_bond;i++)
		{
		 if(bond[i].valid!=2) continue;
		 else if((bond[i].atom_1==atm.id)||
			 (bond[i].atom_2==atm.id))
		 	{rt+=0.500; continue;}
		 else continue;
		}

	return rt;
}

void Ligand::Calculate_HB_Root()
{
	int i,j,tmp,num_nonh;
	float tmpx,tmpy,tmpz;

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

		 tmpx=tmpy=tmpz=0.000; num_nonh=0;
		
		 for(j=0;j<atom[i].num_neib;j++)
			{
			 tmp=atom[i].neib[j];
			 if(atom[tmp-1].type[0]=='H') continue;
			 else
				{
				 tmpx+=atom[tmp-1].coor[0];
				 tmpy+=atom[tmp-1].coor[1];
				 tmpz+=atom[tmp-1].coor[2];
				 num_nonh++;
				}
			}

		 if(num_nonh==0) strcpy(atom[i].hb,"N");
		 else
			{
		 	 tmpx/=num_nonh;
			 tmpy/=num_nonh;
			 tmpz/=num_nonh;
			 atom[i].root[0]=tmpx;
			 atom[i].root[1]=tmpy;
			 atom[i].root[2]=tmpz;
			}
		}

	return;
}

int Ligand::Chemical_Viability_Check()
{
	extern Parameter *parm;
	extern FragLibrary *toxiclib; 

	// first, Lipinski rules 

	if(!strcmp(parm->apply_chemical_rules,"YES"))
	{
	 if(weight>parm->max_weight) return FALSE;
	 else if(weight<parm->min_weight) return FALSE;
	 else if(logp>parm->max_logp) return FALSE;
	 else if(logp<parm->min_logp) return FALSE;
	 else if(num_donor>parm->max_donor) return FALSE;
	 else if(num_donor<parm->min_donor) return FALSE;
	 else if(num_acceptor>parm->max_acceptor) return FALSE;
	 else if(num_acceptor<parm->min_acceptor) return FALSE;
	 else if(binding_score>parm->max_pkd) return FALSE;
	 else if(binding_score<parm->min_pkd) return FALSE;
	 else {}
	}

	// second, toxic substructure search

	if(!strcmp(parm->apply_toxicity_check,"YES"))
		{
		 if(toxiclib->Map_A_Molecule(*this)!=0) return FALSE;
		}

	return TRUE;
}