/***************************************************************************
                          dcd.cpp  -  description
                             -------------------
    begin                : Thu Jul  3 00:58:54 2003
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli@bioc.unizh.ch
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <dcd.h>
#include <fstream>

namespace Almost {
  void rgyr(DCDTrj & dcd, string file){
    ofstream out;
    out.open(file.c_str());
    int natoms = dcd.get_natoms();
 
    for(int i=0;i<dcd.get_nframes();i++){
      float rgyr=0;
      float centr[3] = {0,0,0};
    
      int xnatoms = Coor<float>::DIM*natoms;
      for(int j=0;j<xnatoms;j = j+Coor<float>::DIM){
	centr[0] +=dcd.coor(i).coor[j];
	centr[1] +=dcd.coor(i).coor[j+1];
	centr[2] +=dcd.coor(i).coor[j+2];
      }


      for(int j=0;j<3;j++) centr[j] /= natoms;

      for(int j=0;j<xnatoms;j = j+Coor<float>::DIM){
	rgyr += (dcd.coor(i).coor[j]-centr[0])*(dcd.coor(i).coor[j]-centr[0]);
	rgyr += (dcd.coor(i).coor[j+1]-centr[1])*(dcd.coor(i).coor[j+1]-centr[1]);
	rgyr += (dcd.coor(i).coor[j+2]-centr[2])*(dcd.coor(i).coor[j+2]-centr[2]);
      }

      rgyr /= natoms;
    
      out<<dcd.get_step()*i<<" "<< sqrt(rgyr)<<"\n";
    }
  }



  void hbond_finder(NameSelection & solute, NameSelection & solvent,
		    const HBondFinderOpts& opts,
		    DCDTrj & dcd,
		    string outfile){
    //Tests 
    //	if(solute.molecules()!=solvent.molecules()){
    //		throw "Selections refer to different molecules";
    //	}
    double X = opts.X;
    double Y = opts.Y;
    double Z = opts.Z;
    double d = opts.d;
    double ang = opts.angle;	  
    
    const Molecules & mols = solute.molecules();
    //Setup protein   
    vector<int> N;
    vector<int> H; //ptorein
    vector<int> O;
    vector<int> O2; //water
    vector<int> H1, H2;
	  
    //Protein
    for(int i=0;i<mols.protein_size();i++){
      const Protein & prot = mols.protein(i);
      int b = prot.fragment_offset();
      int e = prot.fragment_size()+b;
		
      for(int j = b; j<e;j++){
	vector<int> atm = prot.fragment_atoms(j);
	int n = -1;
	int o = -1;
	int h = -1;
	for(int aa = 0; aa< atm.size();aa++){
	  if(!solute.is_selected(atm[aa]))continue;
	  if(prot[atm[aa]].name()=="N") n = atm[aa];
	  if(prot[atm[aa]].name()=="O") o = atm[aa];
	  if((prot[atm[aa]].name()=="H")||(prot[atm[aa]].name()=="HN"))
	    h = atm[aa];
				
	}
	if((n!=-1)&&(h!=-1)){ N.push_back(n); H.push_back(h);}
	if(o!=-1){ O.push_back(o);}
      }
    }

	
    //Waters
    for(int i=0;i<mols.solvent_size();i++){
      const Solvent & solv = mols.solvent(i);
      int b = solv.fragment_offset();
      int e = solv.fragment_size()+b;
		
      for(int j = b; j<e;j++){
	vector<int> atm = solv.fragment_atoms(j);
	int oh = -1;
	int h1 = -1;
	int h2 = -1;
	for(int aa = 0; aa< atm.size();aa++){
	  if(!solvent.is_selected(atm[aa]))continue;
	  if(solv[atm[aa]].name()=="OH2") oh = atm[aa];
	  if(solv[atm[aa]].name()=="H1") h1 = atm[aa];
	  if(solv[atm[aa]].name()=="H2") h2 = atm[aa];
	}
	O2.push_back(oh);
	H1.push_back(h1); 
	H2.push_back(h2);
      }
    }	

    cout<<"\t Number of Protein N-H: "<<N.size()<<"\n";
    cout<<"\t Number of Protein O  : "<<O.size()<<"\n";
    cout<<"\t Number of Water   O  : "<<O2.size()<<"\n";
    cout<<"\t Number of Water   H1 : "<<H1.size()<<"\n";
    cout<<"\t Number of Water   H2 : "<<H2.size()<<"\n";
    ofstream out;
    out.open(outfile.c_str());
    for(int f = 0; f< dcd.get_nframes(); f++){
      out<<"Frame "<<f+1<<"\n";
      //Search N-H O hbonds
      for(int h = 0; h< H.size();h++){
	for(int o = 0 ; o < O2.size(); o++){
	  //dist 
	  float dx = dcd.coor(f).coor[H[h]*Coor<float>::DIM] -
	    dcd.coor(f).coor[O2[o]*Coor<float>::DIM];

	  float dy = dcd.coor(f).coor[H[h]*Coor<float>::DIM+1] -
	    dcd.coor(f).coor[O2[o]*Coor<float>::DIM+1];

	  float dz = dcd.coor(f).coor[H[h]*Coor<float>::DIM+2] -
	    dcd.coor(f).coor[O2[o]*Coor<float>::DIM+2];
	  
	  dx = abs(dx);
	  dy = abs(dy);
	  dz = abs(dz);

	  dx = dx - X*floor(dx/X);
	  dy = dy - Y*floor(dy/Y);
	  dz = dz - Z*floor(dz/Z);
	  
	  if(dx>X/2) dx = X -dx;
	  if(dy>Y/2) dy = Y -dy;
	  if(dz>Z/2) dz = Z -dz;
	  float d2 = dx*dx + dy*dy +dz*dz;
	  if(d2<(d*d)){
	    //Angle 
	    float * c1, *c2, *c3;
	    c1 =& dcd.coor(f).coor[N[h]*Coor<float>::DIM];
	    c2 =& dcd.coor(f).coor[H[h]*Coor<float>::DIM];
	    c3 =& dcd.coor(f).coor[O2[o]*Coor<float>::DIM];
	    float angle = 
	      Kernel::Geometry::angle(c1,c2,c3);
	    if(angle>ang*M_PI/180){
	      //Yes
	    out<<"\t"<<mols.atom_name(N[h],Molecules::BASE)<<" "
	       <<mols.atom_name(H[h],Molecules::BASE)<<" "
	       <<mols.atom_name(O2[o],Molecules::BASE)<<" "
	       <<sqrt(d2)<<" "<<angle*180.0/M_PI<<"\n";
	    }
	  }
	}
      }

      //Search O H1 O2 hbonds
      for(int h = 0; h< O.size();h++){
	for(int o = 0 ; o < H1.size(); o++){
	  //dist 
	  float dx = dcd.coor(f).coor[O[h]*Coor<float>::DIM] -
	    dcd.coor(f).coor[H1[o]*Coor<float>::DIM];

	  float dy = dcd.coor(f).coor[O[h]*Coor<float>::DIM+1] -
	    dcd.coor(f).coor[H1[o]*Coor<float>::DIM+1];

	  float dz = dcd.coor(f).coor[O[h]*Coor<float>::DIM+2] -
	    dcd.coor(f).coor[H1[o]*Coor<float>::DIM+2];
	  
	  dx = abs(dx);
	  dy = abs(dy);
	  dz = abs(dz);

	  dx = dx - X*floor(dx/X);
	  dy = dy - Y*floor(dy/Y);
	  dz = dz - Z*floor(dz/Z);
	  
	  if(dx>X/2) dx = X -dx;
	  if(dy>Y/2) dy = Y -dy;
	  if(dz>Z/2) dz = Z -dz;
	  float d2 = dx*dx + dy*dy +dz*dz;
	  if(d2<(d*d)){
	    //Angle 
	    float * c1, *c2, *c3;
	    c1 =& dcd.coor(f).coor[O[h]*Coor<float>::DIM];
	    c2 =& dcd.coor(f).coor[H1[o]*Coor<float>::DIM];
	    c3 =& dcd.coor(f).coor[O2[o]*Coor<float>::DIM];
	    float angle = 
	      Kernel::Geometry::angle(c1,c2,c3);
	    if(angle>ang*M_PI/180){
	      //Yes
	    out<<"\t"<<mols.atom_name(O[h],Molecules::BASE)<<" "
	       <<mols.atom_name(H1[o],Molecules::BASE)<<" "
	       <<mols.atom_name(O2[o],Molecules::BASE)<<" "
	       <<sqrt(d2)<<" "<<angle*180.0/M_PI<<"\n";
	    }
	  }
	}
      }

      //Search O H2 O2 hbonds
      for(int h = 0; h< O.size();h++){
	for(int o = 0 ; o < H2.size(); o++){
	  //dist 
	  float dx = dcd.coor(f).coor[O[h]*Coor<float>::DIM] -
	    dcd.coor(f).coor[H2[o]*Coor<float>::DIM];

	  float dy = dcd.coor(f).coor[O[h]*Coor<float>::DIM+1] -
	    dcd.coor(f).coor[H2[o]*Coor<float>::DIM+1];

	  float dz = dcd.coor(f).coor[O[h]*Coor<float>::DIM+2] -
	    dcd.coor(f).coor[H2[o]*Coor<float>::DIM+2];
	  
	  dx = abs(dx);
	  dy = abs(dy);
	  dz = abs(dz);

	  dx = dx - X*floor(dx/X);
	  dy = dy - Y*floor(dy/Y);
	  dz = dz - Z*floor(dz/Z);
	  
	  if(dx>X/2) dx = X -dx;
	  if(dy>Y/2) dy = Y -dy;
	  if(dz>Z/2) dz = Z -dz;
	  float d2 = dx*dx + dy*dy +dz*dz;
	  if(d2<(d*d)){
	    //Angle 
	    float * c1, *c2, *c3;
	    c1 =& dcd.coor(f).coor[O[h]*Coor<float>::DIM];
	    c2 =& dcd.coor(f).coor[H2[o]*Coor<float>::DIM];
	    c3 =& dcd.coor(f).coor[O2[o]*Coor<float>::DIM];
	    float angle = 
	      Kernel::Geometry::angle(c1,c2,c3);
	    if(angle>ang*M_PI/180){
	      //Yes
	    out<<"\t"<<mols.atom_name(O[h],Molecules::BASE)<<" "
	       <<mols.atom_name(H2[o],Molecules::BASE)<<" "
	       <<mols.atom_name(O2[o],Molecules::BASE)<<" "
	       <<sqrt(d2)<<" "<<angle*180.0/M_PI<<"\n";
	    }
	  }
	}
      }
    }
  }


}