/* ----------------------------------------------------------------------
   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: Andrew Jewett  (jewett.aij at gmail)
------------------------------------------------------------------------- */

#include <cstdlib>
#include <cstring>
#include <cassert>
using namespace std;
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair_lj_softcore2.h"
#define NDEBUG //(disable assert())
#define ENABLE_SOFTCORE2

using namespace LAMMPS_NS;


PairLJSoftcore2::
PairLJSoftcore2(LAMMPS *lmp) : PairLJRemix(lmp)
{
  cut_core_sigma = -1.0;     // <--This dissables the softcore cutoff (default)
  ignore_large_rcore = false;  // Generate error messages if cutoff is too large.
}


PairLJSoftcore2::
~PairLJSoftcore2()
{
  #ifdef ENABLE_SOFTCORE2
  if (allocated) {
    memory->destroy(cut_core_sqd);
    memory->destroy(core_coeff);
    memory->destroy(core_D_coeff);
    memory->destroy(core_offset);
    //fprintf(stderr, "~~~~ invoking ~PairLJSoftcore2 check successful\n");
  }
  #endif
}

void 
PairLJSoftcore2::
allocate() {
  #ifdef ENABLE_SOFTCORE2
  int n = atom->ntypes;
  memory->create(cut_core_sqd,n+1,n+1,"pair:cut_core_sqd");
  memory->create(core_coeff,n+1,n+1,"pair:core_coeff");
  memory->create(core_D_coeff,n+1,n+1,"pair:core_D_coeff");
  memory->create(core_offset,n+1,n+1,"pair:core_offset");
  #endif
}


void PairLJSoftcore2::
settings(int &narg, char **arg)
{
  int i=0;
  while (i < narg)
  {
    int nargs_delete = 0;
    // Is this argument a number?
    char *pc;
    strtod(arg[i], &pc);
    if (strlen(pc) != 0) // then arg[i] is not a number
    { 
      if ((strcmp(arg[i], "rsoftcore") == 0) ||
	  (strcmp(arg[i], "yesshutup") == 0)) {

        char err_msg[1024];
        #ifdef ENABLE_SOFTCORE2
        if (i+1 >= narg) {
          sprintf(err_msg,"The \"%s\" keyword must be followed by a number.\n"
                  "(the ratio of the core cutoff radius divided by sigma)\n",
                  arg[i]);
          error->all(FLERR, err_msg);
        }
        ignore_large_rcore = false;
        if (strcmp(arg[i], "yesshutup") == 0)
          ignore_large_rcore = true;

        cut_core_sigma = force->numeric(arg[i+1]);
	nargs_delete = 2;
        #else
        sprintf(err_msg, "This version of LAMMPS was not compiled with support for \"rcore\".\n"
                "To enable softcore potentials, edit the\n"
                "   \"%s\" file, and\n"
                
                "uncomment the line containing \"#define ENABLE_SOFTCORE2\"\n"
                "(near the beginning of the file) and recompile\n", __FILE__);
        error->all(FLERR, err_msg);
        #endif //#ifdef ENABLE_SOFTCORE2
      } //if ((strcmp(arg[i], "rcore") == 0) || ...
      else
        i++;

      // now delete the argument(s) we just processed
      for (int j = i; j+nargs_delete < narg; j++)
        arg[j] = arg[j+nargs_delete];
      narg -= nargs_delete;
    } //if (strlen(pc) != 0)
    else
      i++;
  } //while (i < narg)

  #ifdef ENABLE_SOFTCORE2
  if ((cut_core_sigma > 0.0) && (coeff_style == COEFF_AB))
    error->all(FLERR, "pair_style arguments clash.  You can not specify a core cutoff distance\n"
               "if you are using the \"ab\" style.  (The core cutoff is in units of sigma.)\n");
  #endif
}


void PairLJSoftcore2::
init_one(int i, int j,
                  double **aa_sigma,
                  double **aa_cut_core_sqd,
                  double **aa_core_coeff,
                  double **aa_core_D_coeff,
                  double **aa_core_offset,
		  double **aa_lj3,
		  double **aa_lj4)
{
  #ifdef ENABLE_SOFTCORE2
  if ((cut_core_sigma > 0.0) && (coeff_style != COEFF_AB)) {
    // rc = cut_core_sigma * sigma
    //
    // The potential and derivative should be continuous at rc
    //                 core_coeff                       lj3            lj4
    //    U(rc)    =   ---------- + core_offset   =   -------     -   ------
    //                    rc^2                          rc^a          rc^b
    //
    // - dU(r)|       2*core_coeff                     a * lj3        b * lj4
    //   -----|    =   ------------               =    --------   -   -------
    //    dr  |rc         rc^3                         rc^(a+1)       rc^(b+1)
    //
    //
    // -->  core_coeff  =  (a/2) * (lj3/rc^(a-2))  -  (b/2) * (lj4/rc^(b-2))
    // -->  core_offset =  lj3/rc^a - lj4/rc^b - core_coeff/rc^2
    double rc;
    bool core_is_large = false;

    rc                       = cut_core_sigma * fabs(aa_sigma[i][j]);
    aa_cut_core_sqd[i][j]    = rc * rc;
    aa_core_coeff[i][j]=((exponent_a/2.0)*aa_lj3[i][j]/pow(rc,(exponent_a-2.0)))
                    + ((-exponent_b/2.0)*aa_lj4[i][j]/pow(rc,(exponent_b-2.0)));
    aa_core_offset[i][j]     = (aa_lj3[i][j]/pow(rc,exponent_a))
                             + (-aa_lj4[i][j]/pow(rc,exponent_b)) 
                             + (-aa_core_coeff[i][j]/(rc*rc));
    aa_core_D_coeff[i][j]    = 2 * aa_core_coeff[i][j];    
    if (aa_core_coeff[i][j] < 0.0) 
      core_is_large = true;

    //   Error checks:

    if (core_is_large && (! ignore_large_rcore))
    {
        char err_msg[1024];
        sprintf(err_msg, "You have an attractive singularity at r=0 (between atom types %d & %d).\n"
                "This is probably because the cutoff radius for your soft core repulsion\n"
                "(\"rcore\" argument) is too large.  (Yours is %g x sigma)\n"
                "Remember that the number after the \"rcore\" keyword should be in units\n"
                "of \"sigma\", not in natural physical units (like Angstroms).\ns"
                "\n"
                "(If you know what you are doing and wish to override this warming\n"
                "then replace \"rcore\" with \"yesshutup\" before the numeric argument.)\n", i, j, cut_core_sigma);
        error->all(FLERR, err_msg);
    }

    // commenting out the next check:
    //   I can't think of any reason the core cutoff should
    //   ever be allowed to exceed any of the other lennard-jones cutoffs.
    //   When this occured, it used to cause an error, but I commented it out.
    //   Either way, the code scales the LJ contribution by the switching
    //   function. I might change my mind later and uncomment this error below:
    //   -Andrew 2012-5-23
    //
    //if (cut_core_sqd[i][j] > cut_lj_innersq) {
    //  char err_msg[2048];
    //  sprintf(err_msg, "Core cuttoff for atom types %d & %d exceeds the inner LJ switch cutoff radius\n", i,j);
    //  error->all(FLERR, err_msg);
    //}

  } //if ((cut_core_sigma > 0.0) && (coeff_style != COEFF_AB))

  else {
    aa_cut_core_sqd[i][j] = -1.0;
  }

  if (cut_core_sigma > 0.0) {
    // i <--> j for params for forces at distances below the "core" cutoff
    aa_cut_core_sqd[j][i]  = aa_cut_core_sqd[i][j];
    aa_core_coeff[j][i]    = aa_core_coeff[i][j];
    aa_core_D_coeff[j][i]  = aa_core_D_coeff[i][j];
    aa_core_offset[j][i]   = aa_core_offset[i][j];
  }
  #endif //#ifdef ENABLE_SOFTCORE2
} //PairLJSoftcore2::init_one(i)


void PairLJSoftcore2::
write_restart_settings(FILE *fp)
{
  #ifdef ENABLE_SOFTCORE2
  fwrite(&cut_core_sigma,sizeof(double),1,fp);
  #endif
}

void PairLJSoftcore2::
read_restart_settings(FILE *fp)
{
  #ifdef ENABLE_SOFTCORE2
  int me = comm->me;
  if (me == 0) {
    fread(&cut_core_sigma,sizeof(double),1,fp);
  }
  MPI_Bcast(&cut_core_sigma,1,MPI_DOUBLE,0,world);
  #endif
}



void PairLJSoftcore2::
write_restart_pair(FILE *fp,
		   int i, int j,
		   double **aa_cut_core_sqd,
		   double **aa_core_coeff,
		   double **aa_core_D_coeff,
		   double **aa_core_offset)
{
  #ifdef ENABLE_SOFTCORE2
  fwrite(&aa_cut_core_sqd[i][j], sizeof(double),1,fp);
  fwrite(&aa_core_coeff[i][j], sizeof(double),1,fp);
  fwrite(&aa_core_D_coeff[i][j], sizeof(double),1,fp);
  fwrite(&aa_core_offset[i][j], sizeof(double),1,fp);
  #endif
}


void PairLJSoftcore2::
read_restart_pair(FILE *fp,
		  int i, int j,
		  double **aa_cut_core_sqd,
		  double **aa_core_coeff,
		  double **aa_core_D_coeff,
		  double **aa_core_offset)
{
  #ifdef ENABLE_SOFTCORE2
  int me = comm->me;
  if (me == 0) {
    fread(&aa_cut_core_sqd[i][j],sizeof(double),1,fp);
    fread(&aa_core_coeff[i][j],sizeof(double),1,fp);
    fread(&aa_core_D_coeff[i][j],sizeof(double),1,fp);
    fread(&aa_core_offset[i][j],sizeof(double),1,fp);
  } //if (me == 0)

  MPI_Bcast(&aa_cut_core_sqd[i][j],1,MPI_DOUBLE,0,world);
  MPI_Bcast(&aa_core_coeff[i][j],1,MPI_DOUBLE,0,world);
  MPI_Bcast(&aa_core_D_coeff[i][j],1,MPI_DOUBLE,0,world);
  MPI_Bcast(&aa_core_offset[i][j],1,MPI_DOUBLE,0,world);
  #endif //#ifdef ENABLE_SOFTCORE2
}