/***************************************************************************
                          alm_mknb.cpp  -  description
                             -------------------
    begin                : Mon Aug 22 13:19:08 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 <stdio.h>
#include <unistd.h>
#include <string>
#include <vector>
#include <iostream>
#include <iterator>
#include <fstream>
#include <map>
#include <set>
#include <sstream>
using namespace std;

static FILE * mknb_output;// = stdout;
//static int mknb_double =1;

enum { ESHIFT, ESWITCH, ETRUNC,
       VSHIFT, VSWITCH, VTRUNC, 
       CDIE, RDIE, ENONE,
       LJ, /*....*/
       ENE,FORCE, ENE_FORCE,
       ATOM_ATOM,WAT_ATOM,
       WAT_WAT,
       PBC,PLAIN,
       DOUBLE,FLOAT
};


void mknb_file_header(const char * file){
  time_t t = time(NULL);
  fprintf(mknb_output,"/***************************************************************************\n");
  fprintf(mknb_output,"                          %s  -  description\n",file);
  fprintf(mknb_output,"                             -------------------\n");
  fprintf(mknb_output,"    begin                : %s",ctime(&t));
  fprintf(mknb_output,"    copyright            : (C) 2005 by Cavalli Andrea\n");
  fprintf(mknb_output,"    email                : cavalli@bioc.unizh.ch\n");
  fprintf(mknb_output,"**************************************************************************/\n");
  fprintf(mknb_output,"/***************************************************************************\n");
  fprintf(mknb_output," *                                                                         *\n");
  fprintf(mknb_output," *   This program is free software; you can redistribute it and/or modify  *\n");
  fprintf(mknb_output," *   it under the terms of the GNU General Public License as published by  *\n");
  fprintf(mknb_output," *   the Free Software Foundation; either version 2 of the License, or     *\n");
  fprintf(mknb_output," *   (at your option) any later version.                                   *\n");
  fprintf(mknb_output," *                                                                         *\n");
  fprintf(mknb_output," ***************************************************************************/\n");



}


struct nb_function_param {
  string name;
  string type;
  bool pointer;
  bool constant;
  nb_function_param(string n,string t,bool c, bool p):
    name(n),type(t),pointer(p),constant(c){}

  string to_string(){
    stringstream ss;
    if(constant)
      ss<<"const ";
    ss<<type<<" ";
    if(pointer) ss <<"*";
    ss<<name;
    return ss.str();
  }
};

struct mknb_function {

  int elec_type;
  int cdie;
  int vdw_type;
  int geom;

  int kind;
  int type;
  //
  int bCompat;
  int ucc; //make all function look the same
  int ni,nj;
  int safe; //deos no check cutoff
  int float_type;
  bool calc_vir;
  mknb_function(){
    ni = nj =1;
  }

  vector<nb_function_param> ene_par;
  void mkenergy_param(){
    {
    nb_function_param par("int","first",true,true);
    ene_par.push_back(par);
    }
//     const int *second,
//       const int *boundary,
//       const int natoms,
//       const double *coor,
//       const double c2,
//       const double c_on2,
//       const double *charge,
//       const int *type_code,
//       const int type_size,
//       const double *vdw_data,
//       double *ene_vdw, double *ene_cul
  }

};




map<string,set<string> > int_vars;

int init_int_vars(){
  ifstream in;
  in.open("mknb_intvars.data");
  istream_iterator<string> iter(in),end;
  
  int count = 0;
  while(iter!=end){
    string f = *iter; iter++;
    string s = *iter; iter++;
    int_vars[f].insert(s);
    ++count;
  }
  cout<<"Read "<<count<<" variables"<<endl;
  return count;
}

static int init_int_vars__ = init_int_vars();

vector<string> int_variables;


mknb_function nb_func;


string spec2name(){
  string n;
  if(nb_func.kind==ENE)
    n = "energy_";
  else if(nb_func.kind==FORCE)
    n = "force_";
  else n = "energy_force_";


  if(nb_func.geom==PBC) n = "pbc"+n;
  if((nb_func.geom==PBC)&&nb_func.calc_vir) n = "v"+n;
  if(nb_func.elec_type==ENONE)
    n += "xx";

  if(nb_func.cdie) n +="c";
  
  if((!nb_func.cdie)&&(nb_func.elec_type!=ENONE)) n +="r";

  if(nb_func.elec_type==ETRUNC) n += "t";
  if(nb_func.elec_type==ESHIFT) n += "s";
  if(nb_func.elec_type==ESWITCH) n += "sw";

  if(nb_func.vdw_type==VTRUNC) n += "t";
  if(nb_func.vdw_type==VSHIFT) n += "s";
  if(nb_func.vdw_type==VSWITCH) n += "sw";


  n += "n";
  if(nb_func.float_type==FLOAT) n+="_f";
  else n+="_d";
  
  return n;
}

void mknb_func_desc(){
  fprintf(mknb_output,"\n\n/***************************************************************************\n\n");
  fprintf(mknb_output,"                     IMPLEMENTATION OF NB ENERGY\n");
  string attr;
  if(nb_func.cdie) attr = "CDIE";
  else if(nb_func.elec_type==ENONE) attr = "NONE";
  else attr = "RDIE";
  fprintf(mknb_output,"                       Elec Functional %s\n",attr.c_str());  
  if(nb_func.elec_type==ETRUNC) attr = "TRUNC";
  if(nb_func.elec_type==ESHIFT) attr = "SHIFT";
  if(nb_func.elec_type==ESWITCH) attr = "SWITCH";
  if(nb_func.elec_type==ENONE) attr = "NONE";
  fprintf(mknb_output,"                       Elec CUTOFF     %s\n",attr.c_str());
  fprintf(mknb_output,"                       VDW  Functional LJ\n");
  if(nb_func.vdw_type==VTRUNC) attr = "TRUNC";
  if(nb_func.vdw_type==VSHIFT) attr = "SHIFT";
  if(nb_func.vdw_type==VSWITCH) attr = "SWITCH";
  fprintf(mknb_output,"                       VDW  CUTOFF     %s\n",attr.c_str());
  fprintf(mknb_output,"                       SOLVATION       NONE\n");
  fprintf(mknb_output,"                       Sol  Functional SolvNone\n\n");
  if(nb_func.geom==PBC) attr = "PBC";
  else attr = "PLAIN";
  fprintf(mknb_output,"                       Geometry        %s\n\n",attr.c_str());
  if(nb_func.calc_vir)
    attr = "yes";
  else attr = "no";
  fprintf(mknb_output,"                       Compute virial  %s\n\n",attr.c_str());
  if(nb_func.float_type==DOUBLE) attr = "double";
  else attr = "float";
  fprintf(mknb_output,"                       Float type      %s\n",attr.c_str());
  fprintf(mknb_output,"                       CPU             CPUPlain\n\n");
  fprintf(mknb_output," ***************************************************************************/\n");
}



int cut_off(mknb_function &nb_func){
  if(nb_func.vdw_type==VTRUNC) return false;
  if((nb_func.elec_type!=ENONE)&&(nb_func.elec_type!=ETRUNC)) return true;
  return true;
}

#define eshift (nb_func.elec_type==ESHIFT)
#define eswitch (nb_func.elec_type==ESWITCH)
#define etrunc (nb_func.elec_type==ETRUNC)
#define rdie   (!nb_func.cdie)

#define vshift (nb_func.vdw_type==VSHIFT)
#define vswitch (nb_func.vdw_type==VSWITCH)
#define vtrunc (nb_func.vdw_type==VTRUNC)

#define CUTOFF cut_off(nb_func)

int force(mknb_function &nb_func){
  if((nb_func.kind == FORCE)||(nb_func.kind == ENE_FORCE))
    return true;
  return false;

}

#define CFORCE force(nb_func)


int energy(mknb_function &nb_func){
  return !(nb_func.kind==FORCE);
}

#define CENE energy(nb_func)

char buff[1024];

void assign(char * a, char * b){
  fprintf(mknb_output,"%s = %s;\n",a,b);
}

void assign(char * a, char * b, char *op, char *c){
  fprintf(mknb_output,"%s = %s %s %s;\n",a,b,op,c);
}

char* array(char * a, char * b){
  sprintf(buff,"%s[%s]",a,b);
  return buff;
}

char * float_type(){
  if(nb_func.float_type==DOUBLE)
    sprintf(buff,"double");
  else sprintf(buff,"float");
  return buff;
}

void mknb_force_function_header(const char *funcname){

  fprintf(mknb_output,"void %s(\n",funcname); 
  fprintf(mknb_output,"\t%s %s,\n","const int *", "first");
  fprintf(mknb_output,"\t%s %s,\n","const int *", "second");
  fprintf(mknb_output,"\t%s %s,\n","const int *", "boundary");
  if(nb_func.calc_vir)
    fprintf(mknb_output,"\t%s %s,\n","const int *", "box");
  fprintf(mknb_output,"\t%s %s,\n","const int", "natoms");
  fprintf(mknb_output,"\t%s %s %s %s,\n","const",float_type(), "*", "coor");
  if(nb_func.geom==PBC)
    fprintf(mknb_output,"\t%s %s %s %s,\n","const",float_type(),
	    "*", "sh_coor");
  if(CFORCE)
    fprintf(mknb_output,"\t%s %s %s,\n",float_type(), "*", "force");

  /* This is cutoff area .... */
  fprintf(mknb_output,"\t%s %s %s,\n","const", float_type(), "c2");

  if((nb_func.ucc==1)||((nb_func.elec_type==ESWITCH)||(nb_func.vdw_type==VSWITCH))){
    fprintf(mknb_output,"\t%s %s %s,\n","const",float_type() , "c_on2");
  }

  if(nb_func.elec_type!=ENONE)
    fprintf(mknb_output,"\t%s %s %s %s,\n","const", float_type(), 
	    "*", "charge");

  fprintf(mknb_output,"\t%s %s,\n","const int *", "type_code");
  fprintf(mknb_output,"\t%s %s,\n","const int", "type_size");

  if(nb_func.calc_vir)
    fprintf(mknb_output,"\t%s %s %s,\n",float_type(), "*", "fvir"); 


  if((nb_func.kind==ENE)||(nb_func.kind==ENE_FORCE)){
    fprintf(mknb_output,"\t%s %s %s %s,\n","const", float_type(), "*", "vdw_data");
    fprintf(mknb_output,"\t%s %s %s,\n",float_type(), "*", "ene_vdw");
    fprintf(mknb_output,"\t%s %s %s)\n",float_type(), "*", "ene_cul");
  }
  else 
    fprintf(mknb_output,"\t%s %s %s %s)\n","const", float_type(), "*", "vdw_data");
}


char * mknb_int(const char * name){
  sprintf(buff,"int %s;\n",name);
  return buff;
}


char * mknb_float(const char * name){
  if(nb_func.float_type==DOUBLE)
    sprintf(buff,"double %s;\n",name);
  else
    sprintf(buff,"float %s;\n",name);
  return buff;
}

char * mknb_float(const char * name, char * init){
  if(nb_func.float_type==DOUBLE)
    sprintf(buff,"double %s = %s;\n",name,init);
  else
    sprintf(buff,"float %s = %s;\n",name,init);
  return buff;
}

char * mknb_const_int(const char * name, char * init){
  sprintf(buff,"const int %s = %s;\n",name,init);
  return buff;
}

char * mknb_const_float(const char * name,const char * init){
  if(nb_func.float_type==DOUBLE)
    sprintf(buff,"const double %s = %s;\n",name,init);
  else
    sprintf(buff,"const float %s = %s;\n",name,init);
  return buff;
}

#define needed(s,v) (int_vars[s].find(v)!=int_vars[n].end())

void 
mknb_force_function_internal_vars(){
  //get name
  string n;
  if(nb_func.kind==ENE) n="energy_";
  else if(nb_func.kind==FORCE) n="force_";
  else n="energy_force_";

  if(nb_func.geom==PBC) n = "pbc"+ n;
  if((nb_func.geom==PBC)&&nb_func.calc_vir) n = "v"+n;
  if(nb_func.cdie) n +="c";
  else {
    if(nb_func.elec_type==ENONE) n+= "xx";
    else n +="r";
  }

  if(nb_func.elec_type==ETRUNC) n+="t";
  if(nb_func.elec_type==ESHIFT) n+="s";
  if(nb_func.elec_type==ESWITCH) n+="sw";

  if(nb_func.vdw_type==VTRUNC) n+="t";
  if(nb_func.vdw_type==VSHIFT) n+="s";
  if(nb_func.vdw_type==VSWITCH) n+="sw";

  n +="n";

  if(nb_func.float_type==DOUBLE) n += "_d";
  else n += "_f";
  int_variables.clear();

  int_variables.push_back(mknb_int("i"));
  int_variables.push_back(mknb_int("j"));
  int_variables.push_back(mknb_int("f"));
  int_variables.push_back(mknb_int("fpos"));
  int_variables.push_back(mknb_int("s"));
  int_variables.push_back(mknb_int("spos"));
  if(needed(n,"bpos"))
    int_variables.push_back(mknb_int("bpos"));
  int_variables.push_back(mknb_int("b1"));
  int_variables.push_back(mknb_int("b2"));
  int_variables.push_back(mknb_int("tc1"));
  int_variables.push_back(mknb_int("tc2"));
  int_variables.push_back(mknb_int("data_pos"));
  
  if(needed(n,"fracX"))
    int_variables.push_back(mknb_float("fracX"));
  if(needed(n,"fracY"))
    int_variables.push_back(mknb_float("fracY"));
  if(needed(n,"fracZ"))
    int_variables.push_back(mknb_float("fracZ"));


  int_variables.push_back(mknb_float("dx,dy,dz"));
  int_variables.push_back(mknb_float("ix,iy,iz"));
  int_variables.push_back(mknb_float("jx,jy,jz"));

  if(needed(n,"tx"))
    int_variables.push_back(mknb_float("tx,ty,tz"));

  if(needed(n,"shX"))
    int_variables.push_back(mknb_float("shX,shY,shZ"));
  

  int_variables.push_back(mknb_float("d2"));

  if(needed(n,"d6"))
    int_variables.push_back(mknb_float("d6"));

  if(needed(n,"dinv"))
    int_variables.push_back(mknb_float("dinv"));


  int_variables.push_back(mknb_float("d2inv"));
  int_variables.push_back(mknb_float("d6inv"));
  int_variables.push_back(mknb_float("d12inv"));

  if(needed(n,"vsh"))
    int_variables.push_back(mknb_float("vsh"));

  if(needed(n,"esh"))
    int_variables.push_back(mknb_float("esh"));

  if(needed(n,"esh2"))
    int_variables.push_back(mknb_float("esh2"));

  if(needed(n,"sh1"))
    int_variables.push_back(mknb_float("sh1"));
  if(needed(n,"sh2"))
    int_variables.push_back(mknb_float("sh2"));
  if(needed(n,"sw"))
    int_variables.push_back(mknb_float("sw"));
  if(needed(n,"sw1"))
    int_variables.push_back(mknb_float("sw1"));


  if(needed(n,"frac"))
    int_variables.push_back(mknb_float("frac"));
  if(needed(n,"fvdw"))
    int_variables.push_back(mknb_float("fvdw"));
  if(needed(n,"felec"))
    int_variables.push_back(mknb_float("felec"));
  if(needed(n,"eelec"))
    int_variables.push_back(mknb_float("eelec"));
  

  if(needed(n,"nb6")){
    int_variables.push_back(mknb_float("nb6"));
    int_variables.push_back(mknb_float("nb12"));
  }

  if(needed(n,"enb"))
    int_variables.push_back(mknb_float("enb"));

  int_variables.push_back(mknb_float("A,B"));

  if(needed(n,"nbelec"))
    int_variables.push_back(mknb_float("nbelec"));
  
  if(needed(n,"q1"))
    int_variables.push_back(mknb_float("q1,q2"));
  

  if(needed(n,"evdw"))
    int_variables.push_back(mknb_float("evdw","0"));
  if(needed(n,"ecul"))
    int_variables.push_back(mknb_float("ecul","0"));
  

  if(needed(n,"c2inv"))
    int_variables.push_back(mknb_const_float("c2inv","1.0/c2"));
  if(needed(n,"c6inv"))
    int_variables.push_back(mknb_const_float("c6inv","c2inv*c2inv*c2inv"));
  if(needed(n,"c12inv"))
    int_variables.push_back(mknb_const_float("c12inv","c6inv*c6inv"));
  if(needed(n,"c18inv"))
    int_variables.push_back(mknb_const_float("c18inv","c12inv*c6inv"));
  
  
  if(needed(n,"c2diffinv"))
    int_variables.push_back(mknb_const_float("c2diffinv","1.0/((c2-c_on2)*(c2-c_on2)*(c2-c_on2))"));

  if(needed(n,"one"))
    int_variables.push_back(mknb_const_float("one","1.0"));

  if(needed(n,"two"))
    int_variables.push_back(mknb_const_float("two","2.0"));

  if(needed(n,"three"))
    int_variables.push_back(mknb_const_float("three","3.0"));

  if(needed(n,"four"))
    int_variables.push_back(mknb_const_float("four","4.0"));

  if(needed(n,"six")){
    int_variables.push_back(mknb_const_float("six","6.0"));
    int_variables.push_back(mknb_const_float("twelve","12.0"));
  }

  if(needed(n,"six_inv")){
    int_variables.push_back(mknb_const_float("six_inv","1.0/6.0"));
    int_variables.push_back(mknb_const_float("twelve_inv","1.0/12.0"));
  }


  for(int i=0;i<int_variables.size();i++){
    fprintf(mknb_output,"%s",int_variables[i].c_str());
  }
}

void mknb_inner_boundary(){
  if(!nb_func.bCompat){
    assign("b1",array("boundary","2*i"));
    assign("b2",array("boundary","2*i+1"));
  }
  else {
    assign("b1",array("boundary","i"));
    assign("b2",array("boundary","i+1"));
  }
}

void 
mknb_comment(char *s)
{
  fprintf(mknb_output,"\n\t/*  %s  */\n",s);
}

void
mknb_start_loop(char *lvar,char *from,char *to){
  fprintf(mknb_output,"for(%s=%s; (%s<%s); %s++)", lvar,from,lvar,to,lvar);
}



void mknb_zero_forces(){
  assign("fracX","0");
  assign("fracY","0");
  assign("fracZ","0");
}


void mknb_load_shifts(){
  assign("shX",array("sh_coor","COOR_DIM*i"));
  assign("shY",array("sh_coor","COOR_DIM*i+1"));
  assign("shZ",array("sh_coor","COOR_DIM*i+2"));
}

void mknb_load_first(){
  assign("f",array("first","i"));
  assign("fpos","COOR_DIM*f");
  if(nb_func.calc_vir)
    assign("bpos",array("3*box","i"));
  if(nb_func.geom!=PBC){
    assign("ix",array("coor","fpos"));
    assign("iy",array("coor","fpos+1"));
    assign("iz",array("coor","fpos+2"));
  }
  else {
    mknb_comment("Load Shift vectors");
    mknb_load_shifts();
    assign("ix",array("coor","fpos"),"+","shX");
    assign("iy",array("coor","fpos+1"),"+","shY");
    assign("iz",array("coor","fpos+2"),"+","shZ");
  }
}

void mknb_load_second(){
  assign("s",array("second","j"));
  assign("spos","COOR_DIM*s");
  assign("jx",array("coor","spos"));
  assign("jy",array("coor","spos+1"));
  assign("jz",array("coor","spos+2"));
}


void mknb_load_par_first(){
  if(nb_func.elec_type!=ENONE)
    assign("q1",array("charge","f"));
  assign("tc1",array("type_code","f"));
}


void mknb_load_par_second(){
  if(nb_func.elec_type!=ENONE)
    assign("q2",array("charge","s"));
  assign("tc2",array("type_code","s"));
  assign("data_pos","4*(type_size*tc1+tc2)");
  assign("A",array("vdw_data","data_pos"));
  assign("B",array("vdw_data","data_pos+1"));
}


void mknb_compute_dist(){
  assign("dx","ix-jx");
  assign("dy","iy-jy");
  assign("dz","iz-jz");
  assign("d2","dx*dx+dy*dy+dz*dz");
}

void mknb_d_powers(){
  if((nb_func.elec_type!=ENONE)&&(nb_func.cdie)){
    if(nb_func.float_type==DOUBLE)
      assign("dinv","invsqrt(d2)");
    else 
      assign("dinv","invsqrtf(d2)");
    
    assign("d2inv","dinv*dinv");
  }
  else
    assign("d2inv","one/d2");
  if(nb_func.vdw_type==VSHIFT){
    assign("d6","d2*d2*d2");
  }
  assign("d6inv","d2inv*d2inv*d2inv");
  assign("d12inv","d6inv*d6inv");    
}


void mknb_ene(){
  if(nb_func.elec_type==ESHIFT){
    assign("esh", "one - d2*c2inv");
    assign("esh2", "esh*esh");
  }

  //elec
  if(nb_func.cdie){
    if(nb_func.elec_type==ETRUNC)
      assign("ecul", "ecul + q1*q2*dinv");
    else if(nb_func.elec_type==ESHIFT){
      assign("ecul", "ecul + q1*q2*esh2*dinv");
    }
  }
  else {
    if(nb_func.elec_type==ETRUNC)
      assign("ecul", "ecul + q1*q2*d2inv");
    else if(nb_func.elec_type==ESHIFT){
      assign("ecul", "ecul + q1*q2*esh2*d2inv");
    }
  } 

  //vdw 
  if(nb_func.vdw_type==VTRUNC)
    assign("evdw", "evdw + A*d12inv - B*d6inv");
  else if(nb_func.vdw_type==VSHIFT){
    assign("vsh", "one - d6*c6inv");
    assign("nb12", "A*(d12inv - c12inv*(one + vsh+vsh))");
    assign("nb6", "B*(d6inv -c6inv*(one + vsh))");
    assign("evdw", "evdw + nb12 - nb6");
  }
  //switch ??

  if((nb_func.elec_type==ESWITCH)||(nb_func.vdw_type==VSWITCH)){
    fprintf(mknb_output,"if(d2<c_on2){\n");
    //elec
    if(nb_func.cdie&&(nb_func.elec_type==ESWITCH))
      assign("ecul", "ecul + q1*q2*dinv");
    if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH))
      assign("ecul", "ecul + q1*q2*d2inv");
    
    //vdw
    if(nb_func.vdw_type==VSWITCH){
      assign("evdw ", "evdw + (A*d12inv - B*d6inv)");
    }
    fprintf(mknb_output,"}\n");
    //cumpute swfunction
    fprintf(mknb_output,"else {\n");
    assign("sw", "c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2)");
    
    //elec
    if(nb_func.cdie&&(nb_func.elec_type==ESWITCH))
      assign("ecul", "ecul + q1*q2*dinv*sw");
    if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH))
      assign("ecul", "ecul + q1*q2*d2inv*sw");
    if(nb_func.vdw_type==VSWITCH)
      assign("evdw ", "evdw + (A*d12inv - B*d6inv)*sw");
    fprintf(mknb_output,"}\n");
  }
}


void mknb_force(){
  if(nb_func.elec_type==ESHIFT){
    assign("esh", "one - d2*c2inv");
  }

  //elec
  if(nb_func.cdie){
    if(nb_func.elec_type==ETRUNC)
      assign("felec", "q1*q2*dinv");
    else if(nb_func.elec_type==ESHIFT)
      assign("felec", "q1*q2*esh*dinv*(esh*d2inv+four*c2inv)");
  }
  else {
    if(nb_func.elec_type==ETRUNC)
      assign("felec", "two*q1*q2*d2inv");
    else if(nb_func.elec_type==ESHIFT)
      assign("felec", "two*q1*q2*esh*d2inv*(esh*d2inv+two*c2inv)");
  }

  //vdw
  if(nb_func.vdw_type==VTRUNC){
    assign("fvdw", "A*d12inv-B*d6inv");
  }
  else if(nb_func.vdw_type==VSHIFT){
    assign("nb6", "B*(d6inv-d6*c12inv);");
    assign("nb12", "A*(d12inv-d6*c18inv)");
    assign("fvdw", "(nb12 - nb6)*d2inv");
  }

  {
    //SWITCH ....
    if((nb_func.elec_type==ESWITCH)||(nb_func.vdw_type==VSWITCH)){
      //elec
      if(nb_func.cdie&&(nb_func.elec_type==ESWITCH)){
	assign("felec", "q1*q2*dinv");
      }
      if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH)){
	assign("felec", "two*q1*q2*d2inv");
      }
      if(nb_func.vdw_type==VSWITCH)
	assign("fvdw", "A*d12inv - B*d6inv");
      //cumulate
      if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSWITCH)){
	assign("frac", "(felec+fvdw)*d2inv");
      }
      else if((nb_func.elec_type==ESWITCH)&&(!(nb_func.vdw_type==VSWITCH))){
	assign("frac", "felec*d2inv");
      }
      else if((!(nb_func.elec_type==ESWITCH))&&((nb_func.vdw_type==VSWITCH))){
	assign("frac", "fvdw*d2inv");
      }
      fprintf(mknb_output,"if(d2>c_on2){\n");
      //compute sw sw1
      assign("sh1", "c2-d2");
      assign("sh2", "(c2+two*d2-three*c_on2)");
      assign("sw", "sh1*sh2");
      assign("sw1", "four*(sh1-sh2)");

      if(!(nb_func.vdw_type==VSWITCH)){
	if(nb_func.cdie&&(nb_func.elec_type==ESWITCH)){
	  assign("frac", "frac*sw - felec*sw1;");
	}
	if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH)){
	  assign("frac", "frac*sw - 0.5*felec*sw1;");
	}
      }
      if(!(nb_func.elec_type==ESWITCH)){
	if(nb_func.vdw_type==VSWITCH)
	  assign("frac", "frac*sw - (twelve_inv*A*d12inv-six_inv*B*d6inv)*sw1");
      }
      //switch switch
      if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSWITCH)){
	if(nb_func.cdie)
	  assign("frac", "frac*sw - (twelve_inv*A*d12inv-six_inv*B*d6inv+felec)*sw1");
	else
	  assign("frac", "frac*sw - (twelve_inv*A*d12inv-six_inv*B*d6inv+0.5*felec)*sw1");
      }
      assign("frac", "frac*c2diffinv*sh1");
      fprintf(mknb_output,"}\n");
    }
  }
  //frac contains scalar forces of switch components so need to add shifted  
  if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VSHIFT)){
    assign("frac", "(felec+fvdw)");
  }
  else if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSHIFT))
    assign("frac", "frac + fvdw");
  else if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VSWITCH))
    assign("frac", "frac + felec");
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "(felec+fvdw)*d2inv");
  }
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VSWITCH)){
    assign("frac", "frac + felec*d2inv");
  }
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VSHIFT)){
    assign("frac", "felec*d2inv + fvdw");
  }
  else if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "frac + fvdw*d2inv");
  }
  else if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "felec + fvdw*d2inv");
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "fvdw*d2inv");
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VSWITCH)){
    //nothing
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VSHIFT)){
    //nothing
    assign("frac", "fvdw;");
  }

}


void mknb_ene_force(){
  if(nb_func.elec_type==ESHIFT){
    assign("esh", "one - d2*c2inv");
  }

  //elec
  if(nb_func.cdie){
    if(nb_func.elec_type==ETRUNC){
      assign("felec", "q1*q2*dinv");
      assign("ecul", "ecul + felec");
    }
    else if(nb_func.elec_type==ESHIFT){
      assign("felec", "q1*q2*dinv*esh");
      assign("ecul", "ecul + felec*esh");
      assign("felec", "felec*(esh*d2inv+four*c2inv)");
    }
  }
  else {
    if(nb_func.elec_type==ETRUNC){
      assign("felec", "q1*q2*d2inv");
      assign("ecul", "ecul + felec");
      assign("felec", "two*felec");
    }
    else if(nb_func.elec_type==ESHIFT){
      assign("felec", "q1*q2*d2inv");
      assign("ecul", "ecul + felec*esh*esh");
      assign("felec", "two*felec*esh*(esh*d2inv+two*c2inv)");
    }
  }

  //vdw
  if(nb_func.vdw_type==VTRUNC){
    assign("nb12", "A*d12inv");
    assign("nb6", "B*d6inv");
    assign("evdw", "evdw + nb12 - nb6");
    assign("fvdw", "twelve*nb12 - six*nb6");
  }
  else if(nb_func.vdw_type==VSHIFT){
    assign("vsh","one - d6*c6inv");
    assign("evdw", "evdw + A*(d12inv - c12inv*(one + vsh+vsh)) - B*(d6inv -c6inv*(one + vsh))");
    assign("nb6", "six*B*(d6inv-d6*c12inv);");
    assign("nb12", "twelve*A*(d12inv-d6*c18inv)");
    assign("fvdw", "(nb12 - nb6)*d2inv");
  }

  {
    //SWITCH ....
    if((nb_func.elec_type==ESWITCH)||(nb_func.vdw_type==VSWITCH)){
      //elec
      if(nb_func.cdie&&(nb_func.elec_type==ESWITCH)){
	assign("eelec","q1*q2*dinv");
	assign("felec","eelec");
      }
      if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH)){
	assign("eelec","q1*q2*d2inv"); 
	assign("felec", "two*q1*q2*d2inv");
      }
      if(nb_func.vdw_type==VSWITCH){
	assign("nb6", "B*d6inv");
	assign("nb12","A*d12inv");
	assign("enb", "nb12-nb6;");
	assign("fvdw", "twelve*nb12-six*nb6");
      }
      //cumulate
      if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSWITCH)){
	assign("frac", "(felec+fvdw)*d2inv");
      }
      else if((nb_func.elec_type==ESWITCH)&&(!(nb_func.vdw_type==VSWITCH))){
	assign("frac", "felec*d2inv");
      }
      else if((!(nb_func.elec_type==ESWITCH))&&((nb_func.vdw_type==VSWITCH))){
	assign("frac", "fvdw*d2inv");
      }
      fprintf(mknb_output,"if(d2>c_on2){\n");
      //compute sw sw1
      assign("sh1", "c2-d2");
      assign("sh2", "(c2+two*d2-three*c_on2)");
      assign("sw", "sh1*sh2");
      assign("sw1", "four*(sh1-sh2)");

      if(!(nb_func.vdw_type==VSWITCH)){
	if(nb_func.cdie&&(nb_func.elec_type==ESWITCH)){
	  assign("frac", "frac*sw - felec*sw1;");
	}
	if((!nb_func.cdie)&&(nb_func.elec_type==ESWITCH)){
	  assign("frac", "frac*sw - 0.5*felec*sw1;");
	}
      }
      if(!(nb_func.elec_type==ESWITCH)){
	if(nb_func.vdw_type==VSWITCH)
	  assign("frac", "frac*sw - (A*d12inv-B*d6inv)*sw1");
      }
      //switch switch
      if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSWITCH)){
	if(nb_func.cdie)
	  assign("frac", "frac*sw - (A*d12inv-B*d6inv+felec)*sw1");
	else
	  assign("frac", "frac*sw - (A*d12inv-B*d6inv+0.5*felec)*sw1");
      }
      assign("sh1", "c2diffinv*sh1");
      assign("frac", "frac*sh1");
      assign("sw", "sw*sh1");
      if(nb_func.elec_type==ESWITCH){
	assign("eelec", "eelec*sw");
      }
      if(nb_func.vdw_type==VSWITCH){
	assign("enb", "enb*sw");
      }
      fprintf(mknb_output,"}\n");            
      if(nb_func.elec_type==ESWITCH){
	assign("ecul", "ecul + eelec");
      }
      if(nb_func.vdw_type==VSWITCH){
	assign("evdw", "evdw+enb");
      }
    }
  }
  //frac contains scalar forces of switch components so need to add shifted  
  if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VSHIFT)){
    assign("frac", "(felec+fvdw)");
  }
  else if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VSHIFT))
    assign("frac", "frac + fvdw");
  else if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VSWITCH))
    assign("frac", "frac + felec");
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "(felec+fvdw)*d2inv");
  }
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VSWITCH)){
    assign("frac", "frac + felec*d2inv");
  }
  else if((nb_func.elec_type==ETRUNC)&&(nb_func.vdw_type==VSHIFT)){
    assign("frac", "felec*d2inv + fvdw");
  }
  else if((nb_func.elec_type==ESWITCH)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "frac + fvdw*d2inv");
  }
  else if((nb_func.elec_type==ESHIFT)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "felec + fvdw*d2inv");
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VTRUNC)){
    assign("frac", "fvdw*d2inv");
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VSWITCH)){
    //nothing
  }
  else if((nb_func.elec_type==ENONE)&&(nb_func.vdw_type==VSHIFT)){
    //nothing
    assign("frac", "fvdw;");
  }

}



void mknb_force_body(){
  if(nb_func.kind==ENE)
    mknb_ene();

  if(nb_func.kind==FORCE)
    mknb_force();  

  if(nb_func.kind==ENE_FORCE)
    mknb_ene_force();  
}

void mknb_assign_first(){
  assign(array("force","fpos"),array("force","fpos"),"+","fracX");
  assign(array("force","fpos+1"),array("force","fpos+1"),"+","fracY");
  assign(array("force","fpos+2"),array("force","fpos+2"),"+","fracZ");
  if(nb_func.calc_vir){
    assign(array("fvir","bpos"),array("fvir","bpos"),"+","fracX");
    assign(array("fvir","bpos+1"),array("fvir","bpos+1"),"+","fracY");
    assign(array("fvir","bpos+2"),array("fvir","bpos+2"),"+","fracZ");
  }

}


void mknb_assign_second(){
  assign("tx","frac*dx");
  assign("ty","frac*dy");
  assign("tz","frac*dz");
  assign("fracX","fracX - tx");
  assign("fracY","fracY - ty");
  assign("fracZ","fracZ - tz");
  assign(array("force","spos"),array("force","spos"),"+","tx");
  assign(array("force","spos+1"),array("force","spos+1"),"+","ty"); 
  assign(array("force","spos+2"),array("force","spos+2"),"+","tz");  
}


void mknb_compute_force(){
  if(CUTOFF){
    fprintf(mknb_output,"if(d2<c2){\n");
  }

  mknb_load_par_second();
  mknb_d_powers();
  mknb_force_body();
  if(nb_func.kind!=ENE)
    mknb_assign_second();


  if(CUTOFF){
    fprintf(mknb_output,"}\n");
  }

}

void mknb_inner_loop(){
  mknb_start_loop("j","b1","b2");
  {
    fprintf(mknb_output,"{\n");

    mknb_comment("Load coor second");
    mknb_load_second();

    mknb_comment("Compute distances");
    mknb_compute_dist();

    mknb_comment("Compute force");
    mknb_compute_force();
    fprintf(mknb_output,"}\n");
  }
}


void mknb_outer_loop() {

  mknb_comment("Start outer loop over neighborlists");
  mknb_start_loop("i","0","natoms");
  {
    fprintf(mknb_output,"{\n");
    mknb_comment("Initialize outer atom");
    mknb_load_first();

    mknb_comment("Load limits for loop over neighbors");
    mknb_inner_boundary();

    if(nb_func.kind!=ENE){
      mknb_comment("Zero forces");
      mknb_zero_forces(); 
    }

    mknb_comment("Load pars");
    mknb_load_par_first();


    mknb_comment("Inner loop");
    mknb_inner_loop();

    if(nb_func.kind!=ENE){
      mknb_comment("Assign to first");
      mknb_assign_first();
    }
    fprintf(mknb_output,"}\n");
  }
}

void generate_function(const char * name){
  static int f = 1;
  {
    printf("\r%i Generating function %s\r",f++,name);
    fflush(stdout);
  }
  //generate function
  mknb_force_function_header(name);

  //init body
  {
    fprintf(mknb_output,"{\n");
    
    mknb_force_function_internal_vars();
    mknb_outer_loop();
    //end body
    if(CENE){
      assign("*ene_vdw","*ene_vdw + evdw");
      if(nb_func.elec_type!=ENONE)
	assign("*ene_cul","*ene_cul+ ecul");
    }
    fprintf(mknb_output,"\n}\n");
  }
}

int main(){

  //Output files

  //ENE files
  FILE * fene_c = fopen("enefunc.c","w");
  mknb_output = fene_c;
  mknb_file_header("enefunc.c");
  fprintf(mknb_output,"#define COOR_DIM 4\n");
  fprintf(mknb_output,"#include <math.h>\n");
  fprintf(mknb_output,"#include <almost.h>\n");


  FILE * fene_h = fopen("enefunc.h","w");
  mknb_output = fene_h;
  mknb_file_header("enefunc.h");
  fprintf(mknb_output,"extern \"C\" {\n");

  FILE * fene_force_c = fopen("eneforcefunc.c","w");
  mknb_output = fene_force_c;
  mknb_file_header("eneforcefunc.c");
  fprintf(mknb_output,"#define COOR_DIM 4\n");
  fprintf(mknb_output,"#include <math.h>\n");
  fprintf(mknb_output,"#include <almost.h>\n");

  FILE * fene_force_h = fopen("eneforcefunc.h","w");
  mknb_output = fene_force_h;
  mknb_file_header("eneforcefunc.h");
  fprintf(mknb_output,"extern \"C\" {\n");


  FILE * fforce_c = fopen("forcefunc.c","w");
  mknb_output = fforce_c;
  mknb_file_header("forcefunc.c");
  fprintf(mknb_output,"#define COOR_DIM 4\n");
  fprintf(mknb_output,"#include <math.h>\n");
  fprintf(mknb_output,"#include <almost.h>\n");

  FILE * fforce_h = fopen("forcefunc.h","w");
  mknb_output = fforce_h;
  mknb_file_header("forcefunc.h");
  fprintf(mknb_output,"extern \"C\" {\n");


  nb_func.ucc = 1;
  nb_func.float_type=DOUBLE;
  nb_func.elec_type = ENONE;
  nb_func.cdie = 0;
  nb_func.kind = ENE;
  nb_func.geom = 0;
  nb_func.bCompat = 0;

  int GEO_KIND[2] = { PLAIN, PBC};
  int VDW_CUT[3] = { VTRUNC, VSHIFT, VSWITCH };
  int ELE_CUT[3] = { ETRUNC, ESHIFT, ESWITCH };
  int ELE_FUN[3] = { ENONE, CDIE, RDIE };
  int FUN_KIND[3] = { ENE, FORCE, ENE_FORCE };

  FILE * HEADER[3] = { fene_h,fforce_h,fene_force_h};
  FILE * CODE[3] = { fene_c,fforce_c,fene_force_c};

  int FTYPE[2] = { DOUBLE, FLOAT};
  {
    printf("Generating headers .......\n");
    //    mknb_output = fopen("nbkernel/almnbkernel.h","w");
    //    mknb_file_header("almnbkernel.h");
    //elec NONE
    for(int g = 0;g<2;g++){
      nb_func.geom = GEO_KIND[g];

      for(int k = 0; k < 3 ; k ++){
	nb_func.kind = FUN_KIND[k];
	nb_func.elec_type = ENONE;
	nb_func.cdie=0;
	mknb_output = HEADER[k];
	for(int vc = 0; vc < 3; vc++){ 
	  nb_func.vdw_type = VDW_CUT[vc];
	  for(int ft = 0; ft < 2; ft++){
	    nb_func.float_type= FTYPE[ft];
	    //generate function
	    nb_func.calc_vir =0;
	    string name = spec2name();
	    mknb_func_desc();
	    mknb_force_function_header(name.c_str());
	    fprintf(mknb_output,";\n\n");
	    for(int v=1;v<2;v++){
	      if(nb_func.geom==PLAIN) continue;
	      nb_func.calc_vir = v;
	      name = spec2name();
	      mknb_func_desc();
	      mknb_force_function_header(name.c_str());
	      fprintf(mknb_output,";\n\n");
	    }
	  }
	}
      }
    }
    

   
    for(int g = 0;g<2;g++){
      nb_func.geom = GEO_KIND[g];
    
      for(int k = 0; k < 3 ; k ++){
	nb_func.kind = FUN_KIND[k];
	mknb_output = HEADER[k];
	for(int et = 1; et < 3 ; et++){
	  if(ELE_FUN[et]==CDIE) nb_func.cdie =1;
	  else nb_func.cdie = 0;
	  for(int ec=0; ec<3; ec++){
	    nb_func.elec_type=ELE_CUT[ec];
	    for(int vc = 0; vc < 3; vc++){ 
	      nb_func.vdw_type = VDW_CUT[vc];
	      for(int ft = 0; ft < 2; ft++){
		nb_func.float_type= FTYPE[ft];
		//generate function
		nb_func.calc_vir =0;
		string name = spec2name();
		mknb_func_desc();
		mknb_force_function_header(name.c_str());
		fprintf(mknb_output,";\n\n");
		for(int v=1;v<2;v++){
		  if(nb_func.geom==PLAIN) continue;
		  nb_func.calc_vir = v;
		  name = spec2name();
		  mknb_func_desc();
		  mknb_force_function_header(name.c_str());
		  fprintf(mknb_output,";\n\n");
		}
	      }
	    }
	  }
	}
      }
      
    }
  }

  {  

    printf("Generating functions .......\n");

    //    mknb_output = fopen("nbkernel/almnbkernel.c","w");

    //    fprintf(mknb_output,"#include <math.h>\n");
    //    fprintf(mknb_output,"#ifndef COOR_DIM\n#define COOR_DIM 4\n#endif\n");
    //    fprintf(mknb_output,"double invsqrt(double);\n");
    //    fprintf(mknb_output,"float invsqrtf(float);\n");
    //elec NONE
    
    for(int g = 0;g<2;g++){
      nb_func.geom = GEO_KIND[g];

      for(int k = 0; k < 3 ; k ++){
	nb_func.kind = FUN_KIND[k];
	mknb_output = CODE[k];
	nb_func.elec_type = ENONE;
	nb_func.cdie=0;
	for(int vc = 0; vc < 3; vc++){ 
	  nb_func.vdw_type = VDW_CUT[vc];
	  for(int ft = 0; ft < 2; ft++){
	    nb_func.float_type= FTYPE[ft];
	    //generate function
	    nb_func.calc_vir =0;
	    string name = spec2name();
	    mknb_func_desc();
	    generate_function(name.c_str());
	    for(int v=1;v<2;v++){
	      if(nb_func.geom==PLAIN) continue;
	      nb_func.calc_vir = v;
	      name = spec2name();
	      mknb_func_desc();
	      generate_function(name.c_str());
	    }
	  }
	}
      }

      
    }

    
    for(int g = 0;g<2;g++){
      nb_func.geom = GEO_KIND[g];

      for(int k = 0; k < 3 ; k ++){
	nb_func.kind = FUN_KIND[k];
	mknb_output = CODE[k];
	for(int et = 1; et < 3 ; et++){
	  if(ELE_FUN[et]==CDIE) nb_func.cdie =1;
	  else nb_func.cdie = 0;
	  for(int ec=0; ec<3; ec++){
	    nb_func.elec_type=ELE_CUT[ec];
	    for(int vc = 0; vc < 3; vc++){ 
	      nb_func.vdw_type = VDW_CUT[vc];
	      for(int ft = 0; ft < 2; ft++){
		nb_func.float_type= FTYPE[ft];
		//generate function
		nb_func.calc_vir =0;
		string name = spec2name();
		mknb_func_desc();
		generate_function(name.c_str());
		for(int v=1;v<2;v++){
		  if(nb_func.geom==PLAIN) continue;
		  nb_func.calc_vir = v;
		  name = spec2name();
		  mknb_func_desc();
		  generate_function(name.c_str());
		}
	      }
	    }
	  }
	  
	}
      }
    }
    //    fclose(mknb_output);
  }

  //
  fclose(fene_force_c);
  fprintf(fene_force_h,"}\n");
  fclose(fene_force_h);

  fclose(fforce_c);
  fprintf(fforce_h,"}\n");
  fclose(fforce_h);

  fclose(fene_c);
  fprintf(fene_h,"}\n");
  fclose(fene_h);




  printf("\n\n");
  fprintf(stderr,"//Optimize: (energy) (c2 - d2) and one - d2 * c2inv in mixed shift, switch\n");
  fprintf(stderr,"//Optimize: (force) c12inv*d6 and c18inv*d6\n");
  fprintf(stderr,"//Optimize: (force) switch switch\n");
  fprintf(stderr,"//Optimize: (force) shift shift\n");
  fprintf(stderr,"//Optimize: (ene force) shift shift\n");
  printf("\n"); 
}