/***************************************************************************
                          md_common.cpp  -  description
                             -------------------
    begin                : Sun Jan  4 11:45:27 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 <almost.h>
#include <coor.h>
#include <mdb.h>
#include <molecules/molecules.h>
#include <kernel/rgyr.h>
#include <md/md.h>
#include <sstream>
#include <ffdoc.h>
#include <rundoc.h>
using namespace Almost;

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




/* template<typename float_type>  */
/* MDOptions<float_type> c19_md_options(){ */
/*   return MDOptions<float_type>(MDOptions<float_type>::C19); */
/* } */

template<typename float_type> 
MDOptions<float_type> c22_md_options(int FF,int algo){
  MDOptions<float_type> opts(MDOptions<float_type>::C22);
  if(algo==MDOptions<float_type>::LANGEVIN){
    opts.algorithm = MDOptions<float_type>::LANGEVIN;
  }
  else
    opts.algorithm = MDOptions<float_type>::BERENDSEN;
  if(FF==MDOptions<float_type>::EEF1){
    opts.elec = MDOptions<float_type>::RDIE;
    opts. elec_cutoff = MDOptions<float_type>::SWITCH;
    opts.vdw = MDOptions<float_type>::LJ;
    opts.vdw_cutoff = MDOptions<float_type>::SWITCH;
    opts.eps = 1.0;
    opts.solv = MDOptions<float_type>::EEF1;
    opts.cut_on = 7.0;
    opts.cut_off =9.0;
    opts.cut_nb = 10.0;
  }
  else opts.solv = MDOptions<float_type>::NONE;
  return opts;
}

template<typename float_type> 
MDOptions<float_type> amber_md_options(int FF,int algo){
  MDOptions<float_type> opts(MDOptions<float_type>::AMBER);
  if(algo==MDOptions<float_type>::LANGEVIN){
    opts.algorithm = MDOptions<float_type>::LANGEVIN;
  }
  else
    opts.algorithm = MDOptions<float_type>::BERENDSEN;
  opts.solv = MDOptions<float_type>::NONE;
  return opts;
}


template<typename float_type> 
MDOptions<float_type> c19_defaults(int FF, int algo){
  MDOptions<float_type> opts;
  if(FF==MDOptions<float_type>::EEF1){
    opts.elec = MDOptions<float_type>::RDIE;
    opts. elec_cutoff = MDOptions<float_type>::SWITCH;
    opts.vdw = MDOptions<float_type>::LJ;
    opts.vdw_cutoff = MDOptions<float_type>::SWITCH;
    opts.eps = 1.0;
    opts.solv = MDOptions<float_type>::EEF1;
    opts.cut_on = 7.0;
    opts.cut_off = 9.0;
    opts.cut_nb = 10.0;
  }
  else if(FF!=MDOptions<float_type>::SASA)
    opts.solv = MDOptions<float_type>::NONE;

  if(algo==MDOptions<float_type>::LANGEVIN){
    opts.algorithm = MDOptions<float_type>::LANGEVIN;
  }
  else
    opts.algorithm = MDOptions<float_type>::BERENDSEN;

  return opts;
}

template<typename DynaType>
DynaType const_md(ConstraintCollection & cc){
  return DynaType(cc);
}


template<int i>
struct OperatorHelperElec;

template<>
struct OperatorHelperElec<0> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF,MDOptionsType>  & data){
    //    cout<<"NO ELEC"<<endl;
    return 0;
  }

  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF,SMDOptionsType>  & data){
    //    cout<<"NO ELEC"<<endl;
    return 0;
  }

};

template<>
struct OperatorHelperElec<1> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"ELEC"<<endl;
    return data.ff_data.elec_energy;
  }
  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"ELEC"<<endl;
    return data.ff_data.elec_energy;
  }

};



template<int i>
struct OperatorHelperSolv;

template<>
struct OperatorHelperSolv<0> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF,MDOptionsType>  & data){
    //    cout<<"NO SOLV"<<endl;
    return 0;
  }

  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF,SMDOptionsType>  & data){
    //    cout<<"NO SOLV"<<endl;
    return 0;
  }

};

template<>
struct OperatorHelperSolv<1> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }
  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }

};

template<>
struct OperatorHelperSolv<2> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }
  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }

};

template<>
struct OperatorHelperSolv<3> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }
  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }

};

template<>
struct OperatorHelperSolv<5> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }
  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"SOLV"<<endl;
    return data.ff_data.solv_energy;
  }

};


template<int i>
struct OperatorHelperUb;

template<>
struct OperatorHelperUb<0> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"NO UB"<<endl;
    return 0;
  }

  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"NO UB"<<endl;
    return 0;
  }
};

template<>
struct OperatorHelperUb<1> {
  template<typename FF, typename MDOptionsType>
  static typename FF::FloatType helper(MDImplData<FF, MDOptionsType>  & data){
    //    cout<<"UB"<<endl;
    return data.ff_data.ub_energy;
  }

  template<typename FF, typename SMDOptionsType>
  static typename FF::FloatType helper(SMDImplData<FF, SMDOptionsType>  & data){
    //    cout<<"UB"<<endl;
    return data.ff_data.ub_energy;
  }
};


template<typename FF, typename MDOptionsType>
DynaLogStream & operator<<(DynaLogStream &log,  
			   MDImplData<FF, MDOptionsType>  & data){
  char buff[256];
  typename FF::FloatType time = ((typename FF::FloatType)data.step)*data.time_step;
  string prmpt=Almost::DynaStatPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time","energy","rmsd","rgyr","const");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"---------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.current_energy,
	  data.rmsd,
	  data.rgyr,
	  data.ff_data.const_energy
	  );
  log.stream()<<buff;
  prmpt=Almost::InternalPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time",
	  "bond","angle","dihe","imph");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.ff_data.bond_energy,
	  data.ff_data.angle_energy,
	  data.ff_data.dihe_energy,
	  data.ff_data.imph_energy);
  log.stream()<<buff;
  prmpt=Almost::NonBondPr;
  log.stream()<<prmpt; 
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time",
	  "elec","vdw","solv","ub");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  OperatorHelperElec<FF::ELEC>::helper(data),
	  data.ff_data.vdw_energy,
	  OperatorHelperSolv<FF::SOLV>::helper(data),
	  OperatorHelperUb<FF::UB>::helper(data));
  log.stream()<<buff;
  prmpt=Almost::MDStatPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s\n",
	  "time","ekin","temp","etot");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.kinetic_energy,
	  data.current_temp,
	  data.kinetic_energy+data.current_energy);
  log.stream()<<buff;
  return log;
}




template<typename FF, typename MDOptionsType>
DynaDataStream & operator<<(DynaDataStream & dat, 
			    MDImplData<FF, MDOptionsType> & data){
  typename FF::FloatType time = ((typename FF::FloatType)data.step)*data.time_step;
  dat.stream()<<time<<"  "
	      <<data.current_energy<<" "
	      <<data.rmsd<<" "
	      <<data.rgyr<<" "
	      <<data.ff_data.const_energy<<" "
	      <<data.ff_data.bond_energy<<" "
	      <<data.ff_data.angle_energy<<" "
	      <<data.ff_data.dihe_energy<<" "
	      <<data.ff_data.imph_energy<<" "
	      <<OperatorHelperElec<FF::ELEC>::helper(data)<<" "
	      <<data.ff_data.vdw_energy<<" "
	      <<OperatorHelperSolv<FF::SOLV>::helper(data)<<" "
	      <<OperatorHelperUb<FF::UB>::helper(data)<<" "
	      <<data.kinetic_energy<<" "
	      <<data.current_temp<<" "
	      <<data.kinetic_energy+data.current_energy<<endl;
  return dat;
}

template<typename FF, typename SMDOptionsType>
DynaLogStream & operator<<(DynaLogStream &log,  
			   SMDImplData<FF, SMDOptionsType>  & data){
  char buff[256];
  typename FF::FloatType time = ((typename FF::FloatType)data.step)*data.time_step;
  string prmpt=Almost::DynaStatPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time","energy","rmsd","rgyr","const");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"---------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.current_energy,
	  data.rmsd,
	  data.rgyr,
	  data.ff_data.const_energy
	  );
  log.stream()<<buff;
  prmpt=Almost::InternalPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time",
	  "bond","angle","dihe","imph");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.ff_data.bond_energy,
	  data.ff_data.angle_energy,
	  data.ff_data.dihe_energy,
	  data.ff_data.imph_energy);
  log.stream()<<buff;
  prmpt=Almost::NonBondPr;
  log.stream()<<prmpt; 
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s | %-10s\n","time",
	  "elec","vdw","solv","ub");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  OperatorHelperElec<FF::ELEC>::helper(data),
	  data.ff_data.vdw_energy,
	  OperatorHelperSolv<FF::SOLV>::helper(data),
	  OperatorHelperUb<FF::UB>::helper(data));
  log.stream()<<buff;
  prmpt=Almost::MDStatPr;
  log.stream()<<prmpt;
  sprintf(buff,"%-10s | %-10s | %-10s | %-10s\n",
	  "time","ekin","temp","etot");
  log.stream()<<buff;
  log.stream()<<prmpt;
  log.stream()<<"--------------------------------------------------------\n";
  log.stream()<<prmpt;
  sprintf(buff,"%-10.2f   %-10.2f   %-10.2f   %-10.2f\n\n",
	  time,
	  data.kinetic_energy,
	  data.current_temp,
	  data.kinetic_energy+data.current_energy);
  log.stream()<<buff;
  return log;
}




template<typename FF, typename SMDOptionsType>
DynaDataStream & operator<<(DynaDataStream & dat, 
			    SMDImplData<FF, SMDOptionsType> & data){
  typename FF::FloatType time = ((typename FF::FloatType)data.step)*data.time_step;
  dat.stream()<<time<<"  "
	      <<data.current_energy<<" "
	      <<data.rmsd<<" "
	      <<data.rgyr<<" "
	      <<data.ff_data.const_energy<<" "
	      <<data.ff_data.bond_energy<<" "
	      <<data.ff_data.angle_energy<<" "
	      <<data.ff_data.dihe_energy<<" "
	      <<data.ff_data.imph_energy<<" "
	      <<OperatorHelperElec<FF::ELEC>::helper(data)<<" "
	      <<data.ff_data.vdw_energy<<" "
	      <<OperatorHelperSolv<FF::SOLV>::helper(data)<<" "
	      <<OperatorHelperUb<FF::UB>::helper(data)<<" "
	      <<data.kinetic_energy<<" "
	      <<data.current_temp<<" "
	      <<data.kinetic_energy+data.current_energy<<endl;
  return dat;
}

#define REGATTR(c,s,w) def_attribute(#s,&c::s,w::s)
#define REGATTR2(c,s,w) def_attribute(#s,&c::s,w)

#define REGISTER_MDC19_OPTIONS(type) \
    Class<MDOptions< type > >(mod.self(),"mdoptions") \
      .REGATTR(MDOptions< type >,kind,FFDoc)\
      .REGATTR(MDOptions< type >,elec,FFDoc)\
      .REGATTR(MDOptions< type >,elec_cutoff,FFDoc)\
      .REGATTR(MDOptions< type >,vdw,FFDoc) \
      .REGATTR(MDOptions< type >,vdw_cutoff,FFDoc) \
      .REGATTR(MDOptions< type >,solv,FFDoc) \
      .REGATTR(MDOptions< type >,solv_ignore_h,FFDoc) \
      .REGATTR(MDOptions< type >,cut_on,FFDoc) \
      .REGATTR(MDOptions< type >,cut_off,FFDoc) \
      .REGATTR(MDOptions< type >,cut_nb,FFDoc) \
      .REGATTR(MDOptions< type >,e_zero,FFDoc) \
      .REGATTR(MDOptions< type >,e_14,FFDoc) \
      .REGATTR(MDOptions< type >,eps,FFDoc) \
      .REGATTR(MDOptions< type >,r_water,FFDoc) \
      .REGATTR(MDOptions< type >,gamma_water,RunDoc) \
      .REGATTR(MDOptions< type >,alpha,FFDoc)	     \
      .REGATTR(MDOptions< type >,vect,RunDoc) \
      .REGATTR(MDOptions< type >,algorithm,RunDoc) \
      .REGATTR(MDOptions< type >,geometry,RunDoc) \
      .REGATTR(MDOptions< type >,Lx,RunDoc) \
      .REGATTR(MDOptions< type >,Ly,RunDoc) \
      .REGATTR(MDOptions< type >,Lz,RunDoc) \
      .REGATTR(MDOptions< type >,granularity,RunDoc) \
      .REGATTR(MDOptions< type >,trj,RunDoc) \
      .REGATTR(MDOptions< type >,vtrj,RunDoc)\
      .REGATTR(MDOptions< type >,rst,RunDoc)\
      .REGATTR(MDOptions< type >,log,RunDoc)\
      .REGATTR(MDOptions< type >,dat,RunDoc)\
      .REGATTR(MDOptions< type >,steps,RunDoc)\
      .REGATTR(MDOptions< type >,nprint,RunDoc)\
      .REGATTR(MDOptions< type >,nsave,RunDoc)\
      .REGATTR(MDOptions< type >,nupdate,RunDoc)\
      .REGATTR(MDOptions< type >,nrottrans,RunDoc)\
      .REGATTR(MDOptions< type >,shake_mode,RunDoc)\
      .REGATTR(MDOptions< type >,shake_tol,RunDoc)\
      .REGATTR(MDOptions< type >,temp,RunDoc)\
      .REGATTR(MDOptions< type >,time_step,RunDoc)\
      .REGATTR(MDOptions< type >,tau,RunDoc)\
      .REGATTR(MDOptions< type >,seed,RunDoc)\
      .REGATTR(MDOptions< type >,save_fitted,RunDoc)\
      .REGATTR(MDOptions< type >,selection,RunDoc)\
      .REGATTR2(MDOptions< type >,rgbmax,"Undocumented")\
      .REGATTR2(MDOptions< type >,fsmax,"Undocumented")\
      .REGATTR2(MDOptions< type >,eps_water,"Undocumented")\
      .REGATTR2(MDOptions< type >,surftens,"Undocumented")\
      .REGATTR2(MDOptions< type >,offset,"Undocumented")\
      .REGATTR2(MDOptions< type >,nsmooth,"Undocumented")\
      .REGATTR2(MDOptions< type >,memb_width,"Undocumented")\
      .REGATTR2(MDOptions< type >,elec_strength,"Undocumented");  



#define REGISTER_CONST(type)\
    mod\
      .def_const("C19",(int)MDOptions< type >::C19)\
      .def_const("C22",(int)MDOptions< type >::C22)\
      .def_const("NONE",(int)MDOptions< type >::NONE)\
      .def_const("LJ",(int)MDOptions< type >::LJ)\
      .def_const("HS",(int)MDOptions< type >::HS)\
      .def_const("SOFT",(int)MDOptions< type >::SOFT)\
      .def_const("SASA",(int)MDOptions< type >::SASA)\
      .def_const("EEF1",(int)MDOptions< type >::EEF1)\
      .def_const("IMM1",(int)MDOptions< type >::IMM1)\
      .def_const("GB",(int)MDOptions< type >::GB)\
      .def_const("GB2",(int)MDOptions< type >::GB2)\
      .def_const("BERENDSEN",(int)MDOptions< type >::BERENDSEN)\
      .def_const("LANGEVIN",(int)MDOptions< type >::LANGEVIN)\
      .def_const("BERE",(int)MDOptions< type >::BERENDSEN)\
      .def_const("LANG",(int)MDOptions< type >::LANGEVIN)\
      .def_const("OFF",(int)MDOptions< type >::OFF)\
      .def_const("RDIE",(int)MDOptions< type >::RDIE)\
      .def_const("CDIE",(int)MDOptions< type >::CDIE)\
      .def_const("LJ",(int)MDOptions< type >::LJ)\
      .def_const("SHIFT",(int)MDOptions< type >::SHIFT)\
      .def_const("SWITCH",(int)MDOptions< type >::SWITCH)\
      .def_const("TRUNC",(int)MDOptions< type >::TRUNC)\
      .def_const("PLAIN",(int)MDOptions< type >::PLAIN)\
      .def_const("CUBIC",(int)MDOptions< type >::CUBIC)\
      .def_const("NONE",(int)Shake< type >::NONE)\
      .def_const("BONDH",(int)Shake< type >::BONDH)\
      .def_const("BONDALL",(int)Shake< type >::BONDALL)\
      .def_const("ANGLEALL",(int)Shake< type >::ANGLEALL)\
      ;