#include "WaterHbonds.h"

double GLOBAL_lower_thresh, GLOBAL_upper_thresh ;
double GLOBAL_Summary, GLOBAL_Printout;

double AtomDistSq(atom_record * a1, atom_record * a2)
{
  double dx = a1->x - a2->x ;
  double dy = a1->y - a2->y ;
  double dz = a1->z - a2->z ;
  return (dx*dx + dy*dy + dz*dz);
}

void PrintWatHbonded
(atom_record *acc, atom_record *don, double Dist)
{
  printf("#  %4d %1c %3s %3s ... %3s %4d %1c %3s  %4.2lf\n",
	 don->resnum, don->chain, don->resname, don->name,
	 acc->name, acc->resnum, acc->chain, acc->resname,
	 Dist);
}

#define ISWAT(a) ((STREQ(a->resname,"HOH") || STREQ(a->resname,"DOD")) &&\
		  (a->name[0] == 'O'))


void CalcWatHbonds (file_records *f) 
{
  int ii,jj ;
  int nwat = 0;
  fprintf(stderr,"%3s %3s %2s %2s %2s %4s %4s %4s\n",
	 "nwat","nres","ww","wp","wwwp","wwd","wpd","wwwpd");
  for (ii=0; ii<f->atomnum; ii++) {
    atom_record * a = f->atoms[ii];
    if (ISWAT(a) && a->chain == ' ') {
      int wwwp = 0;
      int ww = 0;
      int wp = 0;
      double wwwpd = 0;
      double wwd = 0;
      double wpd = 0;
      nwat++;

      for (jj=0; jj<f->atomnum; jj++) {
	atom_record * b = f->atoms[jj];
	if ((b->name[0] == 'N' || b->name[0] == 'O') && a != b) {
	  double d = AtomDistSq(a,b);
	  if (GLOBAL_lower_thresh < d && d < GLOBAL_upper_thresh) {
	    if (GLOBAL_Printout) {
	      PrintWatHbonded(b,a,sqrt(d));
	    }
	    wwwp++; wwwpd += (b->occ + a->occ)/2.0 ;
	    if (ISWAT(b)) {
	      ww++; wwd += (b->occ + a->occ)/2.0 ;
	    } else {
	      wp++; wpd += (b->occ + a->occ)/2.0 ;

	    }
	  }
	}
      }
      if (GLOBAL_Summary) {
	printf("%3d %3d %2d %2d %2d %4.2f %4.2f %4.2f\n",nwat,
	       a->resnum,ww,wp,wwwp,wwd,wpd,wwwpd);
      }
    }
  }
}

void VDWCont1to2(file_records *fa, file_records * fb, double thresh) 
{
  int ii,jj ;
  fprintf(stderr,"Structure A has %d atoms and structure B has %d atoms\n",
	  fa->atomnum,fb->atomnum);
  for (ii=0; ii<fa->atomnum; ii++) {
    atom_record * a = fa->atoms[ii];
    for (jj=0; jj<fb->atomnum; jj++) {
      atom_record * b = fb->atoms[jj];
      double d = sqrt(AtomDistSq(a,b));
      if (d - (b->definition->vdw + a->definition->vdw) < thresh) {
	write_pdb_record(stdout,a,0);
	break ;
      }
    }
  }
}


void PairCloserThanThresh(file_records *fa, file_records * fb, double thresh) 
{
  int ii,jj ;
  fprintf(stderr,"Structure A has %d atoms and structure B has %d atoms\n",
	  fa->atomnum,fb->atomnum);
  for (ii=0; ii<fa->atomnum; ii++) {
    atom_record * a = fa->atoms[ii];
    for (jj=0; jj<fb->atomnum; jj++) {
      atom_record * b = fb->atoms[jj];
      double d = sqrt(AtomDistSq(a,b));
      if (d <= thresh) {
	char bufa[BUFSIZ],bufb[BUFSIZ];
	sprintfAtomShortForm (bufa,a);
	sprintfAtomShortForm (bufb,b);
	printf("Distance-is  %8.5f  from  %10s  to  %10s\n",d,bufa,bufb);
	break ;
      }
    }
  }
}




void Hbonds1to2 (file_records *fa, file_records *fb) 
{
  int ii,jj ;
  int nwat = 0;
  fprintf(stderr,"%3s %3s %2s %2s %2s %4s %4s %4s\n",
	 "nwat","nres","ww","wp","wwwp","wwd","wpd","wwwpd");
  for (ii=0; ii<fa->atomnum; ii++) {
    atom_record * a = fa->atoms[ii];
    int wwwp = 0;
    int ww = 0;
    int wp = 0;
    double wwwpd = 0;
    double wwd = 0;
    double wpd = 0;
    nwat++;
    for (jj=0; jj<fb->atomnum; jj++) {
      atom_record * b = fb->atoms[jj];
      if ((b->name[0] == 'N' || b->name[0] == 'O') && a != b) {
	double d = AtomDistSq(a,b);
	if (GLOBAL_lower_thresh < d && d < GLOBAL_upper_thresh) {
	  if (GLOBAL_Printout) {
	    PrintWatHbonded(b,a,sqrt(d));
	  }
	  wwwp++; wwwpd += (b->occ + a->occ)/2.0 ;
	  if (ISWAT(b)) {
	    ww++; wwd += (b->occ + a->occ)/2.0 ;
	  } else {
	    wp++; wpd += (b->occ + a->occ)/2.0 ;
	    
	  }
	}
      }
    }
    if (GLOBAL_Summary) {
      printf("%3d %3d %2d %2d %2d %4.2f %4.2f %4.2f\n",nwat,
	     a->resnum,ww,wp,wwwp,wwd,wpd,wwwpd);
    }
  }
}