/***************************************************************************
                          gsl_mod.cpp  -  description
                             -------------------
    begin                : Thu Aug 26 22:37:34 2004
    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 <module.h>
#include <class.h>
#include <zextclass.h>
#include <zobj.hpp>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <protein.h>
#include <time.h>
#include <almost.h>
#include <stdlib.h>
#include <coor.h>
#include <molecules.h>
#include <moveset.h>
//#include <forcefield/hardsphere.h>
#include <coormanager/coormanager.h>
#include <forcefield/const/peconst.h>
using namespace Almost;


// class GSLVector {
//   gsl_vector * v_;
// public:
//   GSLVector(int n){
//    v_ =  gsl_vector_alloc(n);
//   }

//   GSLVector(const GSLVector & v){
//     v_ =  gsl_vector_alloc(v.size());
//     gsl_vector_memcpy(v_,v.v_);
//   }

//   ~GSLVector(){
//     gsl_vector_free(v_);
//   }

//   gsl_vector * gsl_vector_p(){ return v_;}

//   size_t size() const { return v_->size;}

//   //members
//   double get(int i) const {
//     return gsl_vector_get(v_,i);
//   }

//   void set(int i, double d){
//     gsl_vector_set(v_,i,d);
//   }
   
//   //init
//   void set_all(double d){
//     gsl_vector_set_all(v_,d);
//   }

//   void set_zero(){
//     gsl_vector_set_zero(v_);
//   }

//   void set_basis(int i){
//     gsl_vector_set_basis(v_,i);
//   }
 
//   //read -write- 
//   void print(){
//     gsl_vector_fprintf(stdout,v_,"%g");
//   }

//   void read(string file,int s){
//     if(s!=size()) return;
//     FILE *f=fopen(file.c_str(),"r");
//     gsl_vector_fscanf(f,v_);
//     fclose(f);
//   }

//   //views ??
// };

template<>
inline string to_string(const GSLVector &v){
  return "<gsl::vector>";
}



// class GSLMatrix {
//   gsl_matrix * m_;
// public:
//   GSLMatrix(int rows, int cols){
//    m_ =  gsl_matrix_alloc(rows,cols);
//   }

//   GSLMatrix(const GSLMatrix & m){
//     m_ =  gsl_matrix_alloc(m.rows(),m.cols());
//     gsl_matrix_memcpy(m_,m.m_);
//   }

//   ~GSLMatrix(){
//     gsl_matrix_free(m_);
//   }

//   gsl_matrix * gsl_matrix_p(){ return m_;}

//   int rows() const { return m_->size1;}
//   int cols() const { return m_->size2;}

//   double get(int i, int j) const {
//     return gsl_matrix_get(m_,i,j);
//   }

//  void set(int i, int j, double d){
//     gsl_matrix_set(m_,i,j,d);
//   }
   
//   //init
//   void set_all(double d){
//     gsl_matrix_set_all(m_,d);
//   }

//   void set_identity(){
//     gsl_matrix_set_identity(m_);
//   }

//   void set_zero(){
//     gsl_matrix_set_zero(m_);
//   }

//   void print(){
//     gsl_matrix_fprintf(stdout,m_,"%g");
//   }

//   void read(string file,int s1,int s2){
//     if(s1!=rows()) return;
//     if(s2!=cols()) return;
//     FILE *f=fopen(file.c_str(),"r");
//     gsl_matrix_fscanf(f,m_);
//     fclose(f);
//   }


//   void print_row(int i){
//     for(int j=0;j<m_->size1;j++){
//       cout<<gsl_matrix_get(m_,i,j)<<" ";
//     }
//     cout<<"\n";
//   }

//   void print_col(int i){
//     for(int j=0;j<m_->size1;j++){
//       cout<<gsl_matrix_get(m_,j,i)<<" ";
//     }
//     cout<<"\n";
//   }

// };

template<>
inline string to_string(const GSLMatrix &m){
  return "<gsl::matrix>";
}

void eigenv(GSLMatrix m, GSLVector & v){
  int rows  = m.rows();
  int cols  = m.cols();
  if(rows!=cols) throw "matrix is not square!";
  if(v.size()!=rows) throw "not enougth space to store EW!";
  //alloc workspace
  gsl_eigen_symm_workspace * ws =  gsl_eigen_symm_alloc(rows);
  int status = gsl_eigen_symm(m.gsl_matrix_p(),v.gsl_vector_p(),ws);
  gsl_eigen_symm_free(ws);
}

void eigenvv(GSLMatrix m,  GSLVector & v,GSLMatrix & vv){
  int rows  = m.rows();
  int cols  = m.cols();
  if(rows!=cols) throw "matrix is not square!";
  if(v.size()!=rows) throw "not enougth space to store EW!";
  //alloc workspace
  gsl_eigen_symmv_workspace * ws =  gsl_eigen_symmv_alloc(rows);
  int status = gsl_eigen_symmv(m.gsl_matrix_p(),
			       v.gsl_vector_p(),
			       vv.gsl_matrix_p(),
			       ws);
  gsl_eigen_symmv_sort( v.gsl_vector_p(),
			vv.gsl_matrix_p(),
			GSL_EIGEN_SORT_VAL_DESC);
  gsl_eigen_symmv_free(ws);
}



GSLMatrix cont_map(const Protein &p, double cut_off){
  int size = p.fragment_size();
  cout<<"\t>> Protin contat map size "<<size<<"x"<<size<<"\n";
  GSLMatrix m(size,size);
  m.set_zero();
  double cut = cut_off*cut_off;
  // fragmets
  for(int i=0;i<size;i++){
    vector<int> tmp = p.fragment_atoms(i);
    vector<int> heavy_atoms;
    for(int a = 0;a<tmp.size();a++){
      if(p.atom_name(tmp[a],Protein::SHORT)[0]!='H')
	heavy_atoms.push_back(tmp[a]);
    }
    //    cout<<"DEB heavy atoms fag "<<i<<" : "<<heavy_atoms.size()<<"\n";
    for(int j = i+3;j<size;j++){
      vector<int> tmp2 = p.fragment_atoms(j);
      vector<int> heavy_atoms2;
      for(int a = 0;a<tmp2.size();a++){
	if(p.atom_name(tmp2[a],Protein::SHORT)[0]!='H')
	  heavy_atoms2.push_back(tmp2[a]);
      }
      //   cout<<"DEB heavy atoms fag "<<j<<" : "<<heavy_atoms2.size()<<"\n";
      //loop on atoms
      double d2;
      bool done = false;
      for(int i1 = 0;i1<heavy_atoms.size();i1++){
	if(done) break;
	int a1 = heavy_atoms[i1];
	for(int i2 = 0;i2<heavy_atoms2.size();i2++){
	  if(done) break;
	  int a2 = heavy_atoms2[i2];
	  d2 = p.distanceP2(a1,a2);
	  if(d2<=cut){
	    done = true;
	  }
	}
      }
      if(done){
	m.set(i,j,1);
	m.set(j,i,1);
      }
    }
  }
  return m;
}


void randomize_cont_map(GSLMatrix &m, double tr){
  int size = m.rows();
  for(int i=0;i<size;i++){
    for(int j=i+3;j<size;j++){
      double r = drand48();
      if(r<tr){
	m.set(i,j,1-m.get(i,j));
	m.set(j,i,m.get(i,j));
      }
    }
  }
}



void assign_cont_map(GSLMatrix & dest,GSLMatrix &source){
  int size = source.rows();
  for(int i=0;i<size;i++){
    for(int j=i+3;j<size;j++){
      dest.set(i,j,source.get(i,j));
      dest.set(j,i,source.get(i,j));
    }
  }
}


double energy_cont_map(GSLVector & pe,  GSLMatrix & cm){
  int size = pe.size();
  GSLMatrix vv(size,size);
  GSLVector ev(size);
  eigenvv(cm, ev, vv);

  double n1 = 0;
  double n2 = 0;
  double ene = 0;
  for(int i=0;i<size;i++){
    n1 += pe.get(i)*pe.get(i);
    n2 += vv.get(i,0)*vv.get(i,0);
  }

  //  cout<<"DEBUG norm: "<<n1<<" "<<n2<<endl;

  for(int i=0;i<size;i++){
    double t = abs(pe.get(i))-abs(vv.get(i,0));
    ene +=t*t;
    //ene +=abs(t);
  }

  return ene/(n1*n2);
}

double energy_cont_map(GSLVector & pe, double orig_ev, GSLMatrix & cm){
  int size = pe.size();
  GSLMatrix vv(size,size);
  GSLVector ev(size);
  eigenvv(cm, ev, vv);

  double n1 = 0;
  double n2 = 0;
  double ene = 0;
  for(int i=0;i<size;i++){
    n1 += pe.get(i)*pe.get(i);
    n2 += vv.get(i,0)*vv.get(i,0);
  }

  //  cout<<"DEBUG norm: "<<n1<<" "<<n2<<endl;

  for(int i=0;i<size;i++){
    double t = abs(pe.get(i))-abs(vv.get(i,0));
    ene +=t*t;
    //ene +=abs(t);
  }

  return ene/(n1*n2) +(orig_ev-ev.get(0))*(orig_ev-ev.get(0));
}


double energy_cont_map(GSLVector & pe,GSLMatrix & cm,
		       Coor<double> & coor,
		       vector<vector<int> > &ha,
		       double c2){
 int size = pe.size();
 GSLMatrix vv(size,size);
 GSLVector ev(size);  

 //compute cm
 cm.set_zero();
 
 for(int i=0;i<size;i++){
   int s1 = ha[i].size();
   for(int j=i+3;j<size;j++){
     int s2 = ha[j].size();
     double d2;
     bool done = false;
     for(int ii=0;ii<s1;ii++){
       if(done) break;
       int a1 = Coor<double>::DIM*ha[i][ii];
       for(int jj = 0;jj<s2;jj++){
	 int a2 = Coor<double>::DIM*ha[j][jj];
	 d2 = 
	   (coor.coor[a1]-coor.coor[a2])*(coor.coor[a1]-coor.coor[a2])+
	   (coor.coor[a1+1]-coor.coor[a2+1])*(coor.coor[a1+1]-coor.coor[a2+1])+
	   (coor.coor[a1+2]-coor.coor[a2+2])*(coor.coor[a1+2]-coor.coor[a2+2]);
	 if(d2<c2){
	   done = true;
	   cm.set(i,j,1);
	   cm.set(j,i,1);
	 }
       }
     }
   }
 }

 eigenvv(cm, ev, vv);
 double n1 = 0;
 double n2 = 0;
 double ene = 0;
 for(int i=0;i<size;i++){
   n1 += pe.get(i)*pe.get(i);
   n2 += vv.get(i,0)*vv.get(i,0);
 }

 // cout<<"DEBUG norm: "<<n1<<" "<<n2<<endl;

 for(int i=0;i<size;i++){
   double t = abs(pe.get(i))-abs(vv.get(i,0));
   ene +=t*t;
   //ene +=abs(t);
 }

 return ene/(n1*n2);

}



double energy_cont_map_gauss(GSLMatrix & cm,
			     Coor<double> & coor,
			     vector<vector<int> > &ha,
			     double c2){
  int size = cm.rows();
  double ene = 0;
  for(int i=0;i<size;i++){
    int s1 = ha[i].size();
    for(int j=i+3;j<size;j++){
      //      if(cm.get(i,j)==0) continue;
      int s2 = ha[j].size();
      double d2;
      double d2_min=1000000;
      for(int ii=0;ii<s1;ii++){
	int a1 = Coor<double>::DIM*ha[i][ii];
	for(int jj = 0;jj<s2;jj++){
	  int a2 = Coor<double>::DIM*ha[j][jj];
	  d2 = 
	    (coor.coor[a1]-coor.coor[a2])*(coor.coor[a1]-coor.coor[a2])+
	    (coor.coor[a1+1]-coor.coor[a2+1])*(coor.coor[a1+1]-coor.coor[a2+1])+
	    (coor.coor[a1+2]-coor.coor[a2+2])*(coor.coor[a1+2]-coor.coor[a2+2]);
	  if(d2<d2_min){
	    d2_min = d2;
	  }
	}	
      }
      if(cm.get(i,j)==1){
	if(d2_min>c2){
	  ene +=abs(d2_min-c2);
	}
      }
      else{
	if(d2_min<c2){
	  ene +=abs(d2_min-c2);
	}
      }
	
    }
  } 
  return ene;
}


// GSLMatrix mc_find(GSLVector & pe, double tmp, int print,int steps, double tr){
//   int size = pe.size();
//   GSLMatrix cm(size,size);
//   cm.set_zero();

//   GSLMatrix cm_buffer(size,size);
//   cm_buffer.set_zero();

//   double ene;
//   double current_ene;

//   int accept = 0;
//   int reject = 0;
//   srand48(clock());
//   randomize_cont_map(cm_buffer,0.5);
//   assign_cont_map(cm,cm_buffer);
//   ene = energy_cont_map(pe,cm);

//   for(int i=0;i<steps;i++){
//     randomize_cont_map(cm_buffer,tr);
//     current_ene = energy_cont_map(pe,cm_buffer);
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     if(r<exp(-de)){
//       //accept
//       assign_cont_map(cm,cm_buffer);
//       ene = current_ene;
//       ++accept;
//     }
//     else{
//       assign_cont_map(cm_buffer,cm);
//       ++reject;
//     }
//     if((i+1)%print == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//     }
//   }
//   return cm;
// }


// GSLMatrix mc_find_cont(GSLMatrix & start,GSLVector & pe, 
// 		       double tmp, int print,int steps, double tr){
//   //  cout<<"DEB cont "<<endl;
//   int size = pe.size();
//   GSLMatrix cm(size,size);
//   cm.set_zero();

//   GSLMatrix cm_buffer(size,size);
//   cm_buffer.set_zero();

//   double ene;
//   double current_ene;

//   int accept = 0;
//   int reject = 0;
//   srand48(clock());
//   //  randomize_cont_map(cm_buffer,0.5);
//   assign_cont_map(cm,start);
//   assign_cont_map(cm_buffer,start);
//   ene = energy_cont_map(pe,cm);

//   for(int i=0;i<steps;i++){
//     randomize_cont_map(cm_buffer,tr);
//     current_ene = energy_cont_map(pe,cm_buffer);
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     if(r<exp(-de)){
//       //accept
//       assign_cont_map(cm,cm_buffer);
//       ene = current_ene;
//       ++accept;
//     }
//     else{
//       assign_cont_map(cm_buffer,cm);
//       ++reject;
//     }
//     if((i+1)%print == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//     }
//   }
//   return cm;  
// }


// GSLMatrix mc_find_ev(GSLVector & pe, double ev,
// 		     double tmp, int print,int steps, double tr){
//   int size = pe.size();
//   GSLMatrix cm(size,size);
//   cm.set_zero();

//   GSLMatrix cm_buffer(size,size);
//   cm_buffer.set_zero();

//   double ene;
//   double current_ene;

//   int accept = 0;
//   int reject = 0;
//   srand48(clock());
//   randomize_cont_map(cm_buffer,tr);
//   assign_cont_map(cm,cm_buffer);
//   ene = energy_cont_map(pe,ev,cm);

//   for(int i=0;i<steps;i++){
//     randomize_cont_map(cm_buffer,tr);
//     current_ene = energy_cont_map(pe,ev,cm_buffer);
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     if(r<exp(-de)){
//       //accept
//       assign_cont_map(cm,cm_buffer);
//       ene = current_ene;
//       ++accept;
//     }
//     else{
//       assign_cont_map(cm_buffer,cm);
//       ++reject;
//     }
//     if((i+1)%print == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//     }
//   }
//   return cm;
// }




// struct DummyData {
//   int nsteps;
//   int nsave;
//   int save_fitted;
//   void write(ofstream & rstOut){
//     rstOut.write((char *)this, sizeof(DummyData));
//   }

//   void read(ifstream & rstIn){
//     rstIn.read(( char *)this, sizeof(DummyData));
//   }

// };


// struct PeOptions {
//   double cut_on;
//   double cut_off;
//   int steps;
//   int nprint;
//   int nsave;
//   double temp;
//   int seed;
//   string move_set;
//   string trj;
//   string rst;
//   double cm_cut_off;

//   double rand_temp;
//   int rand_steps;

//   double lambda;

//   PeOptions(){
//   cut_on = 6.5;
//   cut_off = 7.0;
//   steps= 1000;
//   nprint = 1000;
//   nsave = 1000;
//   temp = 0.01;
//   seed = 30869;
//   move_set = "moveset.def";
//   trj = "guess.dcd";
//   rst = "guess.rst";
//   cm_cut_off = 4.5;

//   rand_temp = 1;
//   rand_steps = 1000;

//   lambda = 0.0001;
//   }
// };


// ostream & operator<<(ostream & out,const PeOptions & opts){
//   out<<"\t>> Cut on:                      "
//      <<opts.cut_on<<"\n";
//   out<<"\t>> Cut off:                     "
//      <<opts.cut_off<<"\n";
//   out<<"\t>> Move set definitions:        "
//     +opts.move_set+"\n";
//   out<<"\t>> Trajectory:                  "
//     +opts.trj+"\n";
//   out<<"\t>> Restart:                     "
//     +opts.rst+"\n";
//   out<<"\t>> Number of steps:             "
//      <<opts.steps<<"\n";
//   out<<"\t>> Print frequency:             "
//      <<opts.nprint<<"\n";
//   out<<"\t>> Save frequency:              "
//      <<opts.nsave<<"\n";
//   out<<"\t>> Temperature:                 "
//      <<opts.temp<<"\n";
//   out<<"\t>> CM Cut off                   "
//      <<opts.cm_cut_off<<"\n";
//   out<<"\t>> RAND Number of steps:        "
//      <<opts.rand_steps<<"\n";
//   out<<"\t>> RAND Temperature:            "
//      <<opts.rand_temp<<"\n";
//   out<<"\t>> Lambda:                      "
//      <<opts.lambda<<"\n";

//   return out;
// }

// inline string to_string(const PeOptions & options){
//   stringstream s;
//   s<<options;
//   return s.str();
// }

// GSLMatrix mc_find_ff(GSLVector & pe,
// 		     const Molecules &molecules,
// 		     const MDB  mdb,
// 		     const PeOptions & options,
// 		     bool  linearize){
//   SC19HSFFOptions<double> opts;
//   opts.nbond_options.nblist_options.CUTON = options.cut_on;
//   opts.nbond_options.nblist_options.CUTOFF = options.cut_off;

//   SC19HSForceField<double> ff(molecules,mdb,opts);
//   MoveSet<double> move_set(molecules,options.move_set);
  
//   DummyData dd;
//   dd.nsteps = options.steps;
//   dd.nsave = options.nsave;
//   dd.save_fitted = 0;

//   int size = molecules.fragment_size();
//   //step 1
//   GSLMatrix cm(size,size);

//   //step 2
//   CoorManagerOptions<double> cm_opts;
//   cm_opts.trj = options.trj;
//   cm_opts.rst = options.rst;
//   Molecules start_m = molecules;
//   if(linearize){
//     for(int i=0;i<start_m.protein_size();i++){
//       start_m.protein(i).linearize();
//     }
//   }
//   CoorManager<double> coor_manager(start_m,cm_opts);
//   coor_manager.write_header(dd);
//   //step 3 local data
//   vector<vector<int> > Ha;
//   Ha.resize(size);
//   {
//     for(int i=0;i<size;i++){
//       vector<int> frag_atoms = molecules.fragment_atoms(i);
//       for(int a = 0; a< frag_atoms.size();a++){
// 	if(molecules[frag_atoms[a]].name()[0]!='H'){
// 	  Ha[i].push_back(frag_atoms[a]);
// 	}
//       }
//     }
//   }

//   double c2 = options.cm_cut_off*options.cm_cut_off;
//   double ff_ene = ff.energy(coor_manager.atomCoor);
//   double const_ene = energy_cont_map(pe,
// 				     cm,
// 				     coor_manager.atomCoor,
// 				     Ha,
// 				     c2);

//   //do monte carlo
//   srand48(clock());
//   double lambda = options.lambda;
//   double ene = lambda*ff_ene + const_ene;
//   cout<<"\t>> Initial energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<endl;

//   {
//     cout<<"\t>> Doing randomization\n";
//     int accept = 0;
//     int reject = 0;

//     double current_ff_ene;
//     double current_const_ene;
//     double current_ene;

  
//     for(int i=0;i<options.rand_steps;i++){
//       move_set.apply(coor_manager.atomCoor);
//       current_ff_ene = ff.energy(coor_manager.atomCoor);
//       current_const_ene =  energy_cont_map(pe,
// 					   cm,
// 					   coor_manager.atomCoor,
// 					   Ha,
// 					   c2);
//       current_ene = lambda*current_ff_ene + current_const_ene;
//       double r = drand48();
//       double de = (current_ene-ene)/options.rand_temp;
//       //accept
//       if(r<exp(-de)){
// 	coor_manager.accept();
// 	ene = current_ene;
// 	ff_ene = current_ff_ene;
// 	const_ene = current_const_ene;
// 	++accept;
//       }
//       else{
// 	coor_manager.reject();
// 	++reject;
//       }
//       if((i+1)%1000 == 0){
// 	cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<"\n";
// 	cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//       }

//     }
//   }
//   int accept = 0;
//   int reject = 0;

//   double current_ff_ene;
//   double current_const_ene;
//   double current_ene;
//   double tmp = options.temp;
//   int steps = options.steps;
//   cout<<"\t>> Starting ......."<<endl;
//   for(int i=0;i<steps;i++){
//     move_set.apply(coor_manager.atomCoor);
//     current_ff_ene = ff.energy(coor_manager.atomCoor);
//     current_const_ene =  energy_cont_map(pe,
// 					 cm,
// 					 coor_manager.atomCoor,
// 					 Ha,
// 					 c2);
//     current_ene = lambda*current_ff_ene + current_const_ene;
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     //accept
//     if(r<exp(-de)){
//       coor_manager.accept();
//       ene = current_ene;
//       ff_ene = current_ff_ene;
//       const_ene = current_const_ene;
//       ++accept;
//     }
//     else{
//       coor_manager.reject();
//       ++reject;
//     }
//     if((i+1)%options.nprint == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<"\n";
// //       cout<<"\t>> step "<<i+1<<" energy: "<<ff.data().internal_energy<<" "<<
// // 	ff.data().vdw_energy<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/(accept+reject)<<"\n";
//     }

//     if((i+1)%options.nsave == 0){
//       coor_manager.save_trj(dd);
//       coor_manager.write_rst(dd);
// //       if(((double)accept)/(accept+reject)<0.5) tmp =1.1*tmp;
// //       else
// // 	tmp =0.9*tmp;
// //       cout<<"\t>> Temperature set to: "<<tmp<<"\n"; 
//       accept = reject = 0;
//     }

//   }
//   return cm;

    
// }


// GSLMatrix mc_find_ff_cont(GSLVector & pe,
// 			  const Molecules &molecules,
// 			  const MDB  mdb,
// 			  const PeOptions & options,
// 			  string rst_file){
//   SC19HSFFOptions<double> opts;
//   opts.nbond_options.nblist_options.CUTON = options.cut_on;
//   opts.nbond_options.nblist_options.CUTOFF = options.cut_off;

//   SC19HSForceField<double> ff(molecules,mdb,opts);
//   MoveSet<double> move_set(molecules,options.move_set);
  
//   DummyData dd;
//   dd.nsteps = options.steps;
//   dd.nsave = options.nsave;
//   dd.save_fitted = 0;

//   int size = molecules.fragment_size();
//   //step 1
//   GSLMatrix cm(size,size);

//   //step 2
//   CoorManagerOptions<double> cm_opts;
//   cm_opts.trj = options.trj;
//   cm_opts.rst = options.rst;
//   CoorManager<double> coor_manager(molecules,cm_opts,true);
//   coor_manager.write_header(dd);
//   coor_manager.read_rst(rst_file,dd);
//   //step 3 local data
//   vector<vector<int> > Ha;
//   Ha.resize(size);
//   {
//     for(int i=0;i<size;i++){
//       vector<int> frag_atoms = molecules.fragment_atoms(i);
//       for(int a = 0; a< frag_atoms.size();a++){
// 	if(molecules[frag_atoms[a]].name()[0]!='H'){
// 	  Ha[i].push_back(frag_atoms[a]);
// 	}
//       }
//     }
//   }

//   double c2 = options.cm_cut_off*options.cm_cut_off;
//   double ff_ene = ff.energy(coor_manager.atomCoor);
//   double const_ene = energy_cont_map(pe,
// 				     cm,
// 				     coor_manager.atomCoor,
// 				     Ha,
// 				     c2);

//   //do monte carlo
//   srand48(clock());
//   double lambda = options.lambda;
//   double ene = lambda*ff_ene + const_ene;
//   cout<<"\t>> Initial energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<endl;

//   {
//     cout<<"\t>> Doing randomization\n";
//     int accept = 0;
//     int reject = 0;

//     double current_ff_ene;
//     double current_const_ene;
//     double current_ene;

  
//     for(int i=0;i<options.rand_steps;i++){
//       move_set.apply(coor_manager.atomCoor);
//       current_ff_ene = ff.energy(coor_manager.atomCoor);
//       current_const_ene =  energy_cont_map(pe,
// 					   cm,
// 					   coor_manager.atomCoor,
// 					   Ha,
// 					   c2);
//       current_ene = lambda*current_ff_ene + current_const_ene;
//       double r = drand48();
//       double de = (current_ene-ene)/options.rand_temp;
//       //accept
//       if(r<exp(-de)){
// 	coor_manager.accept();
// 	ene = current_ene;
// 	ff_ene = current_ff_ene;
// 	const_ene = current_const_ene;
// 	++accept;
//       }
//       else{
// 	coor_manager.reject();
// 	++reject;
//       }
//       if((i+1)%1000 == 0){
// 	cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<"\n";
// 	cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//       }

//     }
//   }
//   int accept = 0;
//   int reject = 0;

//   double current_ff_ene;
//   double current_const_ene;
//   double current_ene;
//   double tmp = options.temp;
//   int steps = options.steps;
//   cout<<"\t>> Starting ......."<<endl;
//   for(int i=0;i<steps;i++){
//     move_set.apply(coor_manager.atomCoor);
//     current_ff_ene = ff.energy(coor_manager.atomCoor);
//     current_const_ene =  energy_cont_map(pe,
// 					 cm,
// 					 coor_manager.atomCoor,
// 					 Ha,
// 					 c2);
//     current_ene = lambda*current_ff_ene + current_const_ene;
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     //accept
//     if(r<exp(-de)){
//       coor_manager.accept();
//       ene = current_ene;
//       ff_ene = current_ff_ene;
//       const_ene = current_const_ene;
//       ++accept;
//     }
//     else{
//       coor_manager.reject();
//       ++reject;
//     }
//     if((i+1)%options.nprint == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<"\n";
// //       cout<<"\t>> step "<<i+1<<" energy: "<<ff.data().internal_energy<<" "<<
// // 	ff.data().vdw_energy<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/(accept+reject)<<"\n";
//     }

//     if((i+1)%options.nsave == 0){
//       coor_manager.save_trj(dd);
//       coor_manager.write_rst(dd);
//       accept = reject = 0;
//     }

//   }
//   return cm;

    
// }



// GSLMatrix rex_find_ff(GSLVector & pe,
// 		      const Molecules &molecules,
// 		      const MDB  mdb,
// 		      const vector<PeOptions> & options,
// 		      bool  linearize,
// 		      int rexfq){
//   SC19HSFFOptions<double> opts;
//   opts.nbond_options.nblist_options.CUTON = options[0].cut_on;
//   opts.nbond_options.nblist_options.CUTOFF = options[0].cut_off;

//   SC19HSForceField<double> ff(molecules,mdb,opts);
//   MoveSet<double> move_set(molecules,options[0].move_set);
  
//   DummyData dd;
//   dd.nsteps = options[0].steps;
//   dd.nsave = options[0].nsave;
//   dd.save_fitted = 0;

//   int size = molecules.fragment_size();
//   //step 1
//   GSLMatrix cm(size,size);

//   //step 2
//   Molecules start_m = molecules;
//   if(linearize){
//     for(int i=0;i<start_m.protein_size();i++){
//       start_m.protein(i).linearize();
//     }
//   }
//   vector<CoorManager<double>* > coor_manager;
//   coor_manager.resize(options.size());
//   for(int i=0;i<coor_manager.size();i++){
//     CoorManagerOptions<double> cm_opts;
//     cm_opts.trj = options[i].trj;
//     cm_opts.rst = options[i].rst;
//     coor_manager[i] = new CoorManager<double>(start_m,cm_opts);
//     coor_manager[i]->write_header(dd);
//   }
  
//   //step 3 local data
//   vector<vector<int> > Ha;
//   Ha.resize(size);
//   {
//     for(int i=0;i<size;i++){
//       vector<int> frag_atoms = molecules.fragment_atoms(i);
//       for(int a = 0; a< frag_atoms.size();a++){
// 	if(molecules[frag_atoms[a]].name()[0]!='H'){
// 	  Ha[i].push_back(frag_atoms[a]);
// 	}
//       }
//     }
//   }

//   double c2 = options[0].cm_cut_off*options[0].cm_cut_off;
//   vector<double> ff_ene; ff_ene.resize(coor_manager.size());
//   vector<double> const_ene; const_ene.resize(coor_manager.size());
//   for(int i=0;i<coor_manager.size();i++){
//     ff_ene[i] = ff.energy(coor_manager[i]->atomCoor);
//     const_ene [i] = energy_cont_map(pe,
// 				    cm,
// 				    coor_manager[i]->atomCoor,
// 				    Ha,
// 				    c2);
//   }
//   //do monte carlo
//   srand48(clock());
//   double lambda = options[0].lambda;
//   vector<double> ene; ene.resize(coor_manager.size());

//   //rex_size
//   int rex_size = coor_manager.size();
//   for(int i = 0;i<rex_size;i++){
//     ene[i] = lambda*ff_ene[i] + const_ene[i];
//     cout<<"\t>> Initial energy: "<<i<<" "
// 	<<ene[i]<<" "<<ff_ene[i]<<" "<<const_ene[i]<<endl;
//   }
//   {
//     cout<<"\t>> Doing randomization\n";
//     vector<int> accept; accept.resize(rex_size);
//     vector<int> reject; reject.resize(rex_size);

//     vector<double> current_ff_ene; current_ff_ene.resize(rex_size);
//     vector<double> current_const_ene; current_const_ene.resize(rex_size);
//     vector<double> current_ene; current_ene.resize(rex_size);

  
//     for(int i=0;i<options[0].rand_steps;i++){
//       for(int rex = 0;rex<rex_size;rex++){
// 	move_set.apply(coor_manager[rex]->atomCoor);
// 	current_ff_ene[rex] = ff.energy(coor_manager[rex]->atomCoor);
// 	current_const_ene[rex] =  energy_cont_map(pe,
// 						  cm,
// 						  coor_manager[rex]->atomCoor,
// 						  Ha,
// 						  c2);
// 	current_ene[rex] = lambda*current_ff_ene[rex] + current_const_ene[rex];
// 	double r = drand48();
// 	double de = (current_ene[rex]-ene[rex])/options[rex].rand_temp;
// 	//accept
// 	if(r<exp(-de)){
// 	  coor_manager[rex]->accept();
// 	  ene[rex] = current_ene[rex];
// 	  ff_ene[rex] = current_ff_ene[rex];
// 	  const_ene[rex] = current_const_ene[rex];
// 	  ++accept[rex];
// 	}
// 	else{
// 	  coor_manager[rex]->reject();
// 	  ++reject[rex];
// 	}
// 	if((i+1)%1000 == 0){
// 	  cout<<"\t>> step "<<i+1<<" energy: "<<rex<<" "<<ene[rex]
// 	      <<" "<<ff_ene[rex]<<" "<<const_ene[rex]<<"\n";
// 	  cout<<"\t>> step "<<i+1<<" accept: "<<rex<<" "
// 	      <<((double)accept[rex])/((double)i+1)<<"\n";
// 	}
//       }
//     }
//   }
//   vector<int> accept; accept.resize(rex_size);
//   vector<int> reject; reject.resize(rex_size);

//   vector<double> current_ff_ene; current_ff_ene.resize(rex_size);
//   vector<double> current_const_ene; current_const_ene.resize(rex_size);
//   vector<double> current_ene; current_ene.resize(rex_size);
//   vector<double> tmp; tmp.resize(rex_size);
//   for(int i=0;i<rex_size;i++) tmp[i] = options[i].temp;
//   int steps = options[0].steps;
//   cout<<"\t>> Starting ......."<<endl;
//   for(int i=0;i<steps;i++){
//     for(int rex = 0;rex<rex_size;rex++){
//       move_set.apply(coor_manager[rex]->atomCoor);
//       current_ff_ene[rex] = ff.energy(coor_manager[rex]->atomCoor);
//       current_const_ene[rex] =  energy_cont_map(pe,
// 						cm,
// 						coor_manager[rex]->atomCoor,
// 						Ha,
// 						c2);
//       current_ene[rex] = lambda*current_ff_ene[rex] + current_const_ene[rex];
//       double r = drand48();
//       double de = (current_ene[rex]-ene[rex])/tmp[rex];
//       //accept
//       if(r<exp(-de)){
// 	coor_manager[rex]->accept();
// 	ene[rex] = current_ene[rex];
// 	ff_ene[rex] = current_ff_ene[rex];
// 	const_ene[rex] = current_const_ene[rex];
// 	++accept[rex];
//       }
//       else{
// 	coor_manager[rex]->reject();
// 	++reject[rex];
//       }

//       if((i+1)%options[rex].nprint == 0){
// 	cout<<"\t>> step "<<i+1<<" energy: "<<rex<<" "<<ene[rex]
// 	    <<" "<<ff_ene[rex]<<" "<<const_ene[rex]<<"\n";
// 	cout<<"\t>> step "<<i+1<<" accept: "<<rex<<" "
// 	    <<accept[rex]<<"\n";
// 	cout<<"\t>> step "<<i+1<<" reject: "<<rex<<" "
// 	    <<reject[rex]<<"\n";
// 	cout<<"\t>> step "<<i+1<<" tem:    "<<rex<<" "<<tmp[rex]<<"\n";
//       }

//       if((i+1)%options[rex].nsave == 0){
// 	coor_manager[rex]->save_trj(dd);
// 	coor_manager[rex]->write_rst(dd);
// 	accept[rex] = reject[rex] = 0;
//       }
//     }
//     //exchange 
//     if((i+1)%rexfq == 0){
//       int q = (i+1)/rexfq;
//       if((q%2)==0){
// 	cout<<"\t>> Trying even ...\n";
// 	for(int k = 0;k+2<=rex_size;k=k+2){
// 	  double p = (ene[k+1]-ene[k])*(1/tmp[k]-1/tmp[k+1]);
// 	  double r = drand48();
// 	  if(r<exp(-p)){
// 	    //exchange
// 	    cout<<"\t>> Exchange: "<<k<<"-"<<k+1<<"\n";
// 	    double t= tmp[k];
// 	    tmp[k] = tmp[k+1];
// 	    tmp[k+1] = t;
// 	  }
// 	}
//       }
//       else {
// 	cout<<"\t>> Trying odd ...\n";
// 	for(int k = 1;k+2<=rex_size;k=k+2){
// 	  double p = (ene[k+1]-ene[k])*(1/tmp[k]-1/tmp[k+1]);
// 	  double r = drand48();
// 	  if(r<exp(-p)){
// 	    //exchange
// 	    cout<<"\t>> Exchange: "<<k<<"-"<<k+1<<"\n";
// 	    double t= tmp[k];
// 	    tmp[k] = tmp[k+1];
// 	    tmp[k+1] = t;
// 	  }
// 	}
//       }
//     }
//   }
//   for(int i=0;i<coor_manager.size();i++){
//     delete coor_manager[i];
//   }
//   return cm;
// }



// Molecules mc_find_struct(GSLMatrix & cm,
// 			 const Molecules &molecules,
// 			 const MDB  mdb,
// 			 const PeOptions & options){
//   SC19HSFFOptions<double> opts;
//   opts.nbond_options.nblist_options.CUTON = options.cut_on;
//   opts.nbond_options.nblist_options.CUTOFF = options.cut_off;

//   SC19HSForceField<double> ff(molecules,mdb,opts);
//   MoveSet<double> move_set(molecules,options.move_set);
  
//   DummyData dd;
//   dd.nsteps = options.steps;
//   dd.nsave = options.nsave;
//   dd.save_fitted = 0;

//   int size = molecules.fragment_size();

//   //step 2
//   CoorManagerOptions<double> cm_opts;
//   cm_opts.trj = options.trj;
//   cm_opts.rst = options.rst;
//   CoorManager<double> coor_manager(molecules,cm_opts);
//   coor_manager.write_header(dd);
//   //step 3 local data
//   vector<vector<int> > Ha;
//   Ha.resize(size);
//   {
//     for(int i=0;i<size;i++){
//       vector<int> frag_atoms = molecules.fragment_atoms(i);
//       for(int a = 0; a< frag_atoms.size();a++){
// 	if(molecules[frag_atoms[a]].name()[0]!='H'){
// 	  Ha[i].push_back(frag_atoms[a]);
// 	}
//       }
//     }
//   }

//   double c2 = options.cm_cut_off*options.cm_cut_off;
//   double ff_ene = ff.energy(coor_manager.atomCoor);
//   double const_ene = energy_cont_map_gauss(cm,
// 					   coor_manager.atomCoor,
// 					   Ha,
// 					   c2);

//   //do monte carlo
//   srand48(clock());
//   double lambda = options.lambda;
//   double ene = lambda*ff_ene + const_ene;
//   cout<<"\t>> Initial energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<endl;

//   {
//     cout<<"\t>> Doing randomization\n";
//     int accept = 0;
//     int reject = 0;

//     double current_ff_ene;
//     double current_const_ene;
//     double current_ene;

//     ff_ene = ff.energy(coor_manager.atomCoor);
//     ene = lambda*ff_ene;
//     for(int i=0;i<options.rand_steps;i++){
//       move_set.apply(coor_manager.atomCoor);
//       current_ff_ene = ff.energy(coor_manager.atomCoor);
// //       current_const_ene =  energy_cont_map_gauss(cm,
// // 						 coor_manager.atomCoor,
// // 						 Ha,
// // 						 c2);
//       current_ene = lambda*current_ff_ene;// + current_const_ene;
//       double r = drand48();
//       double de = (current_ene-ene)/options.rand_temp;
//       //accept
//       if(r<exp(-de)){
// 	coor_manager.accept();
// 	ene = current_ene;
// 	ff_ene = current_ff_ene;
// 	const_ene = current_const_ene;
// 	++accept;
//       }
//       else{
// 	coor_manager.reject();
// 	++reject;
//       }
//       if((i+1)%1000 == 0){
// 	cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<"\n";
// 	cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//       }

//     }
//   }
//   int accept = 0;
//   int reject = 0;

//   double current_ff_ene;
//   double current_const_ene;
//   double current_ene;
//   double tmp = options.temp;
//   int steps = options.steps;
//   cout<<"\t>> Starting ......."<<endl;
//   ff_ene = ff.energy(coor_manager.atomCoor);
//   const_ene = energy_cont_map_gauss(cm,
// 				    coor_manager.atomCoor,
// 				    Ha,
// 				    c2);
//   ene = lambda*ff_ene + const_ene;
//   for(int i=0;i<steps;i++){
//     move_set.apply(coor_manager.atomCoor);
//     current_ff_ene = ff.energy(coor_manager.atomCoor);
//     current_const_ene =  energy_cont_map_gauss(cm,
// 					       coor_manager.atomCoor,
// 					       Ha,
// 					       c2);
//     current_ene = lambda*current_ff_ene + current_const_ene;
//     double r = drand48();
//     double de = (current_ene-ene)/tmp;
//     //accept
//     if(r<exp(-de)){
//       coor_manager.accept();
//       ene = current_ene;
//       ff_ene = current_ff_ene;
//       const_ene = current_const_ene;
//       ++accept;
//     }
//     else{
//       coor_manager.reject();
//       ++reject;
//     }
//     if((i+1)%options.nprint == 0){
//       cout<<"\t>> step "<<i+1<<" energy: "<<ene<<" "<<ff_ene<<" "<<const_ene<<"\n";
// //       cout<<"\t>> step "<<i+1<<" energy: "<<ff.data().internal_energy<<" "<<
// // 	ff.data().vdw_energy<<"\n";
//       cout<<"\t>> step "<<i+1<<" accept: "<<((double)accept)/((double)i+1)<<"\n";
//     }

//     if((i+1)%options.nsave == 0){
//       coor_manager.save_trj(dd);
//       coor_manager.write_rst(dd);
//     }

//   }

//   Molecules res = molecules;
//   for(int i=0;i<res.atom_size();i++){
//     int pos = Coor<double>::DIM*i;
//     res[i].set_x(coor_manager.atomCoor.coor[pos]);
//     res[i].set_y(coor_manager.atomCoor.coor[pos+1]);
//     res[i].set_z(coor_manager.atomCoor.coor[pos+2]);
//   }

//   return res;
// }


//analysis
// class gsl_energy {
//   GSLMatrix cm;
//   GSLVector pe;
//   GSLMatrix l_cm;
//   double    c2;
//   double    cutoff;
//   SC19HSForceField<double> *ff;
//   vector<vector<int> > Ha;
// public:
//   gsl_energy(const Molecules & molecules,const MDB & mdb,
// 	     const GSLMatrix & c_map,
// 	     double cmap_cutoff, double ff_cutoff):
//     cm(c_map.rows(),c_map.cols()),pe(c_map.rows()),
//     l_cm(c_map.rows(),c_map.cols())
//   {
//     c2 = cmap_cutoff*cmap_cutoff;
//     cutoff = ff_cutoff;
//     {
//       int count = 0;
//       for(int i=0;i<cm.rows();i++){
// 	for(int j=0;j<cm.cols();j++){
// 	  cm.set(i,j,c_map.get(i,j));
// 	  if(cm.get(i,j)>0) ++count;
// 	}
//       }
//       cout<<"\t>> Contact map has "<<count/2<<" contacts\n";
//       int size = c_map.rows();
//       GSLMatrix vv(size,size);
//       GSLVector ev(size);
//       eigenvv(cm, ev, vv);
//       for(int i=0;i<pe.size();i++){
// 	pe.set(i,vv.get(i,0));
//       }
//     }

//     //step 3 local data
//     {
//       int size = c_map.rows();
//       Ha.resize(size);
//       {
// 	for(int i=0;i<size;i++){
// 	  vector<int> frag_atoms = molecules.fragment_atoms(i);
// 	  for(int a = 0; a< frag_atoms.size();a++){
// 	    if(molecules[frag_atoms[a]].name()[0]!='H'){
// 	      Ha[i].push_back(frag_atoms[a]);
// 	    }
// 	  }
// 	}
//       }
//     }
//     //ff
//     SC19HSFFOptions<double> opts;
//     opts.nbond_options.nblist_options.CUTON = cutoff;
//     opts.nbond_options.nblist_options.CUTOFF = cutoff;
//     ff = new SC19HSForceField<double>(molecules,mdb,opts);
//   }

//   ~gsl_energy(){ delete ff;}
  
//   double ff_energy(const Coor<float> & coor){
//     Coor<double> l_coor;
//     l_coor.resize(coor.size());
//     for(int i=0;i<coor.size();i++){
//       int pos = Coor<double>::DIM*i;
//       l_coor.coor[pos] = coor.coor[pos];
//       l_coor.coor[pos+1] = coor.coor[pos+1];
//       l_coor.coor[pos+2] = coor.coor[pos+2];
//     }
//     return ff->energy(l_coor);
//   }

//   double const_energy(const Coor<float> & coor){
//     Coor<double> l_coor;
//     l_coor.resize(coor.size());
//     for(int i=0;i<coor.size();i++){
//       int pos = Coor<double>::DIM*i;
//       l_coor.coor[pos] = coor.coor[pos];
//       l_coor.coor[pos+1] = coor.coor[pos+1];
//       l_coor.coor[pos+2] = coor.coor[pos+2];
//     }
//     return energy_cont_map(pe,
// 			   l_cm,
// 			   l_coor,
// 			   Ha,
// 			   c2);
//   }

//   double matches(const Coor<float> & coor){
//     Coor<double> l_coor;
//     int size = cm.rows();
//     l_coor.resize(coor.size());
//     for(int i=0;i<coor.size();i++){
//       int pos = Coor<double>::DIM*i;
//       l_coor.coor[pos] = coor.coor[pos];
//       l_coor.coor[pos+1] = coor.coor[pos+1];
//       l_coor.coor[pos+2] = coor.coor[pos+2];
//     }
//     l_cm.set_zero();
//     for(int i=0;i<size;i++){
//       int s1 = Ha[i].size();
//       for(int j=i+3;j<size;j++){
// 	int s2 = Ha[j].size();
// 	double d2;
// 	bool done = false;
// 	for(int ii=0;ii<s1;ii++){
// 	  if(done) break;
// 	  int a1 = Coor<double>::DIM*Ha[i][ii];
// 	  for(int jj = 0;jj<s2;jj++){
// 	    int a2 = Coor<double>::DIM*Ha[j][jj];
// 	    d2 = 
// 	      (l_coor.coor[a1]-l_coor.coor[a2])*(l_coor.coor[a1]-l_coor.coor[a2])+
// 	      (l_coor.coor[a1+1]-l_coor.coor[a2+1])*(l_coor.coor[a1+1]-l_coor.coor[a2+1])+
// 	      (l_coor.coor[a1+2]-l_coor.coor[a2+2])*(l_coor.coor[a1+2]-l_coor.coor[a2+2]);
// 	    if(d2<c2){
// 	      done = true;
// 	      l_cm.set(i,j,1);
// 	      l_cm.set(j,i,1);
// 	    }
// 	  }
// 	}
//       }
//     }

//     int m = 0;
//     for(int i=0;i<l_cm.rows();i++){
//       for(int j=i+3;j<l_cm.cols();j++){
// 	if(l_cm.get(i,j)==cm.get(i,j)) ++m;
//       }
//     }
//     return m;
//   }

//   int false_pos(const Coor<float> & coor){
//     Coor<double> l_coor;
//     l_coor.resize(coor.size());
//     for(int i=0;i<coor.size();i++){
//       int pos = Coor<double>::DIM*i;
//       l_coor.coor[pos] = coor.coor[pos];
//       l_coor.coor[pos+1] = coor.coor[pos+1];
//       l_coor.coor[pos+2] = coor.coor[pos+2];
//     }
//     int size = cm.rows();
//     l_cm.set_zero();
//     //    int count=0;
//     for(int i=0;i<size;i++){
//       int s1 = Ha[i].size();
//       for(int j=i+3;j<size;j++){
// 	int s2 = Ha[j].size();
// 	double d2;
// 	bool done = false;
// 	for(int ii=0;ii<s1;ii++){
// 	  if(done) break;
// 	  int a1 = Coor<double>::DIM*Ha[i][ii];
// 	  for(int jj = 0;jj<s2;jj++){
// 	    int a2 = Coor<double>::DIM*Ha[j][jj];
// 	    d2 = 
// 	      (l_coor.coor[a1]-l_coor.coor[a2])*(l_coor.coor[a1]-l_coor.coor[a2])+
// 	      (l_coor.coor[a1+1]-l_coor.coor[a2+1])*(l_coor.coor[a1+1]-l_coor.coor[a2+1])+
// 	      (l_coor.coor[a1+2]-l_coor.coor[a2+2])*(l_coor.coor[a1+2]-l_coor.coor[a2+2]);
// 	    if(d2<c2){
// 	      done = true;
// 	      l_cm.set(i,j,1);
// 	      l_cm.set(j,i,1);
// 	      //      ++count;
// 	      //	      cout<<"DEBUG "<<l_cm.get(i,j)<<" "<<cm.get(i,j)<<endl;
// 	    }
// 	  }
// 	}
//       }
//     }
//     cout<<"\t>> Current number of contacts "<<count<<"\n";
//     count = 0;
//     for(int i=0;i<l_cm.rows();i++){
//       for(int j=i+3;j<l_cm.cols();j++){
// 	if(cm.get(i,j)==l_cm.get(i,j)){
// 	  ++count;
// 	}
//       }
//     }
//     cout<<"\t>> Current number of match "<<count<<"\n";
//     int m = 0;
//     for(int i=0;i<l_cm.rows();i++){
//       for(int j=i+3;j<l_cm.cols();j++){
// 	//	cout<<"DEBUG "<<l_cm.get(i,j)<<" "<<cm.get(i,j)<<" "<<((l_cm.get(i,j)>0)&&(cm.get(i,j)==0))<<endl;
// 	if((l_cm.get(i,j)>0)&&(cm.get(i,j)==0)) ++m;
//       }
//     }
//     return m;
//   }

//   int false_neg(const Coor<float> & coor){
//     Coor<double> l_coor;
//     l_coor.resize(coor.size());
//     for(int i=0;i<coor.size();i++){
//       int pos = Coor<double>::DIM*i;
//       l_coor.coor[pos] = coor.coor[pos];
//       l_coor.coor[pos+1] = coor.coor[pos+1];
//       l_coor.coor[pos+2] = coor.coor[pos+2];
//     }
//     l_cm.set_zero();
//     int size = cm.rows();
//     for(int i=0;i<size;i++){
//       int s1 = Ha[i].size();
//       for(int j=i+3;j<size;j++){
// 	int s2 = Ha[j].size();
// 	double d2;
// 	bool done = false;
// 	for(int ii=0;ii<s1;ii++){
// 	  if(done) break;
// 	  int a1 = Coor<double>::DIM*Ha[i][ii];
// 	  for(int jj = 0;jj<s2;jj++){
// 	    int a2 = Coor<double>::DIM*Ha[j][jj];
// 	    d2 = 
// 	      (l_coor.coor[a1]-l_coor.coor[a2])*(l_coor.coor[a1]-l_coor.coor[a2])+
// 	      (l_coor.coor[a1+1]-l_coor.coor[a2+1])*(l_coor.coor[a1+1]-l_coor.coor[a2+1])+
// 	      (l_coor.coor[a1+2]-l_coor.coor[a2+2])*(l_coor.coor[a1+2]-l_coor.coor[a2+2]);
// 	    if(d2<c2){
// 	      done = true;
// 	      l_cm.set(i,j,1);
// 	      l_cm.set(j,i,1);
// 	    }
// 	  }
// 	}
//       }
//     }

//     int m = 0;
//     for(int i=0;i<l_cm.rows();i++){
//       for(int j=i+3;j<l_cm.cols();j++){
// 	if((l_cm.get(i,j)==0)&&(cm.get(i,j)>0)) ++m;
//       }
//     }
//     return m;
//   }

// };
    
template<>
inline string to_string(const Molecules & molecules){
  return "<Molecules>";
}


// inline string to_string(const gsl_energy &v){
//   return "<gsl::energy>";
// }

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

    Module mod = Module("gsl");
    Class<GSLVector,ArgList1<int> >(mod.self(),
				    "vector")
      .def_method("size",&GSLVector::size)
      .def_method("get",&GSLVector::get)
      .def_method("set",&GSLVector::set)
      .def_method("set_all",&GSLVector::set_all)
      .def_method("set_zero",&GSLVector::set_zero)
      .def_method("set_basis",&GSLVector::set_basis)
      .def_method("show",&GSLVector::print)
      .def_method("read",&GSLVector::read);

    Class<GSLMatrix,ArgList2<int,int> >(mod.self(),
					"matrix")
      .def_method("rows",&GSLMatrix::rows)
      .def_method("cols",&GSLMatrix::cols)
      .def_method("get",&GSLMatrix::get)
      .def_method("set",&GSLMatrix::set)
      .def_method("set_all",&GSLMatrix::set_all)
      .def_method("set_zero",&GSLMatrix::set_zero)
      .def_method("set_identity",&GSLMatrix::set_identity)
      .def_method("show",&GSLMatrix::print)
      .def_method("show_row",&GSLMatrix::print_row)
      .def_method("show_col",&GSLMatrix::print_col)
      .def_method("read",&GSLMatrix::read);;

//     Class<PeOptions>(mod.self(),"pe_options")
//       .def_attribute("cut_on",&PeOptions::cut_on)
//       .def_attribute("cut_off",&PeOptions::cut_off)
//       .def_attribute("steps",&PeOptions::steps)
//       .def_attribute("nprint",&PeOptions::nprint)
//       .def_attribute("nsave",&PeOptions::nsave)
//       .def_attribute("temp",&PeOptions::temp)
//       .def_attribute("seed",&PeOptions::seed)
//       .def_attribute("move_set",&PeOptions::move_set)
//       .def_attribute("trj",&PeOptions::trj)
//       .def_attribute("rst",&PeOptions::rst)
//       .def_attribute("cm_cut_off",&PeOptions::cm_cut_off)
//       .def_attribute("rand_temp",&PeOptions::rand_temp)
//       .def_attribute("rand_steps",&PeOptions::rand_steps)
//       .def_attribute("lambda",&PeOptions::lambda)
//       ;


//     Class<vector<PeOptions > >(mod.self(),"rex_pe_options")
//       .def_method("push_back",&vector<PeOptions>::push_back)
//       .def_method("size",&vector<PeOptions>::size);


    mod.def_function("eigenvv",
		     "Computes eigenvalues of symmetric matrix",
		     eigenvv);
    mod.def_function("cont_map",
		     "Computes the contact map of a pretein ad stores the result ina gsl::matrix.",
		     cont_map);

//     mod.def_function("mc_find",
// 		     "Well, i tries ...",
// 		     mc_find);

//     mod.def_function("mc_find_ev",
// 		     "Well, i tries ...",
// 		     mc_find_ev);


//     mod.def_function("mc_find_cont",
// 		     "Well, i tries ...",
// 		     mc_find_cont);

//     mod.def_function("mc_find_ff",
// 		     "Well, i tries ...",
// 		     mc_find_ff);

//     mod.def_function("mc_find_ff_cont",
// 		     "Well, i tries ...",
// 		     mc_find_ff_cont);

//     mod.def_function("mc_find_struct",
// 		     "Well, i tries ...",
// 		     mc_find_struct);


//     mod.def_function("rex_find_ff",
// 		     "Well, i tries ...",
// 		     rex_find_ff);

    //analysis
//     Class<gsl_energy,ArgList5<Molecules,MDB,GSLMatrix,
//       double,double> >(mod.self(),"energy")
//       .def_method("ff_energy",&gsl_energy::ff_energy)
//       .def_method("const_energy",&gsl_energy::const_energy)
//       .def_method("matches",&gsl_energy::matches)
//       .def_method("false_pos",&gsl_energy::false_pos)
//       .def_method("false_neg",&gsl_energy::false_neg);
      

  }
}