/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   https://www.lammps.org/, Sandia National Laboratories
   LAMMPS development team: developers@lammps.org

   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: Shrikrushna Dumbre (IISERB)
------------------------------------------------------------------------- */

#include "pair_soft_cut.h"

#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "fix.h"
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "update.h"

#include <cmath>
#include <cstring>

using namespace LAMMPS_NS;
using namespace MathConst;

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

PairSoftCut::PairSoftCut(LAMMPS *lmp) : Pair(lmp)
{
	writedata=1;
}

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

PairSoftCut::~PairSoftCut()
{
	if(allocated){
		memory->destroy(setflag);
		memory->destroy(cutsq);
		memory->destroy(cut);
		memory->destroy(epsilon);
		memory->destroy(eta1);
		memory->destroy(eta2);
		memory->destroy(mu);
		memory->destroy(nu);
		memory->destroy(v0);
		memory->destroy(r_m);
		memory->destroy(offset);
	}
}

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

void PairSoftCut::allocate()
{
	allocated = 1;
	int n = atom->ntypes + 1;
	
	memory->create(setflag, n, n, "pair:setflag");
	for(int i =1; i<n; i++)
		for(int j = i; j<n; j++) setflag[i][j]=0;
	memory->create(cutsq, n, n, "pair:cutsq");
	memory->create(cut, n, n, "pair:cut");
	memory->create(epsilon, n, n, "pair:epsilon");
	memory->create(eta1, n, n, "pair:eta1");
	memory->create(eta2, n, n, "pair:eta2");
	memory->create(mu, n, n, "pair:mu");
	memory->create(nu, n, n, "pair:nu");
	memory->create(v0, n, n, "pair:v0");
	memory->create(r_m, n, n, "pair:r_m");
	memory->create(offset, n, n, "pair:offset");
}

/* ----------------------------------------------------------------------
   global settings
------------------------------------------------------------------------- */

void PairSoftCut::settings(int narg, char **arg)
{
	if(narg != 1) error->all(FLERR,"Illegal pair_style command");
	
	cut_global = utils::numeric(FLERR, arg[0], false, lmp);
	
	// reset cutoffs that have been explicitly set
	
	if(allocated){
		int i, j;
		for(i = 1; i <= atom->ntypes; i++)
			for(j = 1; j<=atom->ntypes; j++)
				if(setflag[i][j]) cut[i][j] = cut_global;
	}
}

/* ----------------------------------------------------------------------
   set coeffs for one or more type pairs
------------------------------------------------------------------------- */

void PairSoftCut::coeff(int narg, char **arg)
{
if (narg < 9 || narg >10) error->all(FLERR,"Incorrect arguments for pair coefficients");
if(!allocated)allocate();

int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);

double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double v0_one = utils::numeric(FLERR, arg[3], false, lmp);
double r_m_one = utils::numeric(FLERR, arg[4], false, lmp);
double eta1_one = utils::numeric(FLERR, arg[5], false, lmp);
double eta2_one = utils::numeric(FLERR, arg[6], false, lmp);
double mu_one = utils::numeric(FLERR, arg[7],false, lmp);
double nu_one = utils::numeric(FLERR, arg[8], false, lmp);
double cut_one = cut_global;
if(narg ==11) cut_one = utils::numeric(FLERR, arg[9], false, lmp);

int count = 0;
for (int i = ilo; i<= ihi; i++){
	for(int j = MAX(jlo, i); j<= jhi; j++){
		epsilon[i][j] = epsilon_one;
		v0[i][j] = v0_one;
		r_m[i][j] = r_m_one;
		eta1[i][j] = eta1_one;
		eta2[i][j] = eta2_one;
		mu[i][j] = mu_one;
		nu[i][j] = nu_one;
		cut[i][j] = cut_one;
		setflag[i][j] = 1;
		count++;
		}
	}
	if(count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}

/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

double PairSoftCut::init_one(int i, int j)
{
	if(setflag[i][j]==0) error->all(FLERR, "All pair coeffs are not set");
	
	if(offset_flag && cut[i][j]<r_m[i][j]){
		double dr = (cut[i][j]/r_m[i][j]);
		double a1 = pow(dr,eta1[i][j]);
		offset[i][j] = v0[i][j]*pow(1-a1,eta2[i][j]); 
	}
	else
	offset[i][j] = 0.0;
	
	epsilon[j][i] = epsilon[i][j];
	v0[j][i] = v0[i][j];
	r_m[j][i] = r_m[i][j];
	eta1[j][i] = eta1[i][j];
	eta2[j][i] = eta2[i][j];
	mu[j][i] = mu[i][j];
	nu[j][i] = nu[i][j];
	offset[j][i] = offset[i][j];
	
	return cut[i][j];
}

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

void PairSoftCut::compute(int eflag, int vflag)
{
	int i, j, ii, jj, inum, jnum, itype, jtype;
	double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
	double dr, r, rsq, a1, a2, a3, forcesoft, factor_soft;
	int *ilist, *jlist, *numneigh, **firstneigh;
	
	evdwl = 0.0;
	ev_init(eflag, vflag);
	
	double **x = atom->x;
	double **f = atom->f;
	int *type = atom->type;
	int nlocal = atom->nlocal;
	double *special_soft = force->special_lj;
	int newton_pair = force->newton_pair;
	
	inum = list->inum;
	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;
	
	//loop over neighbour of my atoms
	
	for(ii=0; ii< inum; ii++){
		
		i=ilist[ii];
		xtmp = x[i][0];
		ytmp = x[i][1];
		ztmp = x[i][2];
		itype = type[i];
		jlist = firstneigh[i];
		jnum = numneigh[i];
		
		for (jj=0; jj < jnum; jj++){
			j = jlist[jj];
			factor_soft = special_soft[sbmask(j)];
			j &=NEIGHMASK;
			
			delx = xtmp - x[j][0];
			dely = ytmp - x[j][1];
			delz = ztmp - x[j][2];
			rsq = delx * delx + dely * dely + delz * delz;
			jtype = type[j];
			
			if(rsq < cutsq[itype][jtype]){
				r = sqrt(rsq);
				dr = r/r_m[itype][jtype];
				a1 = 1 - pow(dr,eta1[itype][jtype]);
				a2 = pow(dr, eta1[itype][jtype]-1);
				forcesoft = (v0[itype][jtype]*eta1[itype][jtype]*eta2[itype][jtype]/r_m[itype][jtype])*pow(a1,eta2[itype][jtype]-1)*a2;
				fpair = factor_soft * forcesoft;
				f[i][0] += delx * fpair;
				f[i][1] += dely * fpair;
				f[i][2] += delz * fpair;
				if (newton_pair || j < nlocal){
					f[j][0] -= delx * fpair;
					f[j][1] -= dely * fpair;
					f[j][2] -= delz * fpair;
				}
				if(eflag) {
					evdwl=factor_soft*((v0[itype][jtype])*pow(a1,eta2[itype][jtype])-offset[itype][jtype]);
				}
				if(evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
				
				
				}
			}
			}
		
		if (vflag_fdotr) virial_fdotr_compute();
		}
/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairSoftCut::write_restart(FILE *fp)
{
	write_restart_settings(fp);
	int i, j;
	for(i=1; i<= atom->ntypes; i++)
		for(j=i; j<= atom->ntypes; j++){
			fwrite(&setflag[i][j], sizeof(int), 1, fp);
			if(setflag[i][j]){
				fwrite(&epsilon[i][j], sizeof(double), 1, fp);
				fwrite(&v0[i][j], sizeof(double), 1, fp);
				fwrite(&r_m[i][j], sizeof(double), 1, fp);
				fwrite(&eta1[i][j], sizeof(double), 1, fp);
				fwrite(&eta2[i][j], sizeof(double), 1, fp);
				fwrite(&mu[i][j], sizeof(double), 1, fp);
				fwrite(&nu[i][j], sizeof(double), 1, fp);
				fwrite(&cut[i][j], sizeof(double), 1, fp);
			}
		}
}
	
/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairSoftCut::read_restart(FILE *fp)
{
	read_restart_settings(fp);
	allocate();
	
	int i, j;
	int me = comm->me;
	for(i = 1; i<= atom->ntypes; i++){
		for(j=i; j<=atom->ntypes; j++){
			if(me==0)utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
			MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
			if(setflag[i][j]){
				if(me==0){
					utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &v0[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &r_m[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &eta1[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &eta2[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &mu[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &nu[i][j], sizeof(double), 1, fp, nullptr, error);
					utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
				}
				MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&v0[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&r_m[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&eta1[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&eta2[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&mu[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&nu[i][j], 1, MPI_DOUBLE, 0, world);
				MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
			}
		}
	}
}

/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairSoftCut::write_restart_settings(FILE *fp)
{
  fwrite(&cut_global, sizeof(double), 1, fp);
  fwrite(&offset_flag, sizeof(int), 1, fp);
  fwrite(&mix_flag, sizeof(int), 1, fp);
}

/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairSoftCut::read_restart_settings(FILE *fp)
{
  if (comm->me == 0) {
    utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
    utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
    utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
  }
  MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
  MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
  MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
}

/* ----------------------------------------------------------------------
   proc 0 writes to data file
------------------------------------------------------------------------- */

void PairSoftCut::write_data(FILE *fp)
{
	for (int i = 1; i <= atom->ntypes; i++)
    fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, epsilon[i][i], v0[i][i], r_m[i][i], eta1[i][i], eta2[i][i], mu[i][i], nu[i][i]);
}

/* ----------------------------------------------------------------------
   proc 0 writes all pairs to data file
------------------------------------------------------------------------- */

void PairSoftCut::write_data_all(FILE *fp)
{
  for (int i = 1; i <= atom->ntypes; i++)
    for (int j = i; j <= atom->ntypes; j++)
      fprintf(fp, "%d %d %g %g %g %g %g %g %g %g\n", i, j, epsilon[i][j], v0[i][j], eta1[i][j],
              eta2[i][j], mu[i][j], nu[i][j], r_m[i][j], cut[i][j]);
}

void *PairSoftCut::extract(const char *str, int &dim)
{
	dim =2;
	if(strcmp(str, "epsilon")==0) return(void*) epsilon;
	return nullptr;
	
}
double PairSoftCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, double /*factor_coul*/, double factor_soft, double &fforce)
{
	double dr ,r, a1, a2;
	r = sqrt(rsq);
	dr = r/r_m[itype][jtype];
	a1 = 1 - pow(dr,eta1[itype][jtype]);
	a2 = pow(dr, eta1[itype][jtype]-1);
	
	fforce = factor_soft * (v0[itype][jtype]*eta1[itype][jtype]*eta2[itype][jtype]/r_m[itype][jtype])*pow(a1,eta2[itype][jtype]-1)*a2;
	
	return factor_soft * ((v0[itype][jtype]*pow(a1,eta2[itype][jtype])-offset[itype][jtype]));


}