/***************************************************************************
                          shx_mod.cpp  -  description
                             -------------------
    begin                : Fri Apr 22 21:26:11 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 <md/md.h>
#include <md/mdoptions.h>
#include <xmolecules.h>
#include <md_common.h>
#include <fstream>
#include <iterator>
#include <stdlib.h>
#include <sstream>
#include <forcefield/const/base/dihedralconst.h>
#include <forcefield/const/misc/hbconst.h>
#include <forcefield/const/pmf/apmf.h>
#include <kernel/geometry.h>
#include <ff.h>
#include <algorithm/sd_fwd.h>
using namespace Almost;

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

ALMTOSTRING(Molecules);
ALMTOSTRING(MDB);

extern void mol2pdb(const Molecules & mols, string ofile);
extern Coor<double> mols_coor(const Molecules & molecules);
int alm_shx_main(const Molecules & molecules,
		 double * ha,
		 double * h,
		 double * n,
		 double * ca,
		 double * cb,
		 double * c
		 );

int alm_shx_main_fake(const Molecules & molecules,
		      double * ha,
		      double * h,
		      double * n,
		      double * ca,
		      double * cb,
		      double * c
		 );


static string apmf_dir = "../PMF/";

template<typename float_type>
struct SHXOptions {
  //force field options
  enum { NONE = 0, SASA = 1, EEF1 = 2};
  enum { BERENDSEN = 0, LANGEVIN = 1, NEWTON = 2}; 
  enum { C19 = 0, C22 = 1};
  enum { OFF, RDIE, CDIE, LJ, SHIFT, SWITCH };
  enum { PLAIN, CUBIC};
  int kind;
  int elec;
  int elec_cutoff;
  int vdw;
  int vdw_cutoff;

  int  solv;
  bool solv_ignore_h;
  float_type cut_on;
  float_type cut_off;
  float_type cut_nb;
  float_type e_zero;
  float_type e_14;
  float_type eps;
  float_type r_water;
  float_type gamma_water; //NEW
  bool vect;

  //MD
  int algorithm;

  //Geometry
  int geometry;
  float_type Lx;
  float_type Ly;
  float_type Lz;
  int  granularity;

  string trj;
  string vtrj;
  string rst;
  string log;
  string dat;
  int steps;
  int nprint;
  int nsave;
  int nupdate;
  int nrottrans;
  int nsurf; //NEW
  int shake_mode;
  float_type shake_tol;
  float_type begin_temp;
  float_type end_temp;
  int sa_steps;
  int ga_steps;
  float_type time_step;
  float_type tau;
  int seed;

  bool save_fitted;
  string selection;
  //MIP
  int nreplica;
  SHXOptions(int k = C19):kind(k){
    if(kind==C19) c19_defaults();
    else if(kind==C22) c22_defaults();
    else c19_defaults();
    r_water = 1.4;
    gamma_water = 91;
    algorithm = BERENDSEN;
    geometry = PLAIN;
    Lx = Ly = Lz = 3*cut_nb;
    granularity = 1;      
    trj = "run_coor.dcd";
    vtrj = "none";
    rst = "run.rst";
    log = "stdout";
    dat = "none";
    steps = 100000;
    nprint = 1000;
    nsave = 1000;
    nupdate = 10;
    nsurf = 1000;
    nrottrans = 1000;
    shake_mode = Shake<float_type>::BONDH;
    shake_tol = 10E-6;
    begin_temp = 1000.00;
    end_temp = 0.00;
    sa_steps = 100;
    ga_steps = 100;
    time_step = 0.002;
    tau = 5.0;
    seed = 30869;
    save_fitted=0;
    selection = "///CA";
    nreplica = 2;
  }

private:
  void c19_defaults(){
    elec = RDIE;
    elec_cutoff = SHIFT;
    vdw = LJ;
    vdw_cutoff = SHIFT;
    solv = SASA;
    solv_ignore_h =1;
    cut_on = 6.5;
    cut_off = 7.5;
    cut_nb = 8.0;
    e_zero = 332.0716;
    e_14 = 0.4;
    eps = 2.0;
    vect = 0;
  }

  void c22_defaults(){
    elec = RDIE;
    elec_cutoff = SHIFT;
    vdw = LJ;
    vdw_cutoff = SHIFT;
    solv = NONE;
    solv_ignore_h =0;
    cut_on = 10;
    cut_off = 12;
    cut_nb = 14;
    e_zero = 332.0716;
    e_14 = 1;
    eps = 1.0;
    vect = 0;
  }
};
  
template<typename float_type>
inline ostream & operator<<(ostream & out, 
			    const SHXOptions<float_type> & opt){
  out<<"\t>> Force field kind:       ";
  if(opt.kind == SHXOptions<float_type>::C19)
    out<<"     C19\n";
  else
    out<<"     C22\n";
  out<<"\t>> Electrostatic:          ";
  if(opt.elec==SHXOptions<float_type>::RDIE)
    out<<"     RDIE"<<"\n";
  else if(opt.elec==SHXOptions<float_type>::CDIE)
    out<<"     CDIE"<<"\n";
  else
    out<<"     Off"<<"\n";
  out<<"\t>> Elec-Cutoff:            ";
  if(opt.elec_cutoff==SHXOptions<float_type>::SHIFT)
    out<<"     SHIFT"<<"\n";
  else
    out<<"     SWITCH"<<"\n";
  out<<"\t>> VdW-Cutoff:             ";
  if(opt.vdw_cutoff==SHXOptions<float_type>::SHIFT)
    out<<"     SHIFT"<<"\n";
  else
    out<<"     SWITCH"<<"\n";
  out<<"\t>> Solvation:              ";
  if(opt.solv==SHXOptions<float_type>::SASA)
    out<<"     SASA"<<"\n";
  else if(opt.solv==SHXOptions<float_type>::EEF1)
    out<<"     EEF1"<<"\n";
  else
    out<<"     none"<<"\n";
  out<<"\t>> Cut on:                      "
     <<opt.cut_on<<"\n";
  out<<"\t>> Cut off:                     "
     <<opt.cut_off<<"\n";
  out<<"\t>> Cut nb:                      "
     <<opt.cut_nb<<"\n";
  out<<"\t>> E_zero:                      "
     <<opt.e_zero<<"\n";
  out<<"\t>> E_14:                        "
     <<opt.e_14<<"\n";
  out<<"\t>> Eps:                         "
     <<opt.eps<<"\n";
  ;
  out<<"\t>> r_water:                     "<<opt.r_water<<"\n";
  out<<"\t>> Vector code:                 ";
  if(opt.vect)
    out<<"yes"<<"\n";
  else
    out<<"no"<<"\n";
  out<<"\t>> Algorithm:                   ";
  if(opt.algorithm==SHXOptions<float_type>::BERENDSEN)
    out<<"Berendsen\n";
  else if(opt.algorithm==SHXOptions<float_type>::LANGEVIN)
    out<<"Langevin\n";
  else
    out<<"Newton\n";
  out<<"\t>> Geometry:                    ";
  if(opt.geometry==SHXOptions<float_type>::PLAIN)
    out<<"plain\n";
  else
    out<<"cubic";
  out<<"\t>> Box size:                    "
     <<opt.Lx<<"x"<<opt.Ly<<"x"<<opt.Lz<<"\n";
  out<<"\t>> Granularity:                 "
     <<opt.granularity<<"\n";
  out<<"\t>> Trajectory:                  "
     <<opt.trj<<"\n";
  out<<"\t>> Velocity:                    "
     <<opt.vtrj<<"\n";
  out<<"\t>> Restart:                     "
     <<opt.rst<<"\n";
  out<<"\t>> Log    :                     "
     <<opt.log<<"\n";
  out<<"\t>> Data   :                     "
     <<opt.dat<<"\n";
  out<<"\t>> Number of steps:             "
     <<opt.steps<<"\n";
  out<<"\t>> Print frequency:             "
     <<opt.nprint<<"\n";
  out<<"\t>> Save frequency:              "
     <<opt.nsave<<"\n";
  out<<"\t>> NB Update frequency:         "
     <<opt.nupdate<<"\n";
  out<<"\t>> RT stopping frequency:       "
     <<opt.nrottrans<<"\n";
  out<<"\t>> SHAKE mode:                  ";
  switch(opt.shake_mode){
  case Shake<float_type>::NONE:
    out<<"none"<<"\n";
    break;
  case Shake<float_type>::BONDH:
    out<<"bond-H"<<"\n";
    break;
  case Shake<float_type>::BONDALL:
    out<<"all bonds"<<"\n";
    break;
  default:
    out<<"unknown (error)"<<"\n";
  }
  out<<"\t>> SHAKE tolerance:             "<<opt.shake_tol<<"\n";
  out<<"\t>> Begin temperature:           "
     <<opt.begin_temp<<"\n";
  out<<"\t>> End temperature:             "
     <<opt.end_temp<<"\n";
  out<<"\t>> Simulated annealing steps:   "
     <<opt.sa_steps<<"\n";
  out<<"\t>> Genetic algorithm steps:     "
     <<opt.ga_steps<<"\n";
  out<<"\t>> Time step:                   "<<opt.time_step<<"\n";
  out<<"\t>> Tau:                         "<<opt.tau<<"\n";
  out<<"\t>> Seed:                        "
     <<opt.seed<<"\n";
  out<<"\t>> Selection:                   "
     <<opt.selection<<"\n";
  out<<"\t>> Save fitted:                 "
     <<opt.save_fitted<<"\n";
  out<<"\t>> Number  of replicas:         "<<opt.nreplica<<"\n";
  return out;
}

template<typename float_type>
struct SHXMCOptions {
  //force field options
  enum { NONE = 0, SASA = 1, EEF1 = 2};
  enum { BERENDSEN = 0, LANGEVIN = 1, NEWTON = 2}; 
  enum { C19 = 0, C22 = 1};
  enum { OFF, RDIE, CDIE, LJ, SHIFT, SWITCH };
  enum { PLAIN, CUBIC};
  int kind;
  int elec;
  int elec_cutoff;
  int vdw;
  int vdw_cutoff;

  int  solv;
  bool solv_ignore_h;
  float_type cut_on;
  float_type cut_off;
  float_type cut_nb;
  float_type e_zero;
  float_type e_14;
  float_type eps;
  float_type r_water;
  float_type gamma_water; //NEW
  bool vect;

  //MD
  int algorithm;

  //Geometry
  int geometry;
  float_type Lx;
  float_type Ly;
  float_type Lz;
  int  granularity;

  string trj;
  string vtrj;
  string rst;
  string log;
  string dat;
  int steps;
  int nprint;
  int nsave;
  int nupdate;
  int nrottrans;
  int nsurf; //NEW
  int shake_mode;
  float_type shake_tol;
  float_type begin_temp;
  float_type end_temp;
  int sa_steps;
  int mc_steps;
  float_type mc_temp;
  float_type time_step;
  float_type tau;
  int seed;

  bool save_fitted;
  string selection;



  //MIP
  SHXMCOptions(int k = C19):kind(k){
    if(kind==C19) c19_defaults();
    else if(kind==C22) c22_defaults();
    else c19_defaults();
    r_water = 1.4;
    gamma_water = 91;
    algorithm = BERENDSEN;
    geometry = PLAIN;
    Lx = Ly = Lz = 3*cut_nb;
    granularity = 1;      
    trj = "run_coor.dcd";
    vtrj = "none";
    rst = "run.rst";
    log = "stdout";
    dat = "none";
    steps = 100000;
    nprint = 1000;
    nsave = 1000;
    nupdate = 10;
    nsurf = 1000;
    nrottrans = 1000;
    shake_mode = Shake<float_type>::BONDH;
    shake_tol = 10E-6;
    begin_temp = 1000.00;
    end_temp = 0.00;
    sa_steps = 100;
    mc_temp = 1;
    mc_steps = 100;
    time_step = 0.002;
    tau = 5.0;
    seed = 30869;
    save_fitted=0;
    selection = "///CA";
  }

private:
  void c19_defaults(){
    elec = RDIE;
    elec_cutoff = SHIFT;
    vdw = LJ;
    vdw_cutoff = SHIFT;
    solv = SASA;
    solv_ignore_h =1;
    cut_on = 6.5;
    cut_off = 7.5;
    cut_nb = 8.0;
    e_zero = 332.0716;
    e_14 = 0.4;
    eps = 2.0;
    vect = 0;
  }

  void c22_defaults(){
    elec = RDIE;
    elec_cutoff = SHIFT;
    vdw = LJ;
    vdw_cutoff = SHIFT;
    solv = NONE;
    solv_ignore_h =0;
    cut_on = 10;
    cut_off = 12;
    cut_nb = 14;
    e_zero = 332.0716;
    e_14 = 1;
    eps = 1.0;
    vect = 0;
  }
};
  
template<typename float_type>
inline ostream & operator<<(ostream & out, 
			    const SHXMCOptions<float_type> & opt){
  out<<"\t>> Force field kind:       ";
  if(opt.kind == SHXMCOptions<float_type>::C19)
    out<<"     C19\n";
  else
    out<<"     C22\n";
  out<<"\t>> Electrostatic:          ";
  if(opt.elec==SHXMCOptions<float_type>::RDIE)
    out<<"     RDIE"<<"\n";
  else if(opt.elec==SHXMCOptions<float_type>::CDIE)
    out<<"     CDIE"<<"\n";
  else
    out<<"     Off"<<"\n";
  out<<"\t>> Elec-Cutoff:            ";
  if(opt.elec_cutoff==SHXMCOptions<float_type>::SHIFT)
    out<<"     SHIFT"<<"\n";
  else
    out<<"     SWITCH"<<"\n";
  out<<"\t>> VdW-Cutoff:             ";
  if(opt.vdw_cutoff==SHXMCOptions<float_type>::SHIFT)
    out<<"     SHIFT"<<"\n";
  else
    out<<"     SWITCH"<<"\n";
  out<<"\t>> Solvation:              ";
  if(opt.solv==SHXMCOptions<float_type>::SASA)
    out<<"     SASA"<<"\n";
  else if(opt.solv==SHXMCOptions<float_type>::EEF1)
    out<<"     EEF1"<<"\n";
  else
    out<<"     none"<<"\n";
  out<<"\t>> Cut on:                      "
     <<opt.cut_on<<"\n";
  out<<"\t>> Cut off:                     "
     <<opt.cut_off<<"\n";
  out<<"\t>> Cut nb:                      "
     <<opt.cut_nb<<"\n";
  out<<"\t>> E_zero:                      "
     <<opt.e_zero<<"\n";
  out<<"\t>> E_14:                        "
     <<opt.e_14<<"\n";
  out<<"\t>> Eps:                         "
     <<opt.eps<<"\n";
  ;
  out<<"\t>> r_water:                     "<<opt.r_water<<"\n";
  out<<"\t>> Vector code:                 ";
  if(opt.vect)
    out<<"yes"<<"\n";
  else
    out<<"no"<<"\n";
  out<<"\t>> Algorithm:                   ";
  if(opt.algorithm==SHXMCOptions<float_type>::BERENDSEN)
    out<<"Berendsen\n";
  else if(opt.algorithm==SHXMCOptions<float_type>::LANGEVIN)
    out<<"Langevin\n";
  else
    out<<"Newton\n";
  out<<"\t>> Geometry:                    ";
  if(opt.geometry==SHXMCOptions<float_type>::PLAIN)
    out<<"plain\n";
  else
    out<<"cubic";
  out<<"\t>> Box size:                    "
     <<opt.Lx<<"x"<<opt.Ly<<"x"<<opt.Lz<<"\n";
  out<<"\t>> Granularity:                 "
     <<opt.granularity<<"\n";
  out<<"\t>> Trajectory:                  "
     <<opt.trj<<"\n";
  out<<"\t>> Velocity:                    "
     <<opt.vtrj<<"\n";
  out<<"\t>> Restart:                     "
     <<opt.rst<<"\n";
  out<<"\t>> Log    :                     "
     <<opt.log<<"\n";
  out<<"\t>> Data   :                     "
     <<opt.dat<<"\n";
  out<<"\t>> Number of steps:             "
     <<opt.steps<<"\n";
  out<<"\t>> Print frequency:             "
     <<opt.nprint<<"\n";
  out<<"\t>> Save frequency:              "
     <<opt.nsave<<"\n";
  out<<"\t>> NB Update frequency:         "
     <<opt.nupdate<<"\n";
  out<<"\t>> RT stopping frequency:       "
     <<opt.nrottrans<<"\n";
  out<<"\t>> SHAKE mode:                  ";
  switch(opt.shake_mode){
  case Shake<float_type>::NONE:
    out<<"none"<<"\n";
    break;
  case Shake<float_type>::BONDH:
    out<<"bond-H"<<"\n";
    break;
  case Shake<float_type>::BONDALL:
    out<<"all bonds"<<"\n";
    break;
  default:
    out<<"unknown (error)"<<"\n";
  }
  out<<"\t>> SHAKE tolerance:             "<<opt.shake_tol<<"\n";
  out<<"\t>> Begin temperature:           "
     <<opt.begin_temp<<"\n";
  out<<"\t>> End temperature:             "
     <<opt.end_temp<<"\n";
  out<<"\t>> Simulated annealing steps:   "
     <<opt.sa_steps<<"\n";
  out<<"\t>> Monte Carlo temperature:     "
     <<opt.mc_temp<<"\n";
  out<<"\t>> Monte Carlo steps:           "
     <<opt.mc_steps<<"\n";
  out<<"\t>> Time step:                   "<<opt.time_step<<"\n";
  out<<"\t>> Tau:                         "<<opt.tau<<"\n";
  out<<"\t>> Seed:                        "
     <<opt.seed<<"\n";
  out<<"\t>> Selection:                   "
     <<opt.selection<<"\n";
  out<<"\t>> Save fitted:                 "
     <<opt.save_fitted<<"\n";
  return out;
}


class GAConstraints {
  
public:
  struct ResidueConst {
    ResidueConst(){
      phi = psi = omega = 999;
    }

    double phi;
    double psi;
    double omega;
  };
  
  vector< vector<ResidueConst> > const_seeds;
  vector< vector<ResidueConst> > current_const;
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  const Molecules * molecules_;

  GAConstraints(const Molecules & molecules, string file, string shx_file,
		int ng){
    ifstream in;
    in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: no such file "<<file<<endl;
      exit(1);
    }

    parse(molecules,in);
    //Set up current const;
    current_const.resize(ng);
    for(int rep = 0; rep < ng ; rep++)
      current_const[rep].resize(molecules.fragment_size()); 
    ResidueConst rc;
    for(int rep = 0; rep < ng ; rep++){
      for(int i=0;i<const_seeds.size();i++){
	if(const_seeds[i].size()>0){
	  int r = lrand48();
	  r = r%const_seeds[i].size();
	  current_const[rep][i] = const_seeds[i][r];
	}
	else 	current_const[rep][i] = rc;
      }
    }
    molecules_ = &molecules;
    read_shifts(shx_file,HA,HN,N,CA,CB,CO);
    k[0] = 17.53;
    k[1] = 10;
    k[2] = 0.18;
    k[3] = 0.99;
    k[4] = 0.91;
    k[5] = 1.21;
  }


  DihedralConst<double> constraints(int rep, int step){

    DihedralConst<double> dc(*molecules_);
    char buff[256];
    sprintf(buff,"const_%i_%i.cst",rep,step);
    ofstream dbg;
    dbg.open(buff);
    for(int i=0;i<current_const[rep].size();i++){
      int prev = i;
      int curr = i+1;
      int next = i+2;
      sprintf(buff,"//*_%i/",prev);
      string p = string(buff);

      sprintf(buff,"//*_%i/",curr);
      string c = string(buff);

      sprintf(buff,"//*_%i/",next);
      string n = string(buff);

      if(current_const[rep][i].phi!=999){
	dc.assign(p+"C",
		  c+"N",
		  c+"CA",
		  c+"C",
		  current_const[rep][i].phi,
		  10+step);
      }

      if(current_const[rep][i].psi!=999){
	dc.assign(c+"N",
		  c+"CA",
		  c+"C",
		  n+"N",
		  current_const[rep][i].psi,
		  10+step);
      }

      if(current_const[rep][i].omega!=999){
	dc.assign(c+"CA",
		  c+"C",
		  n+"N",
		  n+"CA",
		  current_const[rep][i].omega,
		  20);
      }
      dbg <<i+1<<" "<<current_const[rep][i].phi<<" "
	  <<current_const[rep][i].psi<<" "
	  <<current_const[rep][i].omega<<"\n";
    }

    return dc;
  }


  void eval(const XMolecules & xmols, const vector<double> & vdw){
    map<double, int> score;

    for(int rep = 0; rep< xmols.size();rep++){
      vector<double> ha,hn,n,ca,cb,co;
      double e = 0;
      mol2pdb(xmols[rep],"__.pdb");
      int status = 
	system("shifts_local __.pdb");
      if(status==-1){
	cout<<"FATAL ERROR: could't compute shifts\n";
	exit(-1);
      }
      else {
	read_shifts("__.shx",
		    ha,hn,n,ca,cb,co
		    );
      }

//       for(int i= 0;i<ha.size();i++){
// 	e += k[0]*(ha[i]-HA[i])*(ha[i]-HA[i]);
// 	e += k[1]*(hn[i]-HN[i])*(hn[i]-HN[i]);
// 	e += k[2]*(n[i]-N[i])*(n[i]-N[i]);
// 	e += k[3]*(ca[i]-CA[i])*(ca[i]-CA[i]);
// 	e += k[4]*(cb[i]-CB[i])*(cb[i]-CB[i]);
// 	//	e += k[5]*(co[i]-CO[i])*(co[i]-CO[i]);
//       }
//       e += vdw[rep];
      double c;
      c = correlation(ha,HA);
      cout<<"\t>> Replica "<<rep+1<<" HA correlation "<<c<<"\n";
      e = k[0]*(1-c);

      c = correlation(hn,HN);
      cout<<"\t>> Replica "<<rep+1<<" HN correlation "<<c<<"\n";
      //e += k[1]*(1-c);

      c = correlation(n,N);
      cout<<"\t>> Replica "<<rep+1<<" N  correlation "<<c<<"\n";
      e += k[2]*(1-c);

      c = correlation(ca,CA);
      cout<<"\t>> Replica "<<rep+1<<" CA correlation "<<c<<"\n";
      e += k[3]*(1-c);

      c = correlation(cb,CB);
      cout<<"\t>> Replica "<<rep+1<<" CB correlation "<<c<<"\n";
      e += k[4]*(1-c);

      e *= 100;
//       e += vdw[rep];
      while(score.find(e)!=score.end()) e += 0.0001;
      score[e] = rep;
      cout<<"\t>> Score replica "<<rep+1<<" "<<e<<endl;
    }
    
    // evolve
    vector<vector<ResidueConst> > new_const;
    new_const.resize(xmols.size());
    //Take First randomized
    map<double,int>::iterator iter = score.begin(),end = score.end();
    new_const[0] = current_const[iter->second]; ++iter;
    new_const[1] = current_const[iter->second]; ++iter;
    double rand = drand48();
    if(rand<0.5){
      cross(new_const[0],new_const[1]);
    }
    else {
      shake(new_const[0]);
      shake(new_const[1]);
    }
    //For the oder take other candidate
    ResidueConst rc;
    for(int rep =2; rep<xmols.size();rep++){
      rand = drand48();
      if(rand<0.5){
	new_const[rep] = current_const[iter->second];
	if(rand<0.25)
	  shake(new_const[rep]);
	else{
	  cross_noshake(new_const[rep-1],new_const[rep]);
	}	  
      }
      else {
	for(int i=0;i<const_seeds.size();i++){
	  if(const_seeds[i].size()>0){
	    int r = lrand48();
	    r = r%const_seeds[i].size();
	    new_const[rep].push_back(const_seeds[i][r]);
	  }
	  else 	new_const[rep].push_back(rc);
	}
      }
      ++iter;
    }

    iter = score.begin();
    int pos = 0;
    while(iter!=end){
      current_const[iter->second]  = new_const[pos];
      ++iter; ++pos;
    }
  }
private:

  double correlation(vector<double> & v1, vector<double> & v2){
    double X,Y,XX,YY,XY;
    X = Y = XX = YY = XY = 0;
    for(int i=0;i<v1.size();i++){
      XY += v1[i]*v2[i];
      X  += v1[i];
      XX += v1[i]*v1[i];
      Y  += v2[i];
      YY += v2[i]*v2[i];
    }
    int n= v1.size();

    return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
  }
  void shake(vector<ResidueConst> & rc){
    
    for(int r = 0;r<rc.size();r++){
      double d = drand48();
      if(d<0.75) continue;
      if(rc[r].phi!=999){
	rc[r].phi = rc[r].phi +10.0*(drand48()-0.5);
      }
      if(rc[r].psi!=999){
	rc[r].psi = rc[r].psi +10.0*(drand48()-0.5);
      }
//       if(rc[r].omega!=999){
// 	rc[r].omega = rc[r].omega +10.0*(drand48()-0.5);
//       }
    }
  }

  void cross(vector<ResidueConst> & rc1,
	     vector<ResidueConst> & rc2){
    int r = lrand48()%rc1.size();
    for(int i=0;i<r;i++){
      ResidueConst tmp = rc1[i];
      //rc1[i] = rc2[i];
      rc2[i] = tmp;
    }
    shake(rc1);
  }

  void cross_noshake(vector<ResidueConst> & rc1,
		     vector<ResidueConst> & rc2){
    int r = lrand48()%rc1.size();
    for(int i=0;i<r;i++){
      ResidueConst tmp = rc1[i];
      //rc1[i] = rc2[i];
      rc2[i] = tmp;
    }
  }

  void parse(const Molecules & molecules, ifstream & in){
    const_seeds.resize(molecules.fragment_size());
    istream_iterator<string> iter(in),end;
    //Format num phi psi omega
    while(iter!=end){
      int num = atoi((*iter).c_str()); ++iter;
      if((num<1)||(num>molecules.fragment_size())){
	cout<<"FATAL ERROR: invalid fragment number "<<num<<endl;
	exit(-1);
      }
      --num;
      double phi = atof((*iter).c_str()); ++iter;
      double psi = atof((*iter).c_str()); ++iter;
      double omega = atof((*iter).c_str()); ++iter;
      ResidueConst rc;
      rc.phi = phi;
      rc.psi = psi;
      rc.omega = omega;
      const_seeds[num].push_back(rc);
    }
  }

  void read_shifts(string file,
		   vector<double> &ha_,
		   vector<double> &hn_,
		   vector<double> &n_,
		   vector<double> &ca_,
		   vector<double> &cb_,
		   vector<double> &co_){
    ifstream in; in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: couldn't open shift file "<<file<<endl;
      exit(-1);
    }
    istream_iterator<string> iter(in),end;
    while(iter!=end){
      double sh;
      sh = atof((*iter).c_str());++iter;
      ha_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      hn_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      n_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      ca_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      cb_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      co_.push_back(sh);
    }
  }
};


class MCConstraints {
  
public:
  struct ResidueConst {
    ResidueConst(){
      phi = psi = omega = 999;
    }

    double phi;
    double psi;
    double omega;
  };
  
  vector< vector<ResidueConst> > const_seeds;
  vector<ResidueConst> current_const;
  vector<ResidueConst> buffer_const;
  double current_energy;
  double tmp_energy;
  double temp;
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  const Molecules * molecules_;
  APMF<double> pmf;
  MCConstraints(const Molecules & molecules, string file, string shx_file)
  :pmf(molecules,apmf_dir+"frag.lib",apmf_dir+"alm.pmf",
       APMF<double>::ALL,0.00181242){
    ifstream in;
    in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: no such file "<<file<<endl;
      exit(1);
    }

    parse(molecules,in);
    //Set up current const;
    current_const.resize(molecules.fragment_size()); 
    ResidueConst rc;
    for(int i=0;i<const_seeds.size();i++){
      if(const_seeds[i].size()>0){
	int r = lrand48();
	r = r%const_seeds[i].size();
	current_const[i] = const_seeds[i][r];
      }
      else 	current_const[i] = rc;
    }

    buffer_const = current_const;
    molecules_ = &molecules;
    read_shifts(shx_file,HA,HN,N,CA,CB,CO);
    k[0] = 17.53;
    k[1] = 10;
    k[2] = 0.18;
    k[3] = 0.99;
    k[4] = 0.91;
    k[5] = 1.21;

    eval(molecules);
    current_energy = tmp_energy;
    current_const = buffer_const;
  }

  void set_temp(double t) { temp = t;}
  DihedralConst<double> constraints(int step, int mc_step){

    DihedralConst<double> dc(*molecules_);
    char buff[256];
    sprintf(buff,"const_%i.cst",mc_step);
    ofstream dbg;
    dbg.open(buff);
    for(int i=0;i<current_const.size();i++){
      int prev = i;
      int curr = i+1;
      int next = i+2;
      sprintf(buff,"//*_%i/",prev);
      string p = string(buff);

      sprintf(buff,"//*_%i/",curr);
      string c = string(buff);

      sprintf(buff,"//*_%i/",next);
      string n = string(buff);

      if(current_const[i].phi!=999){
	dc.assign(p+"C",
		  c+"N",
		  c+"CA",
		  c+"C",
		  current_const[i].phi,
		  10+2*step);
      }

      if(current_const[i].psi!=999){
	dc.assign(c+"N",
		  c+"CA",
		  c+"C",
		  n+"N",
		  current_const[i].psi,
		  10+2*step);
      }

      if(current_const[i].omega!=999){
	dc.assign(c+"CA",
		  c+"C",
		  n+"N",
		  n+"CA",
		  current_const[i].omega,
		  10);
      }
      dbg <<i+1<<" "<<current_const[i].phi<<" "
	  <<current_const[i].psi<<" "
	  <<current_const[i].omega<<"\n";
    }

    return dc;
  }


  void eval(const Molecules & mols){
    vector<double> ha,hn,n,ca,cb,co;
    double e = 0;
    mol2pdb(mols,"__.pdb");
    int status = 
      system("shifts_local __.pdb");
    if(status==-1){
      cout<<"FATAL ERROR: could't compute shifts\n";
      exit(-1);
    }
    else {
      read_shifts("__.shx",
		  ha,hn,n,ca,cb,co
		  );
    }

    double c;
    c = correlation(ha,HA);
    cout<<"\t>> HA correlation "<<c<<"\n";
    e = k[0]*(1-c);

    c = correlation(hn,HN);
    cout<<"\t>> HN correlation "<<c<<"\n";
    //    e += k[1]*(1-c);

    c = correlation(n,N);
    cout<<"\t>> N   correlation "<<c<<"\n";
    //    e += k[2]*(1-c);

    c = correlation(ca,CA);
    cout<<"\t>> CA correlation "<<c<<"\n";
    //    e += k[3]*(1-c);

    c = correlation(cb,CB);
    cout<<"\t>> CB correlation "<<c<<"\n";
    //    e += k[4]*(1-c);

    e *= 100;
    double pmf_ene = pmf.energy(mols_coor(mols));
    cout<<"\t>> Coorelation score: "<<e<<"\n";
    cout<<"\t>> PMF Energy: "<<pmf_ene<<"\n";
    tmp_energy = pmf_ene/(log(1+0.01*e));
    {
      //set current to measured
      vector<int> n;
      vector<int> c;
      vector<int> ca;

      for(int i=0;i<mols.atom_size();i++){
	if(mols[i].name()=="N") n.push_back(i);
	if(mols[i].name()=="C") c.push_back(i);
	if(mols[i].name()=="CA") ca.push_back(i);
      }
      Coor<double> coor(mols);
      const int dim = Coor<double>::DIM;
      current_const[0].psi = 
	(180.0/M_PI)*Kernel::Geometry::dihedral(coor.coor+dim*n[0],
						coor.coor+dim*ca[0],
						coor.coor+dim*c[0],
						coor.coor+dim*n[1]);
					      

      for(int aa = 1; aa<current_const.size()-1;aa++){
	//
	current_const[aa].phi = 
	  (180.0/M_PI)*Kernel::Geometry::dihedral(coor.coor+dim*c[aa-1],
						  coor.coor+dim*n[aa],
						  coor.coor+dim*ca[aa],
						  coor.coor+dim*c[aa]);

	current_const[aa].psi = 
	  (180.0/M_PI)*Kernel::Geometry::dihedral(coor.coor+dim*n[aa],
						  coor.coor+dim*ca[aa],
						  coor.coor+dim*c[aa],
						  coor.coor+dim*n[aa+1]);

      }
      int aa = current_const.size()-1;
      current_const[aa].phi = 
	(180.0/M_PI)*Kernel::Geometry::dihedral(coor.coor+dim*c[aa-1],
						coor.coor+dim*n[aa],
						coor.coor+dim*ca[aa],
						coor.coor+dim*c[aa]);
      
    }
  }


  void move(){
    //MODIF: shake now 50%
    double d = drand48();
    if(d<0.5){
      shake(current_const);
      cout<<"\t>> SHXMC shake const\n";
    }
    else {
      rebuild(current_const);
      cout<<"\t>> SHXMC build const\n";
    }
  }


  bool metropolis(){
    if(tmp_energy<current_energy){
      current_energy = tmp_energy;
      buffer_const = current_const;
      cout<<"\t>> Score "<<current_energy<<endl;
      cout<<"\t>> Accepted\n";
      return true;
    }
    double d = drand48();
    if(d<exp(-(tmp_energy-current_energy)/temp)){
      current_energy = tmp_energy;
      buffer_const = current_const;
      cout<<"\t>> Score "<<current_energy<<endl;
      cout<<"\t>> Accepted\n";
      return true;
    }
    else {
      current_const = buffer_const;
      cout<<"\t>> Score "<<tmp_energy<<" "<<"("<<current_energy<<")"<<endl;
      cout<<"\t>> Rejected\n";
      return false;


    }
  }


private:

  double correlation(vector<double> & v1, vector<double> & v2){
    double X,Y,XX,YY,XY;
    X = Y = XX = YY = XY = 0;
    for(int i=0;i<v1.size();i++){
      XY += v1[i]*v2[i];
      X  += v1[i];
      XX += v1[i]*v1[i];
      Y  += v2[i];
      YY += v2[i]*v2[i];
    }
    int n= v1.size();

    return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
  }
  void shake(vector<ResidueConst> & rc){
    //MODIF: shake delta now 25 deg
    for(int r = 0;r<rc.size();r++){
      double d = drand48();
      if(d<0.75) continue;
      if(rc[r].phi!=999){
	rc[r].phi = rc[r].phi +25.0*(drand48()-0.5);
      }
      if(rc[r].psi!=999){
	rc[r].psi = rc[r].psi +25.0*(drand48()-0.5);
      }
//       if(rc[r].omega!=999){
// 	rc[r].omega = rc[r].omega +10.0*(drand48()-0.5);
//       }
    }
  }


  void rebuild(vector<ResidueConst> & rest){
    ResidueConst rc;
    for(int i=0;i<const_seeds.size();i++){
      if(const_seeds[i].size()>0){
	int r = lrand48();
	r = r%const_seeds[i].size();
	rest[i] = const_seeds[i][r];
      }
      else 	rest[i] = rc;
    }
  }


  void parse(const Molecules & molecules, ifstream & in){
    const_seeds.resize(molecules.fragment_size());
    istream_iterator<string> iter(in),end;
    //Format num phi psi omega
    while(iter!=end){
      int num = atoi((*iter).c_str()); ++iter;
      if((num<1)||(num>molecules.fragment_size())){
	cout<<"FATAL ERROR: invalid fragment number!"<<endl;
	exit(-1);
      }
      --num;
      double phi = atof((*iter).c_str()); ++iter;
      double psi = atof((*iter).c_str()); ++iter;
      double omega = atof((*iter).c_str()); ++iter;
      ResidueConst rc;
      rc.phi = phi;
      rc.psi = psi;
      rc.omega = omega;
      const_seeds[num].push_back(rc);
    }
  }

  void read_shifts(string file,
		   vector<double> &ha_,
		   vector<double> &hn_,
		   vector<double> &n_,
		   vector<double> &ca_,
		   vector<double> &cb_,
		   vector<double> &co_){
    ifstream in; in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: couldn't open shift file "<<file<<endl;
      exit(-1);
    }
    istream_iterator<string> iter(in),end;
    while(iter!=end){
      double sh;
      sh = atof((*iter).c_str());++iter;
      ha_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      hn_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      n_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      ca_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      cb_.push_back(sh);

      sh = atof((*iter).c_str());++iter;
      co_.push_back(sh);
    }
  }
};

 class SHX {
  GAConstraints ga_const;
  vector<double> vdw_energy;
public:
  SHX(const Molecules & molecules, string file, string shx_file, int ng)
    :ga_const(molecules,file,shx_file,ng){
    vdw_energy.resize(ng);
  }

  XMolecules run(const XMolecules & molecules, 
		 const SHXOptions<double> & options,
		 const MDB & mdb){
    srand48(options.seed);
    //Does plain aneeling
    XMolecules m = molecules;
    for(int ga_step = 0; ga_step<options.ga_steps; ga_step++){
      for(int i=0;i<molecules.size();i++){
	for(int a=0;a<molecules[i].atom_size();a++){
	  m[i][a][0] = molecules[i][a][0];
	  m[i][a][1] = molecules[i][a][1];
	  m[i][a][2] = molecules[i][a][2];
	}
      }
      for(int i=0;i<molecules.size();i++){	
	Molecules tmp = sa(m[i],options,mdb,i,ga_step);
	for(int a=0;a<tmp.atom_size();a++){
	  m[i][a][0] = tmp[a][0];
	  m[i][a][1] = tmp[a][1];
	  m[i][a][2] = tmp[a][2];
	}
      }
      ga_const.eval(m,vdw_energy);
    }
  }

private:
  Molecules sa(const Molecules & m, const SHXOptions<double> & options,
	       const MDB & mdb, int rep, int ga_step){
    double b = options.begin_temp;
    double e = options.end_temp;
    double dt = (b-e)/options.sa_steps;
    Molecules res = m;
    for(int a=0;a<m.atom_size();a++){
      res[a][0] = m[a][0];
      res[a][1] = m[a][1];
      res[a][2] = m[a][2];
    }    
    for(int step=0;step<options.sa_steps;step++){
      ConstraintCollection cc;
      DihedralConst<double> dc = ga_const.constraints(rep,step);
      HBConst<double> hb = HBConst<double>(m);
      cc.add(dc);
      cc.add(hb);
      MD<double> md(cc);
      MDOptions<double> opts = getopts(options,rep,step,ga_step,b-step*dt);
      Molecules ress = md.run(res,opts,mdb);
      for(int a=0;a<ress.atom_size();a++){
	res[a][0] = ress[a][0];
	res[a][1] = ress[a][1];
	res[a][2] = ress[a][2];
      }
    }
    EnergyOptions<double> ene_opts = geteneopts(options);
    ConstraintCollection cc;
    DihedralConst<double> dc = ga_const.constraints(rep,options.sa_steps);
//     APMF<double> pmf = APMF<double>(m,"frag.lib","alm.pmf",
// 				    APMF<double>::ALL,
// 				    0.001);
    cc.add(dc);
//     cc.add(pmf);
    EnergyData<double> ene = const_sd(res,ene_opts,cc,mdb,1000,0.01);
    vdw_energy[rep] = ene["const"];
    char buff[256];
    sprintf(buff,"coor_pmf_%i_%i.pdb",rep,ga_step);
    mol2pdb(res,buff);
    return res;
  }

  MDOptions<double> getopts(const SHXOptions<double> & options,
			    int rep, int step, int ga_step,
			    double temp){
    MDOptions<double> opts;
    opts.kind = options.kind;
    opts.elec = options.elec;
    opts.elec_cutoff = options.elec_cutoff;
    opts.vdw = options.vdw;
    opts.vdw_cutoff = options.vdw_cutoff;
    opts.solv = options.solv;
    opts.solv_ignore_h = options.solv_ignore_h;
    opts.cut_on = options.cut_on;
    opts.cut_off = options.cut_off;
    opts.cut_nb = options.cut_nb;
    opts.e_zero = options.e_zero;
    opts.e_14 = options.e_14;
    opts.eps = options.eps;
    opts.r_water = options.r_water;
    opts.gamma_water = options.gamma_water;
    opts.vect = options.vect;
    opts.algorithm = options.algorithm;
    opts.geometry = options.geometry;
    opts.Lx = options.Lx;
    opts.Ly = options.Ly;
    opts.Lz = options.Lz;
    opts.granularity = options.granularity;
    opts.trj = mkname(rep,step,ga_step,options.trj);
    //    opts.vtrj = mkname(rep,step,ga_step,options.vtrj);
    opts.rst = mkname(rep,step,ga_step,options.rst);
    //    opts.log = mkname(rep,step,ga_step,options.log);
    opts.dat = mkname(rep,step,ga_step,options.dat);
    opts.steps = options.steps;
    opts.nprint = options.nprint;
    opts.nsave = options.nsave;
    opts.nupdate = options.nupdate;
    opts.nrottrans = options.nrottrans;
    opts.nsurf = options.nsurf;
    opts.shake_mode = options.shake_mode;
    opts.shake_tol = options.shake_tol;
    opts.temp = temp;
    opts.time_step = options.time_step;
    opts.tau = options.tau;
    opts.seed = options.seed+rep*step;
    opts.save_fitted = options.save_fitted;
    opts.selection = options.selection;
    return opts;
  }


  EnergyOptions<double> geteneopts(const SHXOptions<double> & options){
    EnergyOptions<double> opts;
    opts.kind = options.kind;
    opts.elec = options.elec;
    opts.elec_cutoff = options.elec_cutoff;
    opts.vdw = options.vdw;
    opts.vdw_cutoff = options.vdw_cutoff;
    opts.solv = options.solv;
    opts.solv_ignore_h = options.solv_ignore_h;
    opts.cut_on = options.cut_on;
    opts.cut_off = options.cut_off;
    opts.cut_nb = options.cut_nb;
    opts.e_zero = options.e_zero;
    opts.e_14 = options.e_14;
    opts.eps = options.eps;
    opts.r_water = options.r_water;
    opts.vect = options.vect;
    opts.geometry = options.geometry;
    opts.Lx = options.Lx;
    opts.Ly = options.Ly;
    opts.Lz = options.Lz;
    opts.granularity = options.granularity;
    return opts;
  }
  string mkname(int rep,int step, int ga_step, string s){
    char buff[256];
    sprintf(buff,"_%i_%i_%i",ga_step,rep,step);
    int pos = s.find(".");
    string name;
    if(pos!=string::npos){
      name = string(s.begin(),s.begin()+pos)
	+string(buff)+
	string(s.begin()+pos,s.end());
    }
    else {
      name = s+string(buff);
    }
    cout<<"DEBUG name "<<name<<endl;
    return name;
  }
};


class SHXMC {
  MCConstraints mc_const;
public:
  SHXMC(const Molecules & molecules, string file, string shx_file)
    :mc_const(molecules,file,shx_file){}

  XMolecules run(const Molecules & molecules, 
		 const SHXMCOptions<double> & options,
		 const MDB & mdb){
    srand48(options.seed);

    mc_const.set_temp(options.mc_temp);
    Molecules m = molecules;
    Molecules b = molecules;
    for(int a=0;a<molecules.atom_size();a++){
      m[a][0] = molecules[a][0];
      m[a][1] = molecules[a][1];
      m[a][2] = molecules[a][2];
      b[a][0] = molecules[a][0];
      b[a][1] = molecules[a][1];
      b[a][2] = molecules[a][2];
    }

    for(int i=0;i<options.mc_steps;i++){
      mc_const.move();
      Molecules tmp = sa(m,options,mdb,i);
      for(int a=0;a<tmp.atom_size();a++){
	m[a][0] = tmp[a][0];
	m[a][1] = tmp[a][1];
	m[a][2] = tmp[a][2];
      }
      mc_const.eval(m);
      if(mc_const.metropolis()){
	for(int a=0;a<tmp.atom_size();a++){
	  b[a][0] = m[a][0];
	  b[a][1] = m[a][1];
	  b[a][2] = m[a][2];
	}
      }
      else {
	for(int a=0;a<tmp.atom_size();a++){
	  m[a][0] = b[a][0];
	  m[a][1] = b[a][1];
	  m[a][2] = b[a][2];
	}
      }
    }
  }

private:
  Molecules sa(const Molecules & m, const SHXMCOptions<double> & options,
	       const MDB & mdb,int mc_step){
    double b = options.begin_temp;
    double e = options.end_temp;
    double dt = (b-e)/options.sa_steps;
    Molecules res = m;
    for(int a=0;a<m.atom_size();a++){
      res[a][0] = m[a][0];
      res[a][1] = m[a][1];
      res[a][2] = m[a][2];
    }    
    for(int step=0;step<options.sa_steps;step++){
      ConstraintCollection cc;
      DihedralConst<double> dc = mc_const.constraints(step,mc_step);
      HBConst<double> hb = HBConst<double>(m);
      cc.add(dc);
      cc.add(hb);
      MD<double> md(cc);
      MDOptions<double> opts = getopts(options,step,mc_step,b-step*dt);
      Molecules ress = md.run(res,opts,mdb);
      for(int a=0;a<ress.atom_size();a++){
	res[a][0] = ress[a][0];
	res[a][1] = ress[a][1];
	res[a][2] = ress[a][2];
      }
    }
    EnergyOptions<double> ene_opts = geteneopts(options);
    ConstraintCollection cc;
    DihedralConst<double> dc = mc_const.constraints(options.sa_steps,mc_step);
    cc.add(dc);
    //MODIF setps where 10000
    EnergyData<double> ene = const_sd(res,ene_opts,cc,mdb,1000,0.01);
//     vdw_energy[rep] = ene["vdw"];
    char buff[256];
    sprintf(buff,"coor_%i.pdb",mc_step);
    mol2pdb(res,buff);
    return res;
  }

  MDOptions<double> getopts(const SHXMCOptions<double> & options,
			    int step, int mc_step,
			    double temp){
    MDOptions<double> opts;
    opts.kind = options.kind;
    opts.elec = options.elec;
    opts.elec_cutoff = options.elec_cutoff;
    opts.vdw = options.vdw;
    opts.vdw_cutoff = options.vdw_cutoff;
    opts.solv = options.solv;
    opts.solv_ignore_h = options.solv_ignore_h;
    opts.cut_on = options.cut_on;
    opts.cut_off = options.cut_off;
    opts.cut_nb = options.cut_nb;
    opts.e_zero = options.e_zero;
    opts.e_14 = options.e_14;
    opts.eps = options.eps;
    opts.r_water = options.r_water;
    opts.gamma_water = options.gamma_water;
    opts.vect = options.vect;
    opts.algorithm = options.algorithm;
    opts.geometry = options.geometry;
    opts.Lx = options.Lx;
    opts.Ly = options.Ly;
    opts.Lz = options.Lz;
    opts.granularity = options.granularity;
    opts.trj = mkname(step,mc_step,options.trj);
    //    opts.vtrj = mkname(step,mc_step,options.vtrj);
    opts.rst = mkname(step,mc_step,options.rst);
    //    opts.log = mkname(step,mc_step,options.log);
    opts.dat = mkname(step,mc_step,options.dat);
    opts.steps = options.steps;
    opts.nprint = options.nprint;
    opts.nsave = options.nsave;
    opts.nupdate = options.nupdate;
    opts.nrottrans = options.nrottrans;
    opts.nsurf = options.nsurf;
    opts.shake_mode = options.shake_mode;
    opts.shake_tol = options.shake_tol;
    opts.temp = temp;
    opts.time_step = options.time_step;
    opts.tau = options.tau;
    opts.seed = options.seed+step;
    opts.save_fitted = options.save_fitted;
    opts.selection = options.selection;
    return opts;
  }


  EnergyOptions<double> geteneopts(const SHXMCOptions<double> & options){
    EnergyOptions<double> opts;
    opts.kind = options.kind;
    opts.elec = options.elec;
    opts.elec_cutoff = options.elec_cutoff;
    opts.vdw = options.vdw;
    opts.vdw_cutoff = options.vdw_cutoff;
    opts.solv = options.solv;
    opts.solv_ignore_h = options.solv_ignore_h;
    opts.cut_on = options.cut_on;
    opts.cut_off = options.cut_off;
    opts.cut_nb = options.cut_nb;
    opts.e_zero = options.e_zero;
    opts.e_14 = options.e_14;
    opts.eps = options.eps;
    opts.r_water = options.r_water;
    opts.vect = options.vect;
    opts.geometry = options.geometry;
    opts.Lx = options.Lx;
    opts.Ly = options.Ly;
    opts.Lz = options.Lz;
    opts.granularity = options.granularity;
    return opts;
  }
  string mkname(int step, int mc_step, string s){
    char buff[256];
    sprintf(buff,"_%i_%i",mc_step,step);
    int pos = s.find(".");
    string name;
    if(pos!=string::npos){
      name = string(s.begin(),s.begin()+pos)
	+string(buff)+
	string(s.begin()+pos,s.end());
    }
    else {
      name = s+string(buff);
    }
    cout<<"DEBUG name "<<name<<endl;
    return name;
  }
};


void read_shifts(string file,
		 vector<double> &ha_,
		 vector<double> &hn_,
		 vector<double> &n_,
		 vector<double> &ca_,
		 vector<double> &cb_,
		 vector<double> &co_){
  ifstream in; in.open(file.c_str());
  if(!in){
    cout<<"FATAL ERROR: couldn't open shift file "<<file<<endl;
    exit(-1);
  }
  istream_iterator<string> iter(in),end;
  while(iter!=end){
    double sh;
    sh = atof((*iter).c_str());++iter;
    ha_.push_back(sh);

    sh = atof((*iter).c_str());++iter;
    hn_.push_back(sh);

    sh = atof((*iter).c_str());++iter;
    n_.push_back(sh);

    sh = atof((*iter).c_str());++iter;
    ca_.push_back(sh);

    sh = atof((*iter).c_str());++iter;
    cb_.push_back(sh);

    sh = atof((*iter).c_str());++iter;
    co_.push_back(sh);
  }
}

double correlation(vector<double> & v1, vector<double> & v2){
  double X,Y,XX,YY,XY;
  X = Y = XX = YY = XY = 0;
  for(int i=0;i<v1.size();i++){
    XY += v1[i]*v2[i];
    X  += v1[i];
    XX += v1[i]*v1[i];
    Y  += v2[i];
    YY += v2[i]*v2[i];
  }
  int n= v1.size();
  
  return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
}


double correlation(double * v1, vector<double> & v2){
  double X,Y,XX,YY,XY;
  X = Y = XX = YY = XY = 0;
  for(int i=0;i<v2.size();i++){
    XY += v1[i]*v2[i];
    X  += v1[i];
    XX += v1[i]*v1[i];
    Y  += v2[i];
    YY += v2[i]*v2[i];
  }
  int n= v2.size();
  
  return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
}


double newcorrelation(vector<double> & v1, vector<double> & v2){
  double X,Y,XX,YY,XY,Z;
  X = Y = XX = YY = XY = Z = 0;
  for(int i=0;i<v1.size();i++){
    if(v1[i]!=0){
    if(v2[i]!=0){
    XY += v1[i]*v2[i];
    X  += v1[i];
    XX += v1[i]*v1[i];
    Y  += v2[i];
    YY += v2[i]*v2[i]; 
    Z +=1;
    } 
   }
  }
  double n= Z;

  return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
}


double newcorrelation(double * v1, vector<double> & v2){
  double X,Y,XX,YY,XY,Z;
  X = Y = XX = YY = XY = Z = 0;
  for(int i=0;i<v2.size();i++){
    if(v1[i]!=0){
    if(v2[i]!=0){
    XY += v1[i]*v2[i];
    X  += v1[i];
    XX += v1[i]*v1[i];
    Y  += v2[i];
    YY += v2[i]*v2[i];
    Z+=1;  
    }
   }
  }
  double n= Z;

  return (n*XY-X*Y)/sqrt((n*XX-X*X)*(n*YY-Y*Y));
}



double score(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  k[0] = 17.53;
  k[1] = 10;
  k[2] = 0.18;
  k[3] = 0.99;
  k[4] = 0.91;
  k[5] = 1.21;
  read_shifts(file,HA,HN,N,CA,CB,CO);
  vector<double> ha,hn,n,ca,cb,co;
  mol2pdb(mols,"__.pdb");
  //  int status = 
  system("shifts_local __.pdb");

  read_shifts("__.shx",
	      ha,hn,n,ca,cb,co
	      );


  {
    int size = HA.size();
    for(int i=0;i<size;i++){
      if(HA[i]==0) ha[i] =0;
      if(HN[i]==0) hn[i] =0;
      if(N[i]==0) n[i] =0;
      if(CA[i]==0) ca[i] =0;
      if(CB[i]==0) cb[i] =0;
      if(CO[i]==0) co[i] =0;
    }
  }

  double c;
  double e = 0;
  c = correlation(ha,HA);
  cout<<"HA correlation "<<c<<"\n";
  e = k[0]*(1-c);
  
  c = correlation(hn,HN);
  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);
  
  c = correlation(n,N);
  cout<<"N  correlation "<<c<<"\n";
  e += k[2]*(1-c);

  c = correlation(ca,CA);
  cout<<"CA correlation "<<c<<"\n";
  e += k[3]*(1-c);

  c = correlation(cb,CB);
  cout<<"CB correlation "<<c<<"\n";
  e += k[4]*(1-c);

  e *= 100;
  return e;
}



double score_silent(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  k[0] = 17.53;
  k[1] = 10;
  k[2] = 0.18;
  k[3] = 0.99;
  k[4] = 0.91;
  k[5] = 1.21;
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  read_shifts(file,HA,HN,N,CA,CB,CO);
  alm_shx_main(mols,ha,hn,n,ca,cb,co);
  {
    int size = HA.size();
    for(int i=0;i<size;i++){
      if(HA[i]==0) ha[i] =0;
      if(HN[i]==0) hn[i] =0;
      if(N[i]==0) n[i] =0;
      if(CA[i]==0) ca[i] =0;
      if(CB[i]==0) cb[i] =0;
      if(CO[i]==0) co[i] =0;
    }
  }

  double c;
  double e = 0;
  c = correlation(ha,HA);
  cout<<"HA correlation "<<c<<"\n";
  e = k[0]*(1-c);
  
  c = correlation(hn,HN);
  //  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);
  
  c = correlation(n,N);
  cout<<"N  correlation "<<c<<"\n";
  e += k[2]*(1-c);

  c = correlation(ca,CA);
  cout<<"CA correlation "<<c<<"\n";
  e += k[3]*(1-c);

  c = correlation(cb,CB);
  cout<<"CB correlation "<<c<<"\n";
  e += k[4]*(1-c);

  e *= 100;
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
  //  cout<<"DEBUG shx "<<e<<endl;
  return e;
}

double new_score_silent(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  k[0] = 17.53;
  k[1] = 10;
  k[2] = 10;
  k[3] = 10;
  k[4] = 10;
  k[5] = 10;
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  read_shifts(file,HA,HN,N,CA,CB,CO);
  alm_shx_main(mols,ha,hn,n,ca,cb,co);
  {
    int size = HA.size();
    for(int i=0;i<size;i++){
      if(HA[i]==0) ha[i] =0;
      if(HN[i]==0) hn[i] =0;
      if(N[i]==0) n[i] =0;
      if(CA[i]==0) ca[i] =0;
      if(CB[i]==0) cb[i] =0;
      if(CO[i]==0) co[i] =0;
    }
  }

  double c;
  double e = 0;
  c = newcorrelation(ha,HA);
  cout<<"HA correlation "<<c<<"\n";
  if(c>.911){
  c=.911;
  }
  e = 0;//k[0]*(1-c);

  c = newcorrelation(hn,HN);
  //  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);

  c = newcorrelation(n,N);
  cout<<"N  correlation "<<c<<"\n";
  if(c>.909){
  c=.909;
  }

  e += k[2]*(1-c);

  c = newcorrelation(ca,CA);
  cout<<"CA correlation "<<c<<"\n";
  if(c>.980){
  c=.980;
  }

  e += k[3]*(1-c);

  c = newcorrelation(cb,CB);
  cout<<"CB correlation "<<c<<"\n";
  if(c>.996){
  c=.996;
  }

  e += k[4]*(1-c);

  e *= 100;
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
  //  cout<<"DEBUG shx "<<e<<endl;
  return e;
}



void compute_shx(const  Molecules & mols,
		 vector<double> & HA,
		 vector<double> & H,
		 vector<double> & N,
		 vector<double> & CA,
		 vector<double> & CB,
		 vector<double> & C
		 ){
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  alm_shx_main(mols,ha,hn,n,ca,cb,co);
  for(int i=0;i<fsize;i++){
    HA.push_back(ha[i]);
    H.push_back(hn[i]);
    N.push_back(n[i]);
    CA.push_back(ca[i]);
    CB.push_back(cb[i]);
    C.push_back(co[i]);
  }
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
}

double score_noha(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  k[0] = 17.53;
  k[1] = 10;
  k[2] = 0.18;
  k[3] = 0.99;
  k[4] = 0.91;
  k[5] = 1.21;
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  read_shifts(file,HA,HN,N,CA,CB,CO);
  alm_shx_main(mols,ha,hn,n,ca,cb,co);
  {
    int size = HA.size();
    for(int i=0;i<size;i++){
      if(HA[i]==0) ha[i] =0;
      if(HN[i]==0) hn[i] =0;
      if(N[i]==0) n[i] =0;
      if(CA[i]==0) ca[i] =0;
      if(CB[i]==0) cb[i] =0;
      if(CO[i]==0) co[i] =0;
    }
  }

  double c;
  double e = 0;
  c = correlation(ha,HA);
  cout<<"HA correlation "<<c<<"\n";
  //e = k[0]*(1-c);
  
  c = correlation(hn,HN);
  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);
  
  c = correlation(n,N);
  cout<<"N  correlation "<<c<<"\n";
  e += k[2]*(1-c);

  c = correlation(ca,CA);
  cout<<"CA correlation "<<c<<"\n";
  e += k[3]*(1-c);

  c = correlation(cb,CB);
  cout<<"CB correlation "<<c<<"\n";
  e += k[4]*(1-c);

  e *= 100;
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
  //  cout<<"DEBUG shx "<<e<<endl;
  return e;
}


double score_silent_fib(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  //  k[0] = 17.53;
  //  k[1] = 10;
  k[2] = 0.18;
  k[3] = 0.99;
  k[4] = 0.91;
  k[5] = 1.21;
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  read_shifts(file,HA,HN,N,CA,CB,CO);
  alm_shx_main_fake(mols,ha,hn,n,ca,cb,co);
  {
    int size = HA.size();
    for(int i=0;i<size;i++){
      if(HA[i]==0) ha[i] =0;
      if(HN[i]==0) hn[i] =0;
      if(N[i]==0) n[i] =0;
      if(CA[i]==0) ca[i] =0;
      if(CB[i]==0) cb[i] =0;
      if(CO[i]==0) co[i] =0;
    }
  }

  double c;
  double e = 0;
  //  c = correlation(ha,HA);
  //  cout<<"HA correlation "<<c<<"\n";
  //  e = k[0]*(1-c);
  
  //  c = correlation(hn,HN);
  //  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);
  
  c = correlation(n,N);
  cout<<"N  correlation "<<c<<"\n";
  e = k[2]*(1-c);

  c = correlation(ca,CA);
  cout<<"CA correlation "<<c<<"\n";
  e += k[3]*(1-c);

  c = correlation(cb,CB);
  cout<<"CB correlation "<<c<<"\n";
  e += k[4]*(1-c);

  e *= 100;
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
  //  cout<<"DEBUG shx "<<e<<endl;
  return e;
}


double score_silent_fib2(const Molecules & mols, string file){
  vector<double> HA,HN,N,CA,CB,CO;
  double k[6];
  //  k[0] = 17.53;
  //  k[1] = 10;
  k[2] = 0.18;
  k[3] = 0.99;
  k[4] = 0.91;
  k[5] = 1.21;
  double *ha,*hn,*n,*ca,*cb,*co;
  int fsize = mols.fragment_size();
  ha = new double[fsize];
  hn = new double[fsize];
  n = new double[fsize];
  ca = new double[fsize];
  cb = new double[fsize];
  co = new double[fsize];
  read_shifts(file,HA,HN,N,CA,CB,CO);
  alm_shx_main_fake(mols,ha,hn,n,ca,cb,co);
  int size = HA.size();
  for(int i=0;i<size;i++){
    if(HA[i]==0) ha[i] =0;
    if(HN[i]==0) hn[i] =0;
    if(N[i]==0) n[i] =0;
    if(CA[i]==0) ca[i] =0;
    if(CB[i]==0) cb[i] =0;
    if(CO[i]==0) co[i] =0;
  }


  //  double c;
  double e = 0;
  //  c = correlation(ha,HA);
  //  cout<<"HA correlation "<<c<<"\n";
  //  e = k[0]*(1-c);
  
  //  c = correlation(hn,HN);
  //  cout<<"HN correlation "<<c<<"\n";
  //e += k[1]*(1-c);
  double c1=0,c2=0,c3=0;
  for(int i=0;i<size;i++){
//     cout<<HA[i]<<" "<<ha[i]<<" "
// 	<<N[i]<<" "<<n[i]<<" "
// 	<<CA[i]<<" "<<ca[i]<<endl;
    c1 += k[2]*(N[i]-n[i])*(N[i]-n[i]);
    c2 += k[3]*(CA[i]-ca[i])*(CA[i]-ca[i]);
    c3 += k[4]*(CB[i]-cb[i])*(CB[i]-cb[i]);
  }
  cout<<"N  score "<<c1<<"\n";
  cout<<"CA score "<<c2<<"\n";
  cout<<"CB score "<<c3<<"\n";
  e = c1+c2+c3;
  delete [] ha;
  delete [] hn;
  delete [] n;
  delete [] ca;
  delete [] cb;
  delete [] co;
  //  cout<<"DEBUG shx "<<e<<endl;
  return e;
}



template<typename float_type>
inline string to_string(const SHXOptions<float_type> & c19md_opt){
  stringstream s;
  s<<c19md_opt;
  return s.str();
}

template<typename float_type>
inline string to_string(const SHXMCOptions<float_type> & c19md_opt){
  stringstream s;
  s<<c19md_opt;
  return s.str();
}


//template<typename float_type>
inline string to_string(const SHX &){
  return "<SHX>";
}

//template<typename float_type>
inline string to_string(const SHXMC &){
  return "<SHXMC>";
}

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

extern "C" {
  int    shx_main( int argc, char *argv[],
		   double * ha,
		   double * h,
		   double * n,
		   double * ca,
		   double * cb,
		   double * c
		   );
}


int alm_shx_main(const Molecules & molecules,
		 double * ha,
		 double * h,
		 double * n,
		 double * ca,
		 double * cb,
		 double * c
		 ){
  //Print molecules to tmp
  int pid = getpid();
  //  char buff[256];
  //  sprintf(buff,"/tmp/%i.pdb",pid);


  int argc=4;

  char ** argv = new char*[4];
  for(int i=0;i<4;i++)
    argv[i] = new char[256];
  
  strcpy(argv[0],"almost");
  strcpy(argv[1],"1");
  sprintf(argv[2],"/tmp/%i.pdb",pid);
  sprintf(argv[3],"/tmp/%i.s",pid);
  mol2pdb(molecules,argv[2]);
  shx_main(argc,argv,ha,h,n,ca,cb,c);
  unlink(argv[2]);
  unlink(argv[3]);
  for(int i=0;i<4;i++) delete [] argv[i];
  delete [] argv;
  return 1;
}


void mol2pdb_fake(const Molecules &molecules,char * file){
  ofstream out; out.open(file);
  char buff[180];
  for(int i=0;i<molecules.atom_size();i++){
    int fragd = molecules.find_fragment(i);
    string name = molecules[i].name();
    string frag_name = molecules.fragment_name(fragd,Molecules::SHORT);
    sprintf(buff,
	    "%-6s%5i  %-3s%4s %c %3i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
	    "ATOM",
	    i+1,
	    name.c_str(),
	    frag_name.c_str(),
	    ' ',
	    fragd+1,
	    molecules[i].x(),
	    molecules[i].y(),
	    molecules[i].z(),
	    (double)1,
	    (double)0,
	    "MOL"
	    );
    out<<buff;
  }
  out<<"END"<<endl;
}

int alm_shx_main_fake(const Molecules & molecules,
		      double * ha,
		      double * h,
		      double * n,
		      double * ca,
		      double * cb,
		      double * c
		 ){
  //Print molecules to tmp
  int pid = getpid();
  //  char buff[256];
  //  sprintf(buff,"/tmp/%i.pdb",pid);


  int argc=4;

  char ** argv = new char*[4];
  for(int i=0;i<4;i++)
    argv[i] = new char[256];
  
  strcpy(argv[0],"almost");
  strcpy(argv[1],"1");
  sprintf(argv[2],"/tmp/%i.pdb",pid);
  sprintf(argv[3],"/tmp/%i.s",pid);
  mol2pdb_fake(molecules,argv[2]);
  shx_main(argc,argv,ha,h,n,ca,cb,c);
  unlink(argv[2]);
  unlink(argv[3]);
  for(int i=0;i<4;i++) delete [] argv[i];
  delete [] argv;
  return 1;
}

//#include "Main.c"

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

    Module mod = Module("shx","All about chemical shifts");
    Class<SHX,ArgList4<const Molecules &, string, string, int> >(mod.self(),"shx")
      .def_method("run",&SHX::run);

    Class<SHXMC,ArgList3<const Molecules &, string, string> >(mod.self(),"shx_mc")
      .def_method("run",&SHXMC::run);

#define REGATTR(c,s,w) def_attribute(#s,&c::s,w::s)
    
    Class<SHXOptions<double> >(mod.self(),"gaoptions") 
      .REGATTR(SHXOptions<double>,kind,FFDoc)
      .REGATTR(SHXOptions<double>,elec,FFDoc)
      .REGATTR(SHXOptions<double>,elec_cutoff,FFDoc)
      .REGATTR(SHXOptions<double>,vdw,FFDoc) 
      .REGATTR(SHXOptions<double>,vdw_cutoff,FFDoc) 
      .REGATTR(SHXOptions<double>,solv,FFDoc) 
      .REGATTR(SHXOptions<double>,solv_ignore_h,FFDoc) 
      .REGATTR(SHXOptions<double>,cut_on,FFDoc) 
      .REGATTR(SHXOptions<double>,cut_off,FFDoc) 
      .REGATTR(SHXOptions<double>,cut_nb,FFDoc) 
      .REGATTR(SHXOptions<double>,e_zero,FFDoc) 
      .REGATTR(SHXOptions<double>,e_14,FFDoc) 
      .REGATTR(SHXOptions<double>,eps,FFDoc) 
      .REGATTR(SHXOptions<double>,r_water,FFDoc) 
      .REGATTR(SHXOptions<double>,gamma_water,RunDoc) 
      .REGATTR(SHXOptions<double>,vect,RunDoc) 
      .REGATTR(SHXOptions<double>,algorithm,RunDoc) 
      .REGATTR(SHXOptions<double>,geometry,RunDoc) 
      .REGATTR(SHXOptions<double>,Lx,RunDoc) 
      .REGATTR(SHXOptions<double>,Ly,RunDoc) 
      .REGATTR(SHXOptions<double>,Lz,RunDoc) 
      .REGATTR(SHXOptions<double>,granularity,RunDoc) 
      .REGATTR(SHXOptions<double>,trj,RunDoc) 
      .REGATTR(SHXOptions<double>,vtrj,RunDoc)
      .REGATTR(SHXOptions<double>,rst,RunDoc)
      .REGATTR(SHXOptions<double>,log,RunDoc)
      .REGATTR(SHXOptions<double>,dat,RunDoc)
      .REGATTR(SHXOptions<double>,steps,RunDoc)
      .REGATTR(SHXOptions<double>,nprint,RunDoc)
      .REGATTR(SHXOptions<double>,nsave,RunDoc)
      .REGATTR(SHXOptions<double>,nupdate,RunDoc)
      .REGATTR(SHXOptions<double>,nrottrans,RunDoc)
      .REGATTR(SHXOptions<double>,shake_mode,RunDoc)
      .REGATTR(SHXOptions<double>,shake_tol,RunDoc)
      //      .REGATTR(SHXOptions<double>,temp,RunDoc)
      .REGATTR(SHXOptions<double>,time_step,RunDoc)
      .REGATTR(SHXOptions<double>,tau,RunDoc)
      .REGATTR(SHXOptions<double>,seed,RunDoc)
      .REGATTR(SHXOptions<double>,save_fitted,RunDoc)
      .REGATTR(SHXOptions<double>,selection,RunDoc)
      .REGATTR(SHXOptions<double>,nreplica,RunDoc); 

    Class<SHXMCOptions<double> >(mod.self(),"mcoptions")
      .REGATTR(SHXMCOptions<double>,kind,FFDoc)
      .REGATTR(SHXMCOptions<double>,elec,FFDoc)
      .REGATTR(SHXMCOptions<double>,elec_cutoff,FFDoc)
      .REGATTR(SHXMCOptions<double>,vdw,FFDoc) 
      .REGATTR(SHXMCOptions<double>,vdw_cutoff,FFDoc) 
      .REGATTR(SHXMCOptions<double>,solv,FFDoc) 
      .REGATTR(SHXMCOptions<double>,solv_ignore_h,FFDoc) 
      .REGATTR(SHXMCOptions<double>,cut_on,FFDoc) 
      .REGATTR(SHXMCOptions<double>,cut_off,FFDoc) 
      .REGATTR(SHXMCOptions<double>,cut_nb,FFDoc) 
      .REGATTR(SHXMCOptions<double>,e_zero,FFDoc) 
      .REGATTR(SHXMCOptions<double>,e_14,FFDoc) 
      .REGATTR(SHXMCOptions<double>,eps,FFDoc) 
      .REGATTR(SHXMCOptions<double>,r_water,FFDoc) 
      .REGATTR(SHXMCOptions<double>,gamma_water,RunDoc) 
      .REGATTR(SHXMCOptions<double>,vect,RunDoc) 
      .REGATTR(SHXMCOptions<double>,algorithm,RunDoc) 
      .REGATTR(SHXMCOptions<double>,geometry,RunDoc) 
      .REGATTR(SHXMCOptions<double>,Lx,RunDoc) 
      .REGATTR(SHXMCOptions<double>,Ly,RunDoc) 
      .REGATTR(SHXMCOptions<double>,Lz,RunDoc) 
      .REGATTR(SHXMCOptions<double>,granularity,RunDoc) 
      .REGATTR(SHXMCOptions<double>,trj,RunDoc) 
      .REGATTR(SHXMCOptions<double>,vtrj,RunDoc)
      .REGATTR(SHXMCOptions<double>,rst,RunDoc)
      .REGATTR(SHXMCOptions<double>,log,RunDoc)
      .REGATTR(SHXMCOptions<double>,dat,RunDoc)
      .REGATTR(SHXMCOptions<double>,steps,RunDoc)
      .REGATTR(SHXMCOptions<double>,nprint,RunDoc)
      .REGATTR(SHXMCOptions<double>,nsave,RunDoc)
      .REGATTR(SHXMCOptions<double>,nupdate,RunDoc)
      .REGATTR(SHXMCOptions<double>,nrottrans,RunDoc)
      .REGATTR(SHXMCOptions<double>,shake_mode,RunDoc)
      .REGATTR(SHXMCOptions<double>,shake_tol,RunDoc)
      //      .REGATTR(SHXMCOptions<double>,temp,RunDoc)
      .REGATTR(SHXMCOptions<double>,time_step,RunDoc)
      .REGATTR(SHXMCOptions<double>,tau,RunDoc)
      .REGATTR(SHXMCOptions<double>,seed,RunDoc)
      .REGATTR(SHXMCOptions<double>,save_fitted,RunDoc)
      .REGATTR(SHXMCOptions<double>,selection,RunDoc)
      ;  

    mod
      .def_function("score","Returns shx score",score)
      .def_function("score_silent","Returns shx score",score_silent)
      .def_function("new_score_silent","Returns shx score",new_score_silent)
      .def_function("score_noha","Returns shx score",score_noha)
      .def_function("score_silent_fix","Returns shx score",score_silent_fib)
      .def_function("score_silent_fix2","Returns shx score",score_silent_fib2)
      //      .def_function("shx_score","Returns shx score",alm_shx_main)
      .def_const("apmf_dir",apmf_dir);
      ;
  }

}