/***************************************************************************
                          pdbstat_mod.cpp  -  description
                             -------------------
    begin                : Wed Sep 28 18:42:55 2005
    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 <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <protein.h>
#include <nucleicacid.h>
#include <lipid.h>
#include <solvent.h>
#include <molecules.h>

using namespace Almost;

#define ALMTOSTRING(s) \
template<> \
inline string to_string(const s & ){ \
  return "<" #s ">"; \
}

ALMTOSTRING(Molecules);

void count_neighbors_cog(const Molecules & molecules,  int sk,
			 string file){
  vector<vector<int> > fa;
  fa.resize(molecules.fragment_size());
  for(int i=0;i<fa.size();i++){
    fa[i] = molecules.fragment_atoms(i);
  }
  //compute cog;
  vector<double> X,Y,Z;
  for(int i=0;i<fa.size();i++){
    double x = 0;
    double y = 0;
    double z = 0;
    for(int j=0;j<fa[i].size();j++){
      x += molecules[fa[i][j]][0];
      y += molecules[fa[i][j]][1];
      z += molecules[fa[i][j]][2];
    }
    X.push_back(x/fa[i].size());
    Y.push_back(y/fa[i].size());
    Z.push_back(z/fa[i].size());
  }

  vector<int> count0;//less 7A
  vector<int> count1;// 7-10
  vector<int> count2;//10-12
  vector<int> count3;//12+

  count0.resize(X.size());
  count1.resize(X.size());
  count2.resize(X.size());
  count3.resize(X.size());

  for(int i=0;i<X.size();i++){
    for(int j=i+sk;j<X.size();j++){
      double dx = X[i]-X[j];
      double dy = Y[i]-Y[j];
      double dz = Z[i]-Z[j];
      double d = sqrt(dx*dx+dy*dy+dz*dz);
      if(d<7){
	++count0[i];
	++count0[j];
	continue;
      }
      if((d<10)&&(d>=7)){
	++count1[i];
	++count1[j];
	continue;
      }
      if((d<12)&&(d>=10)){
	++count2[i];
	++count2[j];
	continue;
      }
      if((d>=12)){
	++count3[i];
	++count3[j];
	continue;
      }
    }
  }
  ofstream out;
  out.open(file.c_str());
  for(int i=0;i<X.size();i++){
    string n = molecules.fragment_name(i,Molecules::LONG);
    out<<n<<"\t"
       <<count0[i]<<" "
       <<count1[i]<<" "
       <<count2[i]<<" "
       <<count3[i]<<"\n";
  }  
}


void count_neighbors_pair(const Molecules & molecules,  int sk,
			  string file){
  vector<vector<int> > fa;
  fa.resize(molecules.fragment_size());
  for(int i=0;i<fa.size();i++){
    fa[i] = molecules.fragment_atoms(i);
  }
  //compute cog;
  vector<double> X,Y,Z;
  for(int i=0;i<fa.size();i++){
    double x = 0;
    double y = 0;
    double z = 0;
    for(int j=0;j<fa[i].size();j++){
      x += molecules[fa[i][j]][0];
      y += molecules[fa[i][j]][1];
      z += molecules[fa[i][j]][2];
    }
    X.push_back(x/fa[i].size());
    Y.push_back(y/fa[i].size());
    Z.push_back(z/fa[i].size());
  }

  vector<int> count0;//less 7A
  vector<int> count1;// 7-10
  vector<int> count2;//10-12
  vector<int> count3;//12+

  count0.resize(X.size());
  count1.resize(X.size());
  count2.resize(X.size());
  count3.resize(X.size());

  for(int i=0;i<X.size();i++){
    for(int j=i+sk;j<X.size();j++){
      double dx = X[i]-X[j];
      double dy = Y[i]-Y[j];
      double dz = Z[i]-Z[j];
      double d = sqrt(dx*dx+dy*dy+dz*dz);
      if(d<7){
	++count0[i];
	++count0[j];
	continue;
      }
      if((d<10)&&(d>=7)){
	++count1[i];
	++count1[j];
	continue;
      }
      if((d<12)&&(d>=10)){
	++count2[i];
	++count2[j];
	continue;
      }
      if((d>=12)){
	++count3[i];
	++count3[j];
	continue;
      }
    }
  }
  ofstream out;
  out.open(file.c_str());
  for(int i=0;i<X.size();i++){
    for(int j=i+sk;j<X.size();j++){
      double dx = X[i]-X[j];
      double dy = Y[i]-Y[j];
      double dz = Z[i]-Z[j];
      double d = sqrt(dx*dx+dy*dy+dz*dz);
      string n1 = molecules.fragment_name(i,Molecules::LONG);
      string n2 = molecules.fragment_name(j,Molecules::LONG);
      out<<n1<<" "<<i+1<<" "<<n2<<" "<<" "<<j+1<<" "<<d<<" "<<"\t"
	 <<count0[i]<<" "
	 <<count1[i]<<" "
	 <<count2[i]<<" "
	 <<count3[i]<<" "
	 <<count0[j]<<" "
	 <<count1[j]<<" "
	 <<count2[j]<<" "
	 <<count3[j]<<"\n";
    }  
  }
}


void count_neighbors(const Molecules & molecules, string atm, int sk,
		     string file){
  vector<int> A;
  vector<int> count0;//less 7A
  vector<int> count1;// 7-10
  vector<int> count2;//10-12
  vector<int> count3;//12+
  int size = molecules.atom_size();
  for(int i=0;i<size;i++)
    if(molecules[i].name()==atm){
      A.push_back(i);
      //      cout<<"DEB "<<atm<<" "<<molecules[i].name()<<endl;
    }

  //  cout<<"DEB "<<A.size()<<" atoms of type "<<atm<<" detected\n";
  count0.resize(A.size());
  count1.resize(A.size());
  count2.resize(A.size());
  count3.resize(A.size());
  
  for(int i=0;i<A.size();i++){
    for(int j=i+sk;j<A.size();j++){
      double d = molecules.distance(A[i],A[j]);
      //       {
      // 	string n1 = molecules.atom_name(A[i],Molecules::LONG);
      // 	string n2 = molecules.atom_name(A[j],Molecules::LONG);
      // 	cout<<"DEB "<<n1<<" "<<n2<<" "<<d<<"\n";
      //       }
      if(d<7){
	++count0[i];
	++count0[j];
	continue;
      }
      if((d<10)&&(d>=7)){
	++count1[i];
	++count1[j];
	continue;
      }
      if((d<12)&&(d>=10)){
	++count2[i];
	++count2[j];
	continue;
      }
      if((d>=12)){
	++count3[i];
	++count3[j];
	continue;
      }
    }
  }
  ofstream out;
  out.open(file.c_str());
  for(int i=0;i<A.size();i++){
    string n = molecules.atom_name(A[i],Molecules::LONG);
    out<<n<<"\t"
       <<count0[i]<<" "
       <<count1[i]<<" "
       <<count2[i]<<" "
       <<count3[i]<<"\n";
  }
}


class SecStructGeometry {
  string ss; 
  vector<pair<int,int> > Beta;
  vector<pair<int,int> > Helix;
  vector<int> N,CA,C,O;
  const Molecules * molecules_;
public:
  SecStructGeometry(const Molecules & molecules,
		    string ss_file){
    molecules_ = & molecules;
    //Read ss_file
    ifstream in; in.open(ss_file.c_str());
    if(!in){
      cerr<<"No such file: "<<ss_file<<"\n";
      throw "IO error";
    }
    in>>ss;
    in.close();
    cout<<ss<<"\n";
    //
    for(int i=0;i<ss.size();i++){
      if((ss[i]=='B')||(ss[i]=='E')||(ss[i]=='b')) ss[i] = 'E';
      if((ss[i]=='H')||(ss[i]=='G')||(ss[i]=='I')) ss[i] = 'H';
      if((ss[i]=='C')||(ss[i]=='T')) ss[i] = 'C';
    }
    cout<<ss<<"\n";
    find_E();
    find_H();      

    for(int i=0;i<molecules.atom_size();i++){
      if(molecules[i].name()=="N") N.push_back(i);
      if(molecules[i].name()=="C") C.push_back(i);
      if(molecules[i].name()=="CA") CA.push_back(i);
      if(molecules[i].name()=="O") O.push_back(i);
    }
  }


  void SS_stat(string f){
    ofstream out;
    out.open(f.c_str());
    const Molecules & molecules = * molecules_;
    //hb
    for(int ss1 = 0;ss1<Beta.size();ss1++){
      for(int ss2 = ss1+1 ;ss2<Beta.size();ss2++){
	int loop = Beta[ss2].first - Beta[ss1].second;
	for(int d1=Beta[ss1].first;d1<Beta[ss1].second-1;d1++){
	  double v1[3],co1[3],co2[3];
	  v1[0] = molecules[C[d1+1]][0]-molecules[N[d1]][0];
	  v1[1] = molecules[C[d1+1]][1]-molecules[N[d1]][1];
	  v1[2] = molecules[C[d1+1]][2]-molecules[N[d1]][2];

	  co1[0] = molecules[O[d1]][0] - molecules[C[d1]][0];
	  co1[1] = molecules[O[d1]][1] - molecules[C[d1]][1];
	  co1[2] = molecules[O[d1]][2] - molecules[C[d1]][2];
	  norm(co1);
	  //	  cout<<"Dotp "<<d1<<" "<<dp1<<"\n";
	  co2[0] = molecules[O[d1+1]][0] - molecules[C[d1+1]][0];
	  co2[1] = molecules[O[d1+1]][1] - molecules[C[d1+1]][1];
	  co2[2] = molecules[O[d1+1]][2] - molecules[C[d1+1]][2];
	  norm(co2);
	  //	  cout<<"Dotp "<<d1+1<<" "<<dp2<<"\n";
	  for(int d2=Beta[ss2].first;d2<Beta[ss2].second-1;d2++){
	    double v2[3],r[3],co3[3],co4[3];
	    v2[0] = molecules[C[d2+1]][0]-molecules[N[d2]][0];
	    v2[1] = molecules[C[d2+1]][1]-molecules[N[d2]][1];
	    v2[2] = molecules[C[d2+1]][2]-molecules[N[d2]][2];

	    r[0] = (molecules[N[d2]][0]+v2[0]/2)-
	      (molecules[N[d1]][0]+v1[0]/2);
	    r[1] = (molecules[N[d2]][1]+v2[1]/2)-
	      (molecules[N[d1]][1]+v1[1]/2);
	    r[2] = (molecules[N[d2]][2]+v2[2]/2)-
	      (molecules[N[d1]][2]+v1[2]/2);
	    double dist = norm(r);
	    co3[0] = molecules[O[d2]][0] - molecules[C[d2]][0];
	    co3[1] = molecules[O[d2]][1] - molecules[C[d2]][1];
	    co3[2] = molecules[O[d2]][2] - molecules[C[d2]][2];
	    norm(co3);

	    co4[0] = molecules[O[d2+1]][0] - molecules[C[d2+1]][0];
	    co4[1] = molecules[O[d2+1]][1] - molecules[C[d2+1]][1];
	    co4[2] = molecules[O[d2+1]][2] - molecules[C[d2+1]][2];
	    norm(co4);


	    double dp1 = dotp(r,co1);
	    double dp2 = dotp(r,co2);
	    double dp3 = dotp(r,co3);
	    double dp4 = dotp(r,co4);
	    norm(v1);norm(v2);
	    double theta = acos(dotp(v1,v2));

	    //phi
	    double t1[3],t2[3];

	    vectorProduct(v1,r,t1);
	    vectorProduct(v1,t1,t2);
	    double q1 = dotp(v2,t1);
	    double q2 = dotp(v2,t2);
	    double phi = atan2(q1,q2);
	    
	    if(dist<10){	     
	      out<<loop<<" "
		 <<dist<<" "
		 <<theta*(180/M_PI)<<" "
		 <<phi*(180/M_PI) <<" "
		 <<0.5*(abs(dp1)+abs(dp2)+abs(dp3)+abs(dp4))<<" "
		 <<dp1+dp2<<" "<<dp3+dp4<<"\n";
	    }
	  }
	}
      }
    }
  }

  void HH_stat(string f){
    ofstream out;
    out.open(f.c_str());
    const Molecules & molecules = * molecules_;
    for(int hh1 = 0; hh1<Helix.size();hh1<hh1++){
      for(int hh2 = hh1+1;hh2<Helix.size();hh2++){
	int loop = Helix[hh2].first - Helix[hh1].second;
	for(int d1 = Helix[hh1].first;d1<Helix[hh1].second-3;d1++){
	  //find v1
	  double v1[3],m1[3];
	  helix_axe(d1,v1,m1);//begin 4 HHHH

	  for(int d2 =  Helix[hh2].first;d2<Helix[hh2].second-3;d2++){
	    double v2[3],m2[3],r[3];
	    helix_axe(d2,v2,m2);
	    for(int i=0;i<3;i++) r[i] = m2[i]-m1[i];
	    double dist = norm(r);
	    double theta = acos(dotp(v1,v2));
	    double t1[3],t2[3];
	    vectorProduct(v1,r,t1);
	    vectorProduct(v1,t1,t2);
	    double q1 = dotp(v2,t1);
	    double q2 = dotp(v2,t2);
	    double phi = atan2(q1,q2);
	    if(dist<14){	     
	      out<<loop<<" "
		 <<dist<<" "
		 <<theta*(180/M_PI)<<" "
		 <<phi*(180/M_PI) <<"\n";
	    }
	  }
	}
      }
    }
  }

  void HS_stat(string f){
    ofstream out;
    out.open(f.c_str());
    const Molecules & molecules = * molecules_;
    for(int hh1 = 0; hh1<Helix.size();hh1<hh1++){
      for(int ss1 = 0; ss1<Beta.size();ss1++){
	int loop;
	int hs =0;
	if(Helix[hh1].first<Beta[ss1].first){
	  hs =1;
	  loop = Beta[ss1].first - Helix[hh1].second;
	}
	else {
	  loop = Helix[hh1].first - Beta[ss1].second;
	}

	for(int d1 = Helix[hh1].first;d1<Helix[hh1].second-3;d1++){
	  //find v1
	  double v1[3],m1[3];
	  helix_axe(d1,v1,m1);//begin 4 HHHH
	  for(int d2 = Beta[ss1].first;d2<Beta[ss1].second-1;d2++){
	    double v2[3],m2[3],r[3];
	    for(int i=0;i<3;i++){
	      v2[i] = molecules[C[d2+1]][i]-molecules[N[d2]][i];
	      m2[i] = molecules[N[d2]][i]+0.5*v2[i];
	    }
	    norm(v2);
	    for(int i=0;i<3;i++) r[i] = m2[i]-m1[i];
	    double dist = norm(r);
	    double theta = acos(dotp(v1,v2));
	    double t1[3],t2[3];
	    vectorProduct(v1,r,t1);
	    vectorProduct(v1,t1,t2);
	    double q1 = dotp(v2,t1);
	    double q2 = dotp(v2,t2);
	    double phi = atan2(q1,q2);
	    if(dist<14){	     
	      out<<loop<<" "
		 <<dist<<" "
		 <<theta*(180/M_PI)<<" "
		 <<phi*(180/M_PI) <<"\n";
	    }
	  }
	}
      }
    }
  }
private:
  void helix_axe(int b,double v[3],double m1[3]){
    double t1[3]= {0,0,0}, t2[3] = {0,0,0};
    const Molecules & molecules = * molecules_;
    for(int i=0;i<3;i++){
      t1[i] = 
	molecules[N[b]][i]+ molecules[CA[b]][i] + molecules[C[b]][i] +
	molecules[N[b+1]][i]+ molecules[CA[b+1]][i] + molecules[C[b+1]][i] +
	molecules[N[b+2]][i]+ molecules[CA[b+2]][i] + molecules[C[b+2]][i] +
	molecules[N[b+3]][i]+ molecules[CA[b+3]][i];

      t2[i] = 
	                      molecules[CA[b]][i] + molecules[C[b]][i] +
	molecules[N[b+1]][i]+ molecules[CA[b+1]][i] + molecules[C[b+1]][i] +
	molecules[N[b+2]][i]+ molecules[CA[b+2]][i] + molecules[C[b+2]][i] +
	molecules[N[b+3]][i]+ molecules[CA[b+3]][i] + molecules[C[b+3]][i];

      t1[i] = t1[i]/11.0;
      t2[i] = t2[i]/11.0;
      v[i] = (t2[i]-t1[i]);
      
      m1[i] = t1[i]+0.5*v[i];
      norm(v);
    }
  }

  void find_E(){
    int beta = -1;
    for(int i=0;i<ss.size();i++){
      if(ss[i]=='E'){ if(beta==-1) beta = i;}
      else{
	if(beta>=0){
	  if(i-beta>1)
	    Beta.push_back(pair<int,int>(beta,i));
	}
	beta = -1;
      }
    }
    for(int i=0;i<Beta.size();i++){
      cout<<"Strand from: "<<Beta[i].first+1<<" "<<Beta[i].second<<endl;
    }
  }

  void find_H(){
    int helix = -1;
    for(int i=0;i<ss.size();i++){
      if(ss[i]=='H'){ if(helix==-1) helix = i;}
      else{
	if(helix>=0){
	  if(i-helix>3)
	    Helix.push_back(pair<int,int>(helix,i));
	}
	helix = -1;
      }
    }
    for(int i=0;i<Helix.size();i++){
      cout<<"Helix from: "<<Helix[i].first+1<<" "<<Helix[i].second<<endl;
    }
  }

  double norm(double v[3]){
    double d = 0;
    for(int i=0;i<3;i++) d += v[i]*v[i];
    d = sqrt(d);
    for(int i=0;i<3;i++) v[i] =v[i]/d;
    return d;
  }


  double dotp(double v1[3],double v2[3]){
    double d = 0;
    for(int i=0;i<3;i++)
      d +=v1[i]*v2[i];
    return d;
  }

  inline void vectorProduct(const double * a, const double * b, 
			    double * ab)
  {      
    ab[0] = a[1]*b[2]- a[2]*b[1];
    ab[1] = a[2]*b[0] - a[0]*b[2];
    ab[2] = a[0]*b[1] - a[1]*b[0]; 
  }

};

double norm(double v[3]){
  double d = 0;
  for(int i=0;i<3;i++) d += v[i]*v[i];
  d = sqrt(d);
  for(int i=0;i<3;i++) v[i] =v[i]/d;
  return d;
}


double dotp(double v1[3],double v2[3]){
  double d = 0;
  for(int i=0;i<3;i++)
    d +=v1[i]*v2[i];
  return d;
}

inline void vectorProduct(const double * a, const double * b, 
			  double * ab)
{      
  ab[0] = a[1]*b[2]- a[2]*b[1];
  ab[1] = a[2]*b[0] - a[0]*b[2];
  ab[2] = a[0]*b[1] - a[1]*b[0]; 
}


//double dotscore(

// void sectruct_geometry(const Molecules & molecules,
// 		       string ss_file, string file){
//   //Read ss_file
//   ifstream in; in.open(ss_file.c_str());
//   if(!in){
//     cerr<<"No such file: "<<ss_file<<"\n";
//     throw "IO error";
//   }

//   string ss; in>>ss;
//   in.close();
//   cout<<ss<<"\n";
//   //
//   for(int i=0;i<ss.size();i++){
//     if((ss[i]=='B')||(ss[i]=='E')||(ss[i]=='b')) ss[i] = 'E';
//     if((ss[i]=='H')||(ss[i]=='G')||(ss[i]=='I')) ss[i] = 'H';
//     if((ss[i]=='C')||(ss[i]=='T')) ss[i] = 'C';
//   }
//   cout<<ss<<"\n";
//   //find E
//   vector<pair<int,int> > Beta;
//   vector<pair<int,int> > Helix;
//   int beta = -1;
//   for(int i=0;i<ss.size();i++){
//     if(ss[i]=='E'){ if(beta==-1) beta = i;}
//     else{
//       if(beta>=0){
// 	if(i-beta>1)
// 	  Beta.push_back(pair<int,int>(beta,i));
//       }
//       beta = -1;
//     }
//   }
//   for(int i=0;i<Beta.size();i++){
//     cout<<"Strand from: "<<Beta[i].first+1<<" "<<Beta[i].second<<endl;
//   }

//   int helix = -1;
//   for(int i=0;i<ss.size();i++){
//     if(ss[i]=='H'){ if(helix==-1) helix = i;}
//     else{
//       if(helix>=0){
// 	if(helix-i>3)
// 	  Helix.push_back(pair<int,int>(helix,i));
//       }
//       helix = -1;
//     }
//   }
//   for(int i=0;i<Helix.size();i++){
//     cout<<"Helix from: "<<Helix[i].first+1<<" "<<Helix[i].second<<endl;
//   }

//   //Analize Beta
//   {
//     vector<int> N,C,O;
//     for(int i=0;i<molecules.atom_size();i++){
//       if(molecules[i].name()=="N") N.push_back(i);
//       if(molecules[i].name()=="C") C.push_back(i);
//       if(molecules[i].name()=="O") O.push_back(i);
//     }

//     //hb
//     for(int ss1 = 0;ss1<Beta.size();ss1++){
//       for(int ss2 = ss1+1 ;ss2<Beta.size();ss2++){
// 	cout<<"Analizing strands "<<ss1<<" "<<ss2<<"\n";
// 	for(int d1=Beta[ss1].first;d1<Beta[ss1].second-1;d1++){
// 	  double v1[3],co1[3],co2[3];
// 	  v1[0] = molecules[C[d1+1]][0]-molecules[N[d1]][0];
// 	  v1[1] = molecules[C[d1+1]][1]-molecules[N[d1]][1];
// 	  v1[2] = molecules[C[d1+1]][2]-molecules[N[d1]][2];

// 	  co1[0] = molecules[O[d1]][0] - molecules[C[d1]][0];
// 	  co1[1] = molecules[O[d1]][1] - molecules[C[d1]][1];
// 	  co1[2] = molecules[O[d1]][2] - molecules[C[d1]][2];
// 	  norm(co1);
// 	  //	  cout<<"Dotp "<<d1<<" "<<dp1<<"\n";
// 	  co2[0] = molecules[O[d1+1]][0] - molecules[C[d1+1]][0];
// 	  co2[1] = molecules[O[d1+1]][1] - molecules[C[d1+1]][1];
// 	  co2[2] = molecules[O[d1+1]][2] - molecules[C[d1+1]][2];
// 	  norm(co2);
// 	  //	  cout<<"Dotp "<<d1+1<<" "<<dp2<<"\n";
// 	  for(int d2=Beta[ss2].first;d2<Beta[ss2].second-1;d2++){
// 	    double v2[3],r[3],co3[3],co4[3];
// 	    v2[0] = molecules[C[d2+1]][0]-molecules[N[d2]][0];
// 	    v2[1] = molecules[C[d2+1]][1]-molecules[N[d2]][1];
// 	    v2[2] = molecules[C[d2+1]][2]-molecules[N[d2]][2];

// 	    r[0] = (molecules[N[d2]][0]+v2[0]/2)-
// 	      (molecules[N[d1]][0]+v1[0]/2);
// 	    r[1] = (molecules[N[d2]][1]+v2[1]/2)-
// 	      (molecules[N[d1]][1]+v1[1]/2);
// 	    r[2] = (molecules[N[d2]][2]+v2[2]/2)-
// 	      (molecules[N[d1]][2]+v1[2]/2);
// 	    double dist = norm(r);
// 	    co3[0] = molecules[O[d2]][0] - molecules[C[d2]][0];
// 	    co3[1] = molecules[O[d2]][1] - molecules[C[d2]][1];
// 	    co3[2] = molecules[O[d2]][2] - molecules[C[d2]][2];
// 	    norm(co3);

// 	    co4[0] = molecules[O[d2+1]][0] - molecules[C[d2+1]][0];
// 	    co4[1] = molecules[O[d2+1]][1] - molecules[C[d2+1]][1];
// 	    co4[2] = molecules[O[d2+1]][2] - molecules[C[d2+1]][2];
// 	    norm(co4);


// 	    double dp1 = dotp(r,co1);
// 	    double dp2 = dotp(r,co2);
// 	    double dp3 = dotp(r,co3);
// 	    double dp4 = dotp(r,co4);
// 	    norm(v1);norm(v2);
// 	    double theta = acos(dotp(v1,v2));

// 	    //phi
// 	    double t1[3],t2[3];

// 	    vectorProduct(v1,r,t1);
// 	    vectorProduct(v1,t1,t2);
// 	    double q1 = dotp(v2,t1);
// 	    double q2 = dotp(v2,t2);
// 	    double phi = atan2(q1,q2);
// 	    if(dist<10){	     
// 	      cout<<"Dimer distance: "<<d1<<" "<<d2<<" "<< dist<<"\n";
// 	      cout<<"\t> Theta "<<theta*(180/M_PI)<<"\n";
// 	      cout<<"\t> Phi   "<<phi*(180/M_PI)<<"\n";
// 	      cout<<"\t> Dotp "<<0.5*(abs(dp1)+abs(dp2)+abs(dp3)+abs(dp4))<<" "
// 		  <<dp1+dp2<<" "<<dp3+dp4<<"\n";
// 	    }
// 	  }
// 	}
//       }
//     }
//   }
//}
inline string to_string(const SecStructGeometry &){
  return "<SecStructGeometry>";
}

extern "C" {
  void init_pdbstat(){
    //declarations here

    Module mod = Module("pdbstat","Compute statistics on sets of PDBs");
    mod
      .def_function("count_neighbors","Count neighbors of specified type",count_neighbors)
      .def_function("count_neighbors_cog","Count neighbors of specified type",count_neighbors_cog)
      .def_function("count_neighbors_pair","Count neighbors of specified type",count_neighbors_pair)
      ;

    Class<SecStructGeometry,ArgList2<const Molecules &, 
      string> >(mod.self(),"ss_geomerty")
      .def_method("SS_stat",&SecStructGeometry::SS_stat)
      .def_method("HH_stat",&SecStructGeometry::HH_stat)
      .def_method("HS_stat",&SecStructGeometry::HS_stat)
      ;

  }
}