#include <string>
#include <iostream>
#include <boost/algorithm/string.hpp>
#include <mkff.h>
using namespace boost;
using namespace std;

string file_banner =
  "/***************************************************************************\n"
  "                          __FILE__  -  description\n"
  "                             -------------------\n"
  "    begin                : Thu Jan  6 10:20:31 2005\n"
  "    copyright            : (C) 2009 by Cavalli Andrea\n"
  "    author               : Andrea Cavalli\n"
  "    date                 : $Date$\n"
  "    id                   : $Id$\n"
  "    email                : amc82@cam.ac.uk\n"
  "\n"
  "**************************************************************************/\n"
  "/***************************************************************************\n"
  " *                                                                         *\n"
  " *   This program is free software; you can redistribute it and/or modify  *\n"
  " *   it under the terms of the GNU General Public License as published by  *\n"
  " *   the Free Software Foundation; either version 2 of the License, or     *\n"
  " *   (at your option) any later version.                                   *\n"
  " *                                                                         *\n"
  " ***************************************************************************/\n"
  ;


string ForceField::ff_banner = 
  "  /***************************************************************************\n"
  "\n"
  "                  IMPLEMENTATION OF FORCE FIELD\n"
  "                    ELEC FUNCTIONAL __ELEC__\n"
  "                    ELEC CUTOFF     __ELEC_CUT__\n"
  "                    VDW  FUNCTIONAL __VDW__\n"
  "                    VDW  CUTOFF     __VDW_CUT__\n"
  "                    SOLVATION       __SOLV__\n"
  "                    GEOMETRY        __GEOM__\n"
  "                    CPU             __CPU__\n"
  "\n"
  "  ***************************************************************************/\n"
  ;


string ForceField::ff_options = 
  "  template<typename float_type>\n"
  "  struct __TYPE__FFOptions {\n"
  "    __NBTYPE__Options<float_type> nbond_options;\n"
  "\n"
  "    __TYPE__FFOptions(){\n"
  "      nbond_options.nblist_options.CUTON = 10.0;\n"
  "      nbond_options.nblist_options.CUTOFF = 12.0;\n"
  "      nbond_options.nblist_options.CUTNB = 14.0;\n"
  "      nbond_options.elec_options.E0 = 332.0716;\n"
  "      nbond_options.elec_options.EPS = 1.0;\n"
  "      nbond_options.elec_options.E14 = 1.0;\n"
  "    }\n"
  "  };\n"
  ;


string ForceField::ff_options_pbc = 
  "  template<typename float_type>\n"
  "  struct __TYPE__FFOptions {\n"
  "    __NBTYPE__Options<float_type> nbond_options;\n"
  "\n"
  "    __TYPE__FFOptions(){\n"
  "      nbond_options.nblist_options.CUTON = __CUTON__;\n"
  "      nbond_options.nblist_options.CUTOFF = __CUTOFF__;\n"
  "      nbond_options.nblist_options.CUTNB = __CUTNB__;\n"
  "      nbond_options.nblist_options.LX = __LX__;\n"
  "      nbond_options.nblist_options.LY = __LX__;\n"
  "      nbond_options.nblist_options.LZ = __LX__;\n"
  "      nbond_options.nblist_options.K = __K__;\n"
  "      nbond_options.elec_options.E0 = __E0__;\n"
  "      nbond_options.elec_options.EPS = __EPS__;\n"
  "      nbond_options.elec_options.E14 = __E14__;\n"
  "    }\n"
  "  };\n"
  ;

string ForceField::options_ostream =
  "  template<typename float_type>\n"
  "  inline ostream & operator<<(ostream & out, \n"
  "			      const __TYPE__FFOptions<float_type> & opt){\n"
  "    out<<\"\t>> Elec:                   \";\n"
  "    out<<\"     __ELEC__\"<<\"\\n\";\n"
  "    out<<\"\t>> Elec-cutoff:            \";\n"
  "    out<<\"     __ELEC_CUT__\"<<\"\\n\";\n"
  "    out<<\"\t>> VdW:                    \";\n"
  "    out<<\"     __VDW__\"<<\"\\n\";\n"
  "    out<<\"\t>> VdW-cutoff:             \";\n"
  "    out<<\"     __VDW_CUT__\"<<\"\\n\";\n"
  "    out<<\"\t>> Solvation:              \";\n"
  "    out<<\"     __SOLV__\"<<\"\\n\";\n"
  "    out<<\"\t>> Cut on:                      \"\n"
  "       <<opt.nbond_options.nblist_options.CUTON<<\"\\n\";\n"
  "    out<<\"\t>> Cut off:                     \"\n"
  "       <<opt.nbond_options.nblist_options.CUTOFF<<\"\\n\";\n"
  "    out<<\"\t>> Cut nb:                      \"\n"
  "       <<opt.nbond_options.nblist_options.CUTNB<<\"\\n\";\n"
  "    out<<\"\t>> E_zero:                      \"\n"
  "       <<opt.nbond_options.elec_options.E0<<\"\\n\";\n"
  "    out<<\"\t>> E_14:                        \"\n"
  "       <<opt.nbond_options.elec_options.E14<<\"\\n\";\n"
  "    out<<\"\t>> Eps:                         \"\n"
  "       <<opt.nbond_options.elec_options.EPS<<\"\\n\";\n"
  "    return out;\n"
  "  }\n"
  "\n"
  ;


string ForceField::options_ostream_pbc =
  "  template<typename float_type>\n"
  "  inline ostream & operator<<(ostream & out, \n"
  "			      const __TYPE__FFOptions<float_type> & opt){\n"
  "    out<<\"\t>> Elec:                   \";\n"
  "    out<<\"     __ELEC__\"<<\"\\n\";\n"
  "    out<<\"\t>> Elec-cutoff:            \";\n"
  "    out<<\"     __ELEC_CUT__\"<<\"\\n\";\n"
  "    out<<\"\t>> VdW:                    \";\n"
  "    out<<\"     __VDW__\"<<\"\\n\";\n"
  "    out<<\"\t>> VdW-cutoff:             \";\n"
  "    out<<\"     __VDW_CUT__\"<<\"\\n\";\n"
  "    out<<\"\t>> Solvation:              \";\n"
  "    out<<\"     __SOLV__\"<<\"\\n\";\n"
  "    out<<\"\t>> Cut on:                      \"\n"
  "       <<opt.nbond_options.nblist_options.CUTON<<\"\\n\";\n"
  "    out<<\"\t>> Cut off:                     \"\n"
  "       <<opt.nbond_options.nblist_options.CUTOFF<<\"\\n\";\n"
  "    out<<\"\t>> Cut nb:                      \"\n"
  "       <<opt.nbond_options.nblist_options.CUTNB<<\"\\n\";\n"
  "    out<<\"\t>> Box size:                    \"\n"
  "       <<opt.nbond_options.nblist_options.LX<<\" x \"\n"
  "       <<opt.nbond_options.nblist_options.LY<<\" x \"\n"
  "       <<opt.nbond_options.nblist_options.LZ<<\"\\n\";\n"
  "    out<<\"\t>> Granularity:                 \"\n"
  "       <<opt.nbond_options.nblist_options.K<<\"\\n\";\n"
  "    out<<\"\t>> E_zero:                      \"\n"
  "       <<opt.nbond_options.elec_options.E0<<\"\\n\";\n"
  "    out<<\"\t>> E_14:                        \"\n"
  "       <<opt.nbond_options.elec_options.E14<<\"\\n\";\n"
  "    out<<\"\t>> Eps:                         \"\n"
  "       <<opt.nbond_options.elec_options.EPS<<\"\\n\";\n"
  "    return out;\n"
  "  }\n"
  "\n"
  ;



string ForceField::ff_data = 
  "  template<typename FloatType>\n"
  "  struct __TYPE__Data : \n"
  "    public __KIND__InternalEnergyData<FloatType>,\n"
  "    public __TYPE__NBondEnergy<FloatType>::NBEnergyDataType \n"
  "  {\n"
  "    typedef FloatType float_type;\n"
  "    float_type energy;\n"
  "    float_type const_energy;\n"
  "\n"
  "    void write(ofstream & rstOut){\n"
  "      rstOut.write((char *)this, sizeof(__TYPE__Data));\n"
  "    }\n"
  "\n"
  "    void read(ifstream & rstIn){\n"
  "      rstIn.read(( char *)this, sizeof(__TYPE__Data));\n"
  "    }\n"
  "  };\n"
  ;


string ForceField::ff =
  "  template<typename float_type>\n"
  "  class __TYPE__ForceField __PARENT__ {\n"
  "    C22InternalEnergy<float_type> internal_energy;\n"
  "    PBCCSSNNBondEnergy<float_type>   nbond_energy;\n"
  "    __TYPE__Data<float_type> data_;\n"
  "    __TYPE__FFOptions<float_type> options_;\n"
  "\n"
  "#ifdef ALM_MPI_FF\n"
  "    Coor<float_type> f;\n"
  "#endif\n"
  "    vector<Constraint *>   constraints_;\n"
  "    vector<NBConstraintBase *> nb_constraints_;\n"
  "  public:\n"
  "    typedef typename PBCCSSNNBondEnergy<float_type>::Options nbond_options;\n"
  "    typedef __TYPE__FFOptions<float_type> Options;\n"
  "    typedef __TYPE__Data<float_type> DataType;\n"
  "    typedef float_type FloatType;\n"
  "\n"
  "    enum { BOND  = 1,\n"
  "	   ANGLE = 1,\n"
  "	   DIHE  = 1,\n"
  "	   IMPH  = 1,\n"
  "	   UB    = __UB__,\n"
  "	   ELEC  = 1,\n"
  "	   VDW   = 1,\n"
  "	   SOLV  = 0,\n"
  "	   SSE   = 0,\n"
  "	   PBC   = __PBC__\n"
  "    };\n"
  "\n"
  "    __TYPE__ForceField(const Molecules & molecules,\n"
  "			 const MDB & mdb,\n"
  "			 const Options & options)\n"
  "      :internal_energy(molecules,mdb),\n"
  "       nbond_energy(molecules,mdb,options.nbond_options),\n"
  "       options_(options)\n"
  "#ifdef ALM_MPI_FF\n"
  "      ,f(molecules.atom_size())\n"
  "#endif\n"
  "    {\n"
  "      if(almost::verbose()){\n"
  "	cout<<\"\\t>> __KIND__ force field created with options:\\n\\n\";\n"
  "	cout<<options_<<\"\\n\";\n"
  "	cout<<\"\\t>> Constraints:\\n\";\n"
  "	cout<<\"\\t   \"<<constraints_.size()<<\" constraints added to base force field\\n\\n\";\n"
  "	cout<<\"\\t>> Build details:\\n\\n\";\n"
  "	cout<<nbond_energy.exclist_data()<<\"\\n\\n\";\n"
  "      }\n"
  "    }\n"
  "\n"
  "    __TYPE__ForceField(const Molecules & molecules,\n"
  "			 const MDB & mdb,\n"
  "			 const Options & options,\n"
  "			 ConstraintCollection &cc)\n"
  "      :internal_energy(molecules,mdb),\n"
  "       nbond_energy(molecules,mdb,options.nbond_options),\n"
  "       options_(options)\n"
  "#ifdef ALM_MPI_FF\n"
  "      ,f(molecules.atom_size())\n"
  "#endif\n"
  "    {\n"
  "      copy(constraints_,cc);\n"
  "      if(almost::verbose()){\n"
  "	cout<<\"\\t>> __KIND__ force field created with options:\\n\\n\";\n"
  "	cout<<options_<<\"\\n\";\n"
  "	cout<<\"\\t>> Constraints:\\n\";\n"
  "	cout<<\"\\t   \"<<constraints_.size()<<\" constraints added to base force field\\n\\n\";\n"
  "	cout<<\"\\t>> Build details:\\n\\n\";\n"
  "	cout<<nbond_energy.exclist_data()<<\"\\n\\n\";\n"
  "      }\n"
  "    }\n"
  "\n"
  "    ~__TYPE__ForceField(){\n"
  "      for(int i = 0;i<constraints_.size();i++){\n"
  "	delete constraints_[i];\n"
  "      }\n"
  "    }\n"
  "\n"
  "    template<typename CoorType>\n"
  "    float_type energy(const CoorType & coor, bool update = true)"
  "__ENERGY_STRING__"
  "\n"
  "    template<typename CoorType>\n"
  "    float_type intenergy(const CoorType & coor, bool update = true)"
  "__INTENERGY_STRING__"
  "\n"
  "    void force(const Coor<float_type> & coor, Coor<float_type> & ff, \n"
  "               bool update)"
  "__FORCE_STRING__"
  "\n"
  "    void intforce(const Coor<float_type> & coor, Coor<float_type> & ff)"
  "__INTFORCE_STRING__"
  "\n"
  "    void energy_force(const Coor<float_type> & coor, Coor<float_type> & ff, \n"
  "		      bool update)"
  "__ENERGY_FORCE_STRING__"
  "\n"
  "    const DataType & data(){ return data_;}\n"
  "    void copy(vector<Constraint *> &constraints, \n"
  "	      ConstraintCollection &cc){\n"
  "      int nconst = cc.constraints().size();\n"
  "      for(int i=0;i<nconst;i++){\n"
  "	Constraint * c = cc.constraints()[i]->clone();\n"
  "	constraints.push_back(c);\n"
  "      }\n"
  "    }\n"
  "    \n"
  "  };\n"
  ;

string ForceField::energy_string = 
  "{\n"
  "      data_.internal_energy = internal_energy.energy(coor);\n"
  "      data_.nbond_energy = nbond_energy.energy(coor,update);\n"
  "      data_.energy = data_.internal_energy + data_.nbond_energy;\n"
  "\n"
  "      data_.const_energy = 0;\n"
  "#ifdef ALM_MPI_FF\n"
  "      data_.bond_energy = internal_energy.data().bond_energy;\n"
  "      data_.angle_energy = internal_energy.data().angle_energy;\n"
  "      data_.dihe_energy = internal_energy.data().dihe_energy;\n"
  "      data_.imph_energy = internal_energy.data().imph_energy;\n"
  "      data_.ub_energy = internal_energy.data().ub_energy;\n"
  "      data_.elec_energy = nbond_energy.nb_data().elec_energy;\n"
  "      data_.vdw_energy = nbond_energy.nb_data().vdw_energy;\n"
  "      almost_reduce_energy(data_);\n"
  "#endif\n"
  "      for(int i=0;i<constraints_.size();++i)\n"
  "	data_.const_energy += \n"
  "	  constraints_[i]->energy(coor);\n"
  "#ifndef ALM_MPI_FF\n"
  "      data_.energy  += data_.const_energy;\n"
  "      data_.bond_energy = internal_energy.data().bond_energy;\n"
  "      data_.angle_energy = internal_energy.data().angle_energy;\n"
  "      data_.dihe_energy = internal_energy.data().dihe_energy;\n"
  "      data_.imph_energy = internal_energy.data().imph_energy;\n"
  "\n"
  "      data_.elec_energy = nbond_energy.nb_data().elec_energy;\n"
  "      data_.vdw_energy = nbond_energy.nb_data().vdw_energy;\n"
  "#endif      \n"
  "      return data_.energy;\n"
  "    }\n"
  ;


string ForceField::intenergy_string = 
"{\n"
  "      data_.internal_energy = internal_energy.energy(coor);\n"
  "      data_.nbond_energy = 0;\n"
  "      data_.energy = data_.internal_energy + data_.nbond_energy;\n"
  "\n"
  "      data_.const_energy = 0;\n"
  "#ifndef ALM_MPI_FF\n"
  "      for(int i=0;i<constraints_.size();++i)\n"
  "	data_.const_energy += \n"
  "	  constraints_[i]->energy(coor);\n"
  "#endif\n"
  "\n"
  "      data_.energy  += data_.const_energy;\n"
  "      data_.bond_energy = internal_energy.data().bond_energy;\n"
  "      data_.angle_energy = internal_energy.data().angle_energy;\n"
  "      data_.dihe_energy = internal_energy.data().dihe_energy;\n"
  "      data_.imph_energy = internal_energy.data().imph_energy;\n"
  "\n"
  "      data_.elec_energy = 0;\n"
  "      data_.vdw_energy = 0;\n"
  "      \n"
  "#ifdef ALM_MPI_FF\n"
  "      almost_reduce_energy(data_);\n"
  "      return data_.energy;\n"
  "#endif\n"
  "    }\n"
  ;


string ForceField::force_string = 
  "{\n"
  "#ifdef ALM_MPI_FF\n"
  "      f.clear();\n"
  "      nbond_energy.force(coor,f,update);\n"
  "      internal_energy.force(coor,f);\n"
  "      almost_reduce_forces(f,ff);\n"
  "#else\n"
  "      nbond_energy.force(coor,ff,update);\n"
  "      internal_energy.force(coor,ff);\n"
  "#endif\n"
  "    }\n"
  ;

string ForceField::intforce_string = 
  "{\n"
  "#ifdef ALM_MPI_FF\n"
  "      f.clear();\n"
  "      internal_energy.force(coor,f);\n"
  "      almost_reduce_forces(f,ff);\n"
  "#else\n"
  "      ff.clear();\n"
  "      internal_energy.force(coor,ff);\n"
  "#endif\n"
  "    }\n"
  ;

string ForceField::energy_force_string = 
  "{\n"
  "#ifdef ALM_MPI_FF\n"
  "      f.clear();\n"
  "      nbond_energy.energy_force(coor,f,update);\n"
  "      internal_energy.energy_force(coor,f);\n"
  "\n"
  "      data_.energy = internal_energy.data().internal_energy \n"
  "	+ nbond_energy.nb_data().vdw_energy \n"
  "	+ nbond_energy.nb_data().elec_energy\n"
  "	;\n"
  "\n"
  "\n"
  "      data_.bond_energy = internal_energy.data().bond_energy;\n"
  "      data_.angle_energy = internal_energy.data().angle_energy;\n"
  "      data_.dihe_energy = internal_energy.data().dihe_energy;\n"
  "      data_.imph_energy = internal_energy.data().imph_energy;\n"
  "      data_.ub_energy = internal_energy.data().ub_energy;\n"
  "      data_.elec_energy = nbond_energy.nb_data().elec_energy;\n"
  "      data_.vdw_energy = nbond_energy.nb_data().vdw_energy;\n"
  "      almost_reduce_energy(data_);\n"
  "      almost_reduce_forces(f,ff);\n"
  "#else\n"
  "      nbond_energy.energy_force(coor,ff,update);\n"
  "      internal_energy.energy_force(coor,ff);\n"
  "\n"
  "      data_.energy = internal_energy.data().internal_energy \n"
  "	+ nbond_energy.nb_data().nbond_energy;\n"
  "\n"
  "      data_.bond_energy = internal_energy.data().bond_energy;\n"
  "      data_.angle_energy = internal_energy.data().angle_energy;\n"
  "      data_.dihe_energy = internal_energy.data().dihe_energy;\n"
  "      data_.imph_energy = internal_energy.data().imph_energy;\n"
  "      data_.ub_energy = internal_energy.data().ub_energy;\n"
  "      data_.elec_energy = nbond_energy.nb_data().elec_energy;\n"
  "      data_.vdw_energy = nbond_energy.nb_data().vdw_energy;\n"
  "#endif      \n"
  "    }\n"
  ;


string ForceField::intenergy_force_string = "__NOT_INMPLEMENTED__";
void ForceField::mk_name(){
  string tmp = "";//"PBC";
  if(elec==CDIE) tmp += "C";
  if(elec==RDIE) tmp += "R";
  if(elec==NONE) tmp += "N";

  if(elec_cutoff==SHIFT) tmp+="S";
  if(elec_cutoff==SWITCH) tmp+="Sw";

  if(vdw_cutoff==SHIFT) tmp+="S";
  if(vdw_cutoff==SWITCH) tmp+="Sw";

  if(solv==SASA) tmp+="SASA";
  if(solv==EEF1) tmp+="EEF1";
  if(solv==NO) tmp+="N";

  if(geometry==PBC)
    tmp = "PBC"+tmp;

    
  nb_name = tmp;

  if(kind==C19)
    tmp += "C19";
  else if(kind==C22)
    tmp += "C22";
  else if(kind==AMBER)
    tmp += AMBER;

  name = tmp;
}

void ForceField::mk_banner(){
  char elec_string[3][5] = {"CDIE","RDIE","NONE"};
  char vdw_string[1][5] = { "LJ  "};
  char cutoff_string[2][7] = {"SHIFT ","SWITCH"};
  char solv_string[3][5] = { "SASA","EEF1","NO"};
  char cpu_string[2][6] = {"PLAIN","VECT "};
  char geom_string[2][6] = {"PLAIN","PBC"};
  string tmp = ff_banner;
 
  tmp = replace_all_copy(tmp,"__ELEC__",elec_string[elec]);
  tmp = replace_all_copy(tmp,"__ELEC_CUT__",cutoff_string[elec_cutoff]);
  tmp = replace_all_copy(tmp,"__VDW__",vdw_string[vdw]);
  tmp = replace_all_copy(tmp,"__VDW_CUT__",cutoff_string[vdw_cutoff]);
  tmp = replace_all_copy(tmp,"__SOLV__",solv_string[solv]);
  tmp = replace_all_copy(tmp,"__GEOM__",geom_string[geometry]);
  tmp = replace_all_copy(tmp,"__CPU__",cpu_string[cpu]);

  tmp = replace_all_copy(tmp,"__TYPE__",name);
  cout<<tmp;
}



int main(int argc, char ** argv){
  string banner = file_banner;
  cout<<replace_all_copy(banner,"__FILE__","pbcc22none.h");

  cout<<"#ifndef __PBCC22NONE__\n";
  cout<<"#define __PBCC22NONE__\n";
  cout<<"\n";
  cout<<"#include <almost.h>\n";
  cout<<"#ifdef ALM_MPI_FF\n";
  cout<<"#include <forcefield/ffutils.h>\n";
  cout<<"#endif\n";
  cout<<"namespace Almost {\n";

  for(int elec = 0; elec<2;elec++){
    for(int elec_cut = 0; elec_cut<2;elec_cut++){
      for(int vdw = 0; vdw<1;vdw++){
	for(int vdw_cut = 0; vdw_cut<2;vdw_cut++){
	  
	  ForceField ff;
	  //  ff.name = "PBCCSSNC22";
	  //	  ff.parent = "AlmostForceField";
	  ff.full = true;
	  ff.elec = elec;// ForceField::CDIE;
	  ff.elec_cutoff = elec_cut;//ForceField::SHIFT;
	  ff.vdw = vdw;//ForceField::LJ;
	  ff.vdw_cutoff = vdw_cut;//ForceField::SHIFT;
	  ff.solv = ForceField::NO;
	  ff.cpu = ForceField::PLAIN;
	  ff.geometry = ForceField::PBC;
	  ff.kind = ForceField::C22;
	  ff.mk_name();
	  ff.mk_banner();
	  ff.mk_options();
	  ff.mk_options_ostream();
	  ff.mk_data();
	  ff.mk_ff();
	}
      }
    }
  }


  cout<<"}\n";
  cout<<"#endif\n";
}