/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under 
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Paul Crozier (SNL)
   Changed in September 2012 by Byshkin Maxim (department of Chimestry, University of Salerno,It. email maxbysh@gmail.com) 
------------------------------------------------------------------------- */

#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_gcmc.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "group.h"
#include "domain.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
//#include "bond.h"
//#include  "neigh_list.h"
//#include "neigh_request.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
#define PI 3.14159265358979

/* ---------------------------------------------------------------------- */

FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{

  if (narg < 11) error->all(FLERR,"Illegal fix GCMC command");

  vector_flag = 1;
  size_vector = 6;
  global_freq = 1;
  extvector = 0;
  restart_global = 1;
  time_depend = 1;

  // required args
  
  nevery = atoi(arg[3]);
  nexchanges = atoi(arg[4]);
  nmcmoves = atoi(arg[5]);
  ntype = atoi(arg[6]);
  seed = atoi(arg[7]);
  reservoir_temperature = atof(arg[8]);
  chemical_potential = atof(arg[9]);
  displace = atof(arg[10]);

  if (ntype <= 0 || ntype > atom->ntypes) 
    error->all(FLERR,"Invalid atom type in fix GCMC command");
  if (nexchanges < 0) error->all(FLERR,"Illegal fix GCMC command");
  if (nmcmoves < 0) error->all(FLERR,"Illegal fix GCMC command");
  if (seed <= 0) error->all(FLERR,"Illegal fix GCMC command");
  if (reservoir_temperature < 0.0) error->all(FLERR,"Illegal fix GCMC command");  
  if (displace < 0.0) error->all(FLERR,"Illegal fix GCMC command"); 

  // compute beta, lambda, sigma, and the zz factor
  
  beta = 1.0/(force->boltz*reservoir_temperature);
  double gas_mass = atom->mass[ntype];
  double lambda = sqrt(force->hplanck*force->hplanck/
		       (2.0*MY_PI*gas_mass*force->mvv2e*
			force->boltz*reservoir_temperature));
  sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
  zz = exp(beta*chemical_potential)/(pow(lambda,3));

  // set defaults

  molflag = 1;

  // read options from end of input line

  options(narg-11,&arg[11]);

  // random number generator, same for all procs

  random_equal = new RanPark(lmp,seed);

  // random number generator, not the same for all procs

  random_unequal = new RanPark(lmp,seed);
  
  // compute the number of MC cycles that occur nevery timesteps
  
  ncycles = nexchanges + nmcmoves;
  
  // set up reneighboring

  force_reneighbor = 1;
  next_reneighbor = update->ntimestep + 1;

  // zero out counters
  
  nmove_attempts = 0.0;   
  nmove_successes = 0.0;  
  ndel_attempts = 0.0;    
  ndel_successes = 0.0;   
  ninsert_attempts = 0.0; 
  ninsert_successes = 0.0;

  nmax = 0;
  local_gas_list = NULL;

}

/* ---------------------------------------------------------------------- */

FixGCMC::~FixGCMC()
{
  delete random_equal;
  delete random_unequal;
  memory->sfree(local_gas_list);
}

/* ---------------------------------------------------------------------- */

int FixGCMC::setmask()
{
  int mask = 0;
  mask |= PRE_EXCHANGE;
  return mask;
}

/* ---------------------------------------------------------------------- */

void FixGCMC::init()
{
  // check that no deletable atoms are in atom->firstgroup
  // deleting such an atom would not leave firstgroup atoms first

  int *type = atom->type;
  
  if (atom->firstgroup >= 0) {
    int *mask = atom->mask;
    int nlocal = atom->nlocal;
    int firstgroupbit = group->bitmask[atom->firstgroup];

    int flag = 0;
    for (int i = 0; i < nlocal; i++)
      if ((type[i] == ntype) && (mask[i] && firstgroupbit)) flag = 1;

    int flagall;
    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);

    if (flagall)
      error->all(FLERR,"Cannot do GCMC on atoms in atom_modify first group");
  }

  // if molflag not set, warn if any deletable atom has a mol ID

  if (molflag == 1 && atom->molecule_flag) {
    int *molecule = atom->molecule;
    int *mask = atom->mask;
    int nlocal = atom->nlocal;
    int flag = 0;
    for (int i = 0; i < nlocal; i++)
      if (type[i] == ntype)
        if (molecule[i]) flag = 1;
    int flagall;
    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
    if (flagall && comm->me == 0)
      error->warning(FLERR,"Fix GCMC may delete atom with non-zero molecule ID");
  }

//  if (molflag && atom->molecule_flag == 0)
 //     error->all(FLERR,"Fix GCMC molecule command requires atom attribute molecule");
      
 // if (molflag != 0) error->all(FLERR,"Fix GCMC molecule feature does not yet work"); 
  
  if (force->pair->single_enable == 0) 
    error->all(FLERR,"Fix GCMC incompatible with given pair_style");
}

/* ----------------------------------------------------------------------
   attempt particle insertions and deletions
   done before exchange, borders, reneighbor
   so that ghost atoms and neighbor lists will be correct
------------------------------------------------------------------------- */

void FixGCMC::pre_exchange()
{
  // just return if should not be called on this timestep
  int *num_bond = atom->num_bond;
  int **bond_atom = atom->bond_atom;

 if (next_reneighbor != update->ntimestep) return;

  if (domain->triclinic == 0) {
    xlo = domain->boxlo[0];
    xhi = domain->boxhi[0];
    ylo = domain->boxlo[1];
    yhi = domain->boxhi[1];
    zlo = domain->boxlo[2];
    zhi = domain->boxhi[2];
    sublo = domain->sublo;
    subhi = domain->subhi;
  } else {
    xlo = domain->boxlo_bound[0];
    xhi = domain->boxhi_bound[0];
    ylo = domain->boxlo_bound[1];
    yhi = domain->boxhi_bound[1];
    zlo = domain->boxlo_bound[2];
    zhi = domain->boxhi_bound[2];
    sublo = domain->sublo_lamda;
    subhi = domain->subhi_lamda;
  } 
  
  volume = domain->xprd * domain->yprd * domain->zprd;
  
  // grow local_gas_list array if necessary

  if (atom->nlocal > nmax) {
    memory->sfree(local_gas_list);
    nmax = atom->nmax;
    local_gas_list = (int *) memory->smalloc(nmax*sizeof(int),
					     "GCMC:local_gas_list");
  }
  //molflag- number of types atoms in molecule
  //ntypes  -  ntype+molflag-1 - gas inserted interaction which interact
  //ntype+molflag    -  ntypes+2*molflag-1 gas inserted interaction which donot interact 
  int *type = atom->type;
  ngas_local = 0;
  for (int i = 0; i < atom->nlocal; i++)
    if (type[i] == ntype)
      local_gas_list[ngas_local++] = i;
                                           
  MPI_Allreduce(&ngas_local,&ngas,1,MPI_INT,MPI_SUM,world);
  MPI_Scan(&ngas_local,&ngas_before,1,MPI_INT,MPI_SUM,world);
  ngas_before -= ngas_local;
  ngas_local = 0;
  for (int i = 0; i < atom->nlocal; i++)
    if (atom->type[i]>=ntype&&atom->type[i]<ntype+molflag*2)//molecule[i]>0
      local_gas_list[ngas_local++] = i;
   int maxt=0;
   int mint=1000000;
   int nind=1;
    for(int i=1;i<ngas_local;i++)
	{ 
		if(atom->type[local_gas_list[i]]>maxt)
			maxt=atom->type[local_gas_list[i]];
		if(atom->type[local_gas_list[i]]<mint)
			mint=atom->type[local_gas_list[i]];
		if(atom->molecule[local_gas_list[i]]==atom->molecule[local_gas_list[0]])
			nind++;
	}
	int chgr=0;
        double qq;

    printf("ngas %i ngas_local %i\n",ngas, ngas_local); 
  // perform ncycles MC cycles
   printf("atoms loc%i",atom->nlocal);
  printf("atoms ngh%i\n",atom->nghost);
  printf("gas types from %i to %i,   atoms in malecule %i \n",mint,maxt,nind);
 
  for (int i = 0; i < ncycles; i++) {
    int random_int_fraction = 
      static_cast<int>(random_equal->uniform()*ncycles) + 1;	
    if (random_int_fraction <= nmcmoves) {
      attempt_move();
    } else {
      if (random_equal->uniform() < 0.5) attempt_deletion();
      else attempt_insertion();
    }
  }
  
  next_reneighbor = update->ntimestep + nevery;
} 

/* ----------------------------------------------------------------------
   choose particle randomly across all procs and attempt displacement
------------------------------------------------------------------------- */

void FixGCMC::attempt_move()//was not changed and should not be called with molecules
{  
  int i,iwhichglobal,iwhichlocal;
  double rx,ry,rz;
  double coord[3];
  double **x = atom->x;

  nmove_attempts += 1.0;
  
  int success = 0;
  iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
  if ((iwhichglobal >= ngas_before) && 
      (iwhichglobal < ngas_before + ngas_local)) {
    iwhichlocal = iwhichglobal - ngas_before;
    i = local_gas_list[iwhichlocal];
    double energy_before = energy(ntype,i,x[i]);
    double rsq = 1.1;
    while (rsq > 1.0) {
      rx = 2*random_unequal->uniform() - 1.0;
      ry = 2*random_unequal->uniform() - 1.0;
      if (domain->dimension == 3) rz = 2*random_unequal->uniform() - 1.0;
      else rz = 0.0;
      rsq = rx*rx + ry*ry + rz*rz;
    }
    coord[0] = x[i][0] + displace*rx;
    coord[1] = x[i][1] + displace*ry;
    coord[2] = x[i][2] + displace*rz;
    double energy_after = energy(ntype, i,coord);
    if (random_unequal->uniform() < exp(-beta*(energy_after - energy_before))) {
      x[i][0] = coord[0]; 
      x[i][1] = coord[1]; 
      x[i][2] = coord[2];
      success = 1;
    }
  } 
  
  int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
    nmove_successes += 1.0;  
    comm->borders(); 
  }
  
}

void FixGCMC::attempt_insertion()
{  
  int flag,success;
  int i,iwhichglobal,iwhichlocal;
  double **x = atom->x;
  double lamda[3];
  double *newcoord;
  double coord[10][3];
  int **bond_atom = atom->bond_atom;
  int *num_bond = atom->num_bond;
  int indxs[10];
  ninsert_attempts += 1.0;
  int m;
  //find index of gas in switch on
  for(m=0;m<ngas_local;m++)
	  if(atom->type[local_gas_list[m]]==ntype+molflag)
		  break;
//  printf("gas found");
  if(m==ngas_local)
		  error->all(FLERR,"not enaugh gas for GCMC");
  m=local_gas_list[m];
  coord[0][0] = xlo + random_equal->uniform() * (xhi-xlo);
  coord[0][1] = ylo + random_equal->uniform() * (yhi-ylo);
  coord[0][2] = zlo + random_equal->uniform() * (zhi-zlo);
  //attempt inset m-th atom to this coordinate
  if (domain->triclinic) {
    domain->x2lamda(coord[0],lamda);
    newcoord = lamda;
  } else newcoord = coord[0];
  indxs[0]=m;
  int nind=1;
  int nall = atom->nlocal;// + atom->nghost;
  for(int inh=0;inh<ngas_local;inh++)
  {
	  int in=local_gas_list[inh];
	  if(atom->molecule[in]==atom->molecule[m])
	  {
		  if (in==m) continue;
		  indxs[nind]=in;
		  nind++;
	  }
  }
//  if(nind!=molflag)
//	  printf("nind %i ",nind);
  for (int in=1; in < nind; in++) {
	  int ih=indxs[in];
		coord[in][0]=x[ih][0]-x[m][0]+coord[0][0];
		coord[in][1]=x[ih][1]-x[m][1]+coord[0][1];
		coord[in][2]=x[ih][2]-x[m][2]+coord[0][2];
	}
        double bRMAT[4][4];
	    double C1 = random_unequal->uniform() * PI;
        double C2 = random_unequal->uniform() * 2*PI;
        double C3 = random_unequal->uniform() * 2*PI;
        double CC1 = cos(C1);
        double SC1 = sin(C1);
        double CC2 = cos(C2);
        double SC2 = sin(C2);
        double CC3 = cos(C3);
        double SC3 = sin(C3);
        bRMAT[1][1] = CC2*CC3 - SC2*SC3*CC1;
        bRMAT[2][1] = SC2*CC3 + CC2*SC3*CC1;
        bRMAT[3][1] = SC1*SC3;
        bRMAT[1][2] = -CC2*SC3 - SC2*CC3*CC1;
        bRMAT[2][2] = -SC2*SC3 + CC2*CC3*CC1;
        bRMAT[3][2] = SC1*CC3;
        bRMAT[1][3] = SC1*SC2;
        bRMAT[2][3] = -SC1*CC2;
        bRMAT[3][3] = CC1;
		float XP0[4];
		float XP[4];
		for(int in=1;in<nind;in++)
		{
	    XP0[1]=coord[in][0]-coord[0][0];
		XP0[2]=coord[in][1]-coord[0][1];
		XP0[3]=coord[in][2]-coord[0][2];
    	//random rotation of molecule
        for(int I=1;I<4;I++)
			XP[I] = bRMAT[I][1]*XP0[1] + bRMAT[I][2]*XP0[2] + bRMAT[I][3]*XP0[3];
		coord[in][0]=coord[0][0]+XP[1];
		coord[in][1]=coord[0][1]+XP[2];
		coord[in][2]=coord[0][2]+XP[3];
		}

  flag = 1;
  for(int ih=0;ih<nind;ih++)
  if (!(coord[ih][0] >= sublo[0] && coord[ih][0] < subhi[0] &&
    coord[ih][1] >= sublo[1] && coord[ih][1] < subhi[1] &&
    coord[ih][2] >= sublo[2] && coord[ih][2] < subhi[2])) flag = 0;

  success = 0;
 


  if (flag) {
    double insertion_energy =0;// energy(ntype,m,coord[0]);
	for(int in=0;in<nind;in++)
{
		int ih=indxs[in];
                    atom->q[ih]*=1e+10;
		insertion_energy+=energy((atom->type[ih])-molflag,ih,coord[in]);
                    atom->q[ih]/=1e+10;
       }

    if (random_unequal->uniform() <
        zz*volume*exp(-beta*insertion_energy)/(ngas+1)) {
	  for(int in=0;in<nind;in++)
	
	{
		int ih=indxs[in];
		atom->type[ih]-=molflag;
                atom->q[ih]*=1e+10;
		atom->v[ih][0] = random_unequal->gaussian()*sigma;
		atom->v[ih][1] = random_unequal->gaussian()*sigma;
		atom->v[ih][2] = random_unequal->gaussian()*sigma;
		x[ih][0] = coord[in][0]; 
		x[ih][1] = coord[in][1]; 
		x[ih][2] = coord[in][2];
	  }          
      success = 1;
      printf("inserted atom of %i atoms, insertion energy %f",nind, insertion_energy);
    }
  }

   int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
    ninsert_successes += 1.0;  
    ngas++;
    printf(" now inserted %i\n",ngas);

    comm->borders(); 
  }
  
}
/* ----------------------------------------------------------------------
   attempt particle deletion
------------------------------------------------------------------------- */

void FixGCMC::attempt_deletion()
{ 
  ndel_attempts += 1.0;
  if (ngas == 0) return;
  int **bond_atom = atom->bond_atom;
  int *num_bond = atom->num_bond;
  int indxs[10];
  int i,iwhichglobal,iwhichlocal;
  AtomVec *avec = atom->avec;
  int m;
  int ms=0;
  do
  { 
	  ms++;
	  if(ms%10000==0) printf("nema gasu, ngas=%i",ngas);
	  m = static_cast<int> (ngas_local*random_equal->uniform());
  }
  while (atom->type[local_gas_list[m]]!=ntype);


  //find index of gas in switch off
   m=local_gas_list[m];
   indxs[0]=m;
  int nind=1;
  int nall = atom->nlocal;
  for(int inh=0;inh<ngas_local;inh++)
  {
	  int in=local_gas_list[inh];
	  if(atom->molecule[in]==atom->molecule[m])
	  { 
		  if (in==m) continue;
		  indxs[nind]=in;
		  nind++;
   //                printf("attempt deln %i %i",in,atom->type[in]);
	  }
   }
	  
  // choose particle randomly across all procs and delete it
  // keep ngas, ngas_local, ngas_before, and local_gas_list current
  // after each deletion
 //  if( nind!=molflag)
//	   printf("gcMCerror");
    
    int success = 0;
  

    double deletion_energy = 0;//energy(ntype,m,atom->x[m]);      
	for(int in=0;in<nind;in++)
	{
		int ih=indxs[in];
 		deletion_energy+=energy(atom->type[ih],ih,atom->x[ih]);
	}

    if (random_unequal->uniform() < 
	ngas*exp(beta*deletion_energy)/(zz*volume)) {
	  for(int in=0;in<nind;in++)
            {
		  atom->type[indxs[in]]+=molflag;
                  atom->q[indxs[in]]/=1e+10;
            } 
      printf("deleted  atom of %i atoms, deletion enery %f",nind,deletion_energy);
      success = 1; 
    }
  
  int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
     ndel_successes += 1.0;
     ngas--;
     comm->borders();
    printf("now inserted %i\n",ngas);

   }
}

/* ----------------------------------------------------------------------
   attempt particle insertion
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   compute particle's interaction energy with the rest of the system
------------------------------------------------------------------------- */

double FixGCMC::energy(int ity, int i, double *coord) 
{
  double delx,dely,delz,rsq;

  double **x = atom->x;
  int *type = atom->type;
  int nall = atom->nlocal + atom->nghost;
  pair = force->pair;
  cutsq = force->pair->cutsq;

  double fpair = 0.0;
  double factor_coul = 1.0;
  double factor_lj = 1.0;
  double phicoul=0;
  double total_energy = 0.0;
  for (int j = 0; j < nall; j++) {

   //   if (i == j) continue;
	int jtype = type[j];
	if((jtype>ntype+molflag-1)&&(jtype<ntype+2*molflag)) continue;//no interaction with gas not inserted
        if(atom->molecule[i]==atom->molecule[j]) continue;// no nonbonded iteraction inside molecule 
    delx = coord[0] - x[j][0];
    dely = coord[1] - x[j][1];
    delz = coord[2] - x[j][2];
    rsq = delx*delx + dely*dely + delz*delz;

    if (rsq < cutsq[ity][jtype])
      total_energy += 
	pair->single(i,j,ity,jtype,rsq,0,factor_lj,fpair);
    phicoul += force->qqrd2e * atom->q[i]*atom->q[j]/sqrt(rsq);
  }
  //printf("LJ  %f coul %f\n",total_energy, phicoul);
 total_energy += phicoul;
  return total_energy;
}

/* ----------------------------------------------------------------------
   parse optional parameters at end of input line 
------------------------------------------------------------------------- */

void FixGCMC::options(int narg, char **arg)
{
  if (narg < 0) error->all(FLERR,"Illegal fix GCMC command");

  int iarg = 0;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"molecule") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix GCMC command");
      if (strcmp(arg[iarg+1],"no") == 0) molflag = 1;
      else sscanf(arg[iarg+1],"%i",&molflag);//if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1;
      //else error->all(FLERR,"Illegal fix evaporate command");
      iarg += 2;
    } else error->all(FLERR,"Illegal fix GCMC command");
  }
  printf("%i atom types  in molecules",molflag);
}

/* ----------------------------------------------------------------------
  return acceptance ratios
------------------------------------------------------------------------- */

double FixGCMC::compute_vector(int n)
{
  if (n == 0) return nmove_attempts;
  if (n == 1) return nmove_successes;
  if (n == 2) return ndel_attempts;
  if (n == 3) return ndel_successes;
  if (n == 4) return ninsert_attempts;
  if (n == 5) return ninsert_successes;
  return 0.0;
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based arrays
------------------------------------------------------------------------- */

double FixGCMC::memory_usage()
{
  double bytes = nmax * sizeof(int);
  return bytes;
}

/* ----------------------------------------------------------------------
   pack entire state of Fix into one write 
------------------------------------------------------------------------- */

void FixGCMC::write_restart(FILE *fp)
{
  int n = 0;
  double list[4];
  list[n++] = random_equal->state();
  list[n++] = random_unequal->state();
  list[n++] = next_reneighbor;

  if (comm->me == 0) {
    int size = n * sizeof(double);
    fwrite(&size,sizeof(int),1,fp);
    fwrite(list,sizeof(double),n,fp);
  }
}

/* ----------------------------------------------------------------------
   use state info from restart file to restart the Fix  
------------------------------------------------------------------------- */

void FixGCMC::restart(char *buf)
{
  int n = 0;
  double *list = (double *) buf;

  seed = static_cast<int> (list[n++]);
  random_equal->reset(seed);
  
  seed = static_cast<int> (list[n++]);
  random_unequal->reset(seed);
  
  next_reneighbor = static_cast<int> (list[n++]);
}