/***************************************************************************
                          mdx_mod.cpp  -  description
                             -------------------
    begin                : Tue Nov 30 17:01:41 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 <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <mdx_common.h>
#include <monitor/monitor.h>

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

ALMTOSTRING(ConstraintCollection);


ALMTOSTRING(MonitorSet<double>);
ALMTOSTRING(RMSDMonitor<double>);
ALMTOSTRING(RGYRMonitor<double>);
//ALMTOSTRING(NameSecection);

typedef MDXImpl<double,MDXOptions<double> > * 
(MDX<double>::*mdx_factory)(const Molecules & molecules,
			  const MDXOptions<double> & options,
			  const MDB & mdb);


#if (__GNUC__ >= 3)
template<>
#endif
map<int,MDX<double>::mdx_factory > MDX<double>::factory
#if (__GNUC__ >= 3)
= map<int,MDX<double>::mdx_factory >();
#endif
;

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


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

template<>
inline string to_string(const MDB & mdb){
  return "<MDB>";
}


template<>
inline string to_string(const MDX<double> &){
  return "<MDX>";
}

MDXOptions<double> c19_mdx_options(){
  return MDXOptions<double>(MDXOptions<double>::C19);
}

MDXOptions<double> c22_mdx_options(){
  return MDXOptions<double>(MDXOptions<double>::C22);
}


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

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

  return opts;
}

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


void check_xconst(ConstraintCollection & cc, 
		  XMolecules & xmols,
		  double dx, int steps){
  vector<XConstraint *> xconst;
  xconst = cc.xconstraints();
  int size = xmols.size();
  int asize= xmols[0].atom_size();
  Coor<double> **coor, **coor_buffer,**ff;
  coor = new Coor<double> * [size];
  coor_buffer = new Coor<double> *[size];
  ff = new Coor<double>* [size];
  for(int rep = 0;rep<size;rep++){
    coor[rep] = new Coor<double>(xmols[rep]);
    coor_buffer[rep] = new Coor<double>(xmols[rep]);
    ff[rep] = new Coor<double>;
    ff[rep]->resize(asize);
    ff[rep]->clear();
  }

  double W = 0;
  double Ene = 0;
  double tmp_ene = 0;
  for(int i=0;i<xconst.size();i++){
    Ene += xconst[i]->energy_force(coor,ff,size);
  }
  cout<<"Initial energy "<<Ene<<endl;
  for(int i = 0;i<steps;i++){
    for(int rep = 0;rep<size;rep++){
      coor_buffer[rep]->assign(*coor[rep]);
      ff[rep]->clear();
      for(int j = 0;j<asize;j++){
	int pos = Coor<double>::DIM*j;
	coor[rep]->coor[pos] += dx*(drand48()-0.5);
	coor[rep]->coor[pos+1] += dx*(drand48()-0.5);
	coor[rep]->coor[pos+2] += dx*(drand48()-0.5);
      }
    }
    tmp_ene = 0;
    for(int i=0;i<xconst.size();i++){
      tmp_ene += xconst[i]->energy_force(coor,ff,size);
    }
    cout<<"Current energy "<<tmp_ene<<endl;
    tmp_ene = tmp_ene - Ene;
    for(int rep = 0;rep<size;rep++){
      for(int j = 0;j<asize;j++){
	int pos = Coor<double>::DIM*j;
	W += 
	  (coor[rep]->coor[pos]-coor_buffer[rep]->coor[pos])*ff[rep]->coor[pos];
	pos = pos + 1;
	W += 
	  (coor[rep]->coor[pos]-coor_buffer[rep]->coor[pos])*ff[rep]->coor[pos];
	pos = pos + 1;
	W += 
	  (coor[rep]->coor[pos]-coor_buffer[rep]->coor[pos])*ff[rep]->coor[pos];
      }
    }
    cout<<"Current energy difference "<<tmp_ene<<" "<<W<<endl;
  }

  //cleanup
  for(int rep=0;rep<size;rep++){
    delete coor[rep];
    delete coor_buffer[rep];
    delete ff[rep];
  }
  delete [] coor;
  delete [] coor_buffer;
  delete [] ff;
}

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

    Module mod = Module("mdx","Replica molecular dynamics");
    mod.def_const("NONE",(int)Shake<double>::NONE);
    mod.def_const("BONDH",(int)Shake<double>::BONDH);
    mod.def_const("BONDALL",(int)Shake<double>::BONDALL);

    REGISTER_MDXC19_OPTIONS(double); 
    Class<MDX<double> >(mod.self(),"mdx")
      .def_method("run",&MDX<double>::run)
      .def_method("cont",&MDX<double>::cont)
      .def_method("run_from_mol",&MDX<double>::run_from_mol)
      .def_method("restart",&MDX<double>::restart)
      .def_method("set_monitors",&MDX<double>::set_monitors)
      .def_method("clear_monitors",&MDX<double>::clear_monitors);

    mod
      .def_const("C19",(int)MDXOptions<double>::C19)
      .def_const("C22",(int)MDXOptions<double>::C22)
      .def_const("NONE",(int)MDXOptions<double>::NONE)
      .def_const("SASA",(int)MDXOptions<double>::SASA)
      .def_const("EEF1",(int)MDXOptions<double>::EEF1)
      .def_const("BERENDSEN",(int)MDXOptions<double>::BERENDSEN)
      .def_const("LANGEVIN",(int)MDXOptions<double>::LANGEVIN)
      .def_const("BERE",(int)MDXOptions<double>::BERENDSEN)
      .def_const("LANG",(int)MDXOptions<double>::LANGEVIN)
      .def_const("MDX",(int)MDXOptions<double>::MDX)
      .def_const("REMD",(int)MDXOptions<double>::REMD)
      ;

    mod.def_function("const_mdx","Run a constrained MDX simulation",
		     const_mdx<MDX<double> >);


    Class<MonitorSet<double> >(mod.self(),"monitor_set")
      .def_method("size",&MonitorSet<double>::size)
      .def_method("clear",&MonitorSet<double>::clear)
      .def_method("add_rmsd",&MonitorSet<double>::add_monitor<RMSDMonitor<double> >)
      ;


    Class<RMSDMonitor<double>,ArgList2<const Molecules & , const NameSelection &> >(mod.self(),
										    "rmsd");

    mod
      .def_function("c19_options","Return options with right defaults",c19_mdx_options)
      .def_function("c22_options","Return options with right defaults",c22_mdx_options)
      .def_function("c19_defaults","Return options with right defaults",c19_mdx_defaults)
      .def_function("check_xconst","Check that MDX constraints work",check_xconst);


  }
}