# include "hbond.h"

/*
 * this routines takes three atoms and returns the cosine of 
 * angle between them. atom b is the vertex of the angle.
 */

/*     ==========                                                  */
double FindCosine (atom_record *a , atom_record *b, atom_record *c)
/*     ==========                                                  */
{
   /* find vectors */
  double
    xab = a->x - b->x ,
    yab = a->y - b->y ,
    zab = a->z - b->z ,

    xcb = c->x - b->x ,
    ycb = c->y - b->y ,
    zcb = c->z - b->z ;
    

  /* find dot product */
  double
    dot_product = xab*xcb + yab*ycb + zab*zcb ;

  /* find the square of lengths of vectors */
  double 
    r2ab = xab*xab + yab*yab + zab*zab , 
    r2cb = xcb*xcb + ycb*ycb + zcb*zcb ;
 
  /* evaluate to find cosine and return value */
  double 
    cosine_angle = dot_product/sqrt(r2ab*r2cb) ;

  return cosine_angle ;
}


/*                wrapper function */
/*     =========                                                    */   
double FindAngle (atom_record *a, atom_record * b, atom_record * c)
/*     =========                                                    */
{
# define DEG_ACOS(x) TODEG(acos(x))

  double
    cosine_angle = FindCosine(a, b, c) ;

  return DEG_ACOS(cosine_angle) ;
}


/*   ====================                                                     */
void CalculateHbondAngles
/*   ====================                                                     */
(atom_record *acc, atom_record *don, double *CNH_O_Angle, double *NH_OC_Angle)
{
  AtomList
    *ptr_don_bonded,
    *ptr_acc_bonded;


  /*              checking only the nearest neighbor in both cases          */

  /* finding the C-N-H..O angle, 90 to 150 deg */
  
  ptr_don_bonded = don->Contacts;
  if ( AreCovalentlyConnected(ptr_don_bonded->atom, ptr_don_bonded->Dist) ) 
    *CNH_O_Angle = FindAngle(ptr_don_bonded->atom, don, acc);


  /* finding the N-H..0-C angle, 180 to 110 deg = 110 deg threshhold*/

  ptr_acc_bonded = acc->Contacts;
  if ( AreCovalentlyConnected(ptr_acc_bonded->atom, ptr_acc_bonded->Dist) ) 
    *NH_OC_Angle = FindAngle(don, acc, ptr_acc_bonded->atom);
}



/*                    checks for oxygen only as acceptor                  */

/*   ============                                                  */    
bool IsAnAcceptor(atom_record * ac)
/*   ============                                                  */    
{
  if (ac->name[0] == 'O') 
    return true ; 
 
  else 
    return false ; 
}


/*
 *  checks if donor is nitrogen or oxygen. 
 *  only includes following oxygens as donors:
 *     - gamma oxygen on serine and threonine
 *     - eta oxygen or OH on tyrosines
 *     the following two due to difficulty in distinguishing between
 *       nitrogen and oxygen since they differ by only one electron
 *     - delta oxygen on asparagines
 *     - epsilon oxygen on glutamines
 */

/*   ========                                                  */    
bool IsADonor (atom_record * don)
/*   ========                                                  */    
{
  if (don->name[0] == 'N') 
    return true ;

  else if (don->name[0] == 'O' &&
	   (don->name[1] == 'G' || don->name[1] == 'H' ||
	    ((don->name[1] == 'E' || don->name[1] == 'D') && 
	     don->resname[2] == 'N'))) 
    return true ;
  
  else 
    return false ; 
}


/*
   As for your second question, I get covalent bonds by getting all distances
to a given atom i, sorting by increasing distance and then checking the shortestdistances. 
I take two atoms X as being bonded if the distance between them
is less than 1.7 Angstroms.  If one of the X is S (sulfur) I use 2 A and if
both are S I use 2.8 A.  If one of the X is H, I use 1.2 A.  In spite of 
its simplicity, this scheme has worked for thousands of molecules.
*/

/*   ======================                                 */    
bool AreCovalentlyConnected ( atom_record * cc, double Dist )
/*   ======================                                 */    
{
  /* didn't include other atom since in this code.                  */
  /* it must be either an O for donor or an N or O as and acceptor       */

  /* distance requirements */
  double
    d1S = 2.0,
    d1H = 1.2,
    rest = 1.7;

  
  if(cc->name[0] == 'S' && Dist < d1S)
    return true;
  
  else if(cc->name[0] == 'H' && Dist < d1H)
    return true;

  else if(Dist < rest)
    return true;
  
  else
    return false;
}



/*
Mikes says:

   For hydrogen bonds, I use a criterion based on distances and angles.
There is a rourtine in ~levitt/encad/encadv5 called hbonof.lc that finds
hydrogen bonds.  If there are no explicit hydrogen atoms, it uses the coordinates
of the acceptor (O), the donor (N) and an atom bonded to the donor (C).  The
conditions are d(O..N) < 3.6 A, angle (O..N-C) between 90 and 150 (I think).
Again, the code works very well.  There is another routine in ~levitt/vitrage/lc??
called hbonoh.lc that is more recent and cleaner.
*/

/*
   			from Arthur M. Lesk
I calculate the N-H..O-C angle (for backbone atoms, for instance)

I expect the angles to be distributed around 120 degrees for
a hydrogen bond to a carbonyl group.

In your case, you would expect a different "ideal" angle depending on whether
you were dealing with a quaternary nitrogen (as in lysine) or a tertiary
nitrogen (as in the distal nitrogens in arg), because the C-N-H angle differs
-- even if the actual hydrogen bond   N-H ...O  is linear. 

Currently default angle threshold    = 110. degrees
          default distance threshold = 3.5 A
*/


/*           tests Hbond distance; calculates and tests Hbond angle */

/*   ==========                                                     */
bool AreHbonded
/*   ==========                                                     */
(atom_record *acc, atom_record *don, double Dist, 
 double *CNH_O_Angle, double *NH_OC_Angle, bool flag)
{
# define HBOND_DIST_THRESHOLD 3.6

# define CNH_O_ANGLEMIN 90
# define CNH_O_ANGLEMAX 150

# define NH_CO_ANGLETHRESH 110


  /* someday put in code to test for HETATOMs ie water, carbonate, effectors */

  /* tests Hbonds for mainchain interactions */
  if ( (IsAnAcceptor(acc)) && 
      !(don->resnum == acc->resnum) &&
      !(AreCovalentlyConnected(don, Dist)) &&
       (IsADonor(don)) )  {
    
    /* test Hbond distance */
    if ( Dist > HBOND_DIST_THRESHOLD )
      return false;

    /* test Hbond angle */
    CalculateHbondAngles (acc, don, CNH_O_Angle, NH_OC_Angle);

    if (flag)  {
      if (*NH_OC_Angle >= NH_CO_ANGLETHRESH)
	return true;
      else 
	return false;
    }

    else if (*CNH_O_Angle >= CNH_O_ANGLEMIN && *CNH_O_Angle <= CNH_O_ANGLEMAX)
      return true;

    else 
      return false;
  }

  else
    return false;

# undef HBOND_DIST_THRESHOLD

# undef CNH_O_ANGLEMIN
# undef CNH_O_ANGLEMAX

# undef NH_CO_ANGLETHRESH
}


/*  ---------                                                       */
int SaveHbond
/*  ---------                                                       */
(atom_record * acc, atom_record * don, double Dist, 
 double CNH_O_Angle, double NH_OC_Angle, int n, Hbond_Data **Hbonds)
{
# define MAX_HBONDS 1000

  if (n == 0) 
    *Hbonds = (Hbond_Data *) calloc (MAX_HBONDS,sizeof(Hbond_Data));
 

  if (n >= MAX_HBONDS)  {
    STDERR("Hey! Too Many Hydrogen bonds for the array!");
    exit(1);
  }

  (*Hbonds)[n].acc = acc ; 
  (*Hbonds)[n].don = don ;
  (*Hbonds)[n].Dist = Dist ;
  (*Hbonds)[n].CNH_O_Angle = CNH_O_Angle ;
  (*Hbonds)[n].NH_OC_Angle = NH_OC_Angle ; 

  return ++n ;

# undef MAX_HBONDS
}

/*  --------                                         */
int HbondCmp (const void *v1, const void *v2)
/*  --------                                         */
{
  Hbond_Data *h1 = (Hbond_Data *) v1 ;
  Hbond_Data *h2 = (Hbond_Data *) v2 ;

  return (h1->mid->resnum - h2->mid->resnum) ;
}


/*         ===========                                     */
Hbond_Data *SortHbonds (Hbond_Data *Hbond, int n)
/*         ===========                                     */
{
  int
    count,
    b = 2*n;

  Hbond_Data *Big;

  Big = (Hbond_Data *) calloc (b, sizeof(Hbond_Data));

  for (count=0; count<n; count++)  {
    Big[count].mid = Hbond[count].don;
    Big[count].acc = Hbond[count].acc;
    Big[count].Dist = Hbond[count].Dist ;
    Big[count].CNH_O_Angle = Hbond[count].CNH_O_Angle ;
    Big[count].NH_OC_Angle = Hbond[count].NH_OC_Angle ;
    Big[count].type = 'a';
  }

  for (count; count<b; count++)  {
    Big[count].don = Hbond[count-n].don;
    Big[count].mid = Hbond[count-n].acc;
    Big[count].Dist = Hbond[count-n].Dist ;
    Big[count].CNH_O_Angle = Hbond[count-n].CNH_O_Angle ;
    Big[count].NH_OC_Angle = Hbond[count-n].NH_OC_Angle ;
    Big[count].type = 'd';
  }
  
  qsort((void *) Big,b,sizeof(Hbond_Data),HbondCmp);

  return Big;
}
/*   ===========                                         */
void PrintHbonds (Hbond_Data *Hbond, int n)
/*   ===========                                         */
{
  int
    b = 2*n,
    count,
    xx = 0,
    num = 0;

  char
    nada = ' ',
    *nothing = "";

  Hbond_Data *Big;

  Big = Hbond;
  
  if( xx == 0 )  {
    puts("  Donor              Acc.    or    Don.     Acceptor        Dist  CNH_O   NH_OC");
    puts("  --------------     ------------------     --------------  ----  ------  -----");
  }     /*  Donor                                     Acceptor        Dist  CNH_O   NH_OC */

  for(count=0; count<b; count++)  {
    if(Big[count].type == 'a')  {
     
      printf("#                        %4d %1c %3s %3s ... %3s %4d %1c %3s  %4.2lf  %6.2lf  %5.1lf\n",
	     Big[count].mid->resnum, Big[count].mid->chain, Big[count].mid->resname, Big[count].mid->name,
	     Big[count].acc->name, Big[count].acc->resnum, Big[count].acc->chain, Big[count].acc->resname,
	     Big[count].Dist, Big[count].CNH_O_Angle, Big[count].NH_OC_Angle);
    }

    else if(Big[count].type == 'd')  {
      
      printf("# %4d %1c %3s %3s ... %3s %4d %1c %3s                         %4.2lf  %6.2lf  %5.1lf\n",
	     Big[count].don->resnum, Big[count].don->chain, Big[count].don->resname, Big[count].don->name,
	     Big[count].mid->name, Big[count].mid->resnum, Big[count].mid->chain, Big[count].mid->resname,
	     Big[count].Dist, Big[count].CNH_O_Angle, Big[count].NH_OC_Angle);
    }
  }
  xx++;
}


/*   ============                                                      */
void PrintHbonded
/*   ============                                                      */
(atom_record *acc, atom_record *don,
 double Dist, double CNH_O_Angle, double NH_OC_Angle)
{
  static int x = 0;

  if( x == 0 )  {
    puts("   Donor              Acceptor        Dist  CNH_O   NH_OC");
    puts("   --------------     --------------  ----  ------  -----");
  }

  printf("#  %4d %1c %3s %3s ... %3s %4d %1c %3s  %4.2lf  %6.2lf  %5.1lf\n",
	 don->resnum, don->chain, don->resname, don->name,
	 acc->name, acc->resnum, acc->chain, acc->resname,
	 Dist, CNH_O_Angle, NH_OC_Angle);

  x++;
}



/*   =========================                                         */
void FindAllHbondsInAStructure (file_records * f, bool flag, bool sort)
/*   =========================                                         */
{
# define CUTOFF_DIST 4.0

  int iac, n=0;

  AtomList *ptr;

  Hbond_Data *Hbonds;

  double
    CNH_O_Angle,
    NH_OC_Angle;

  for (iac=0; iac<f->atomnum; iac++)  {
    atom_record * acc = f->atoms[iac];
    
    /*  Tests if an acceptor on the main chain only, no HETATOMs put through */
    if ( IsAnAcceptor(acc) && strcmp(acc->ident, "HETATM"))  {
      for (ptr=acc->Contacts; ptr->Dist < CUTOFF_DIST ; ptr = ptr->NextAtom) {
		
	if (AreHbonded(acc, ptr->atom, ptr->Dist,
		       &CNH_O_Angle, &NH_OC_Angle, flag ))  {
	  if (sort)
	    n = SaveHbond(acc, ptr->atom, ptr->Dist, CNH_O_Angle,
			  NH_OC_Angle, n, &Hbonds);
	  else
	    PrintHbonded(acc, ptr->atom, ptr->Dist, CNH_O_Angle, NH_OC_Angle);
	}
      }                                                
    }
  }

  if (sort) {
    Hbonds = SortHbonds(Hbonds, n);
    PrintHbonds(Hbonds, n);
  }

# undef CUTOFF_DIST
}