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

#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_charmm_coul_charmm_inter.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include <cassert>
using namespace LAMMPS_NS;

#define ENABLE_INTER_DIF
#define ENABLE_SOFTCORE2
#define NDEBUG //(disable assert())

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

PairLJCharmmCoulCharmmInter::
PairLJCharmmCoulCharmmInter(LAMMPS *lmp) : PairLJSoftcore2(lmp)
{
  implicit = 0;

  // default styles:
  mix_flag = ARITHMETIC;

  cut_coul_and_lj_differ = 0;      //(same switching func used for LJ and COUL)
}


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

PairLJCharmmCoulCharmmInter::~PairLJCharmmCoulCharmmInter()
{
  if (allocated) {
    memory->destroy(cutsq);

    // Tables needed for initialization
    memory->destroy(setflag);

    memory->destroy(eps14);
    memory->destroy(sigma14);
    memory->destroy(kappa14);
    memory->destroy(lambda14);

    #ifdef ENABLE_INTER_DIF
    memory->destroy(epsDIF);
    memory->destroy(sigmaDIF);
    memory->destroy(kappaDIF);
    memory->destroy(lambdaDIF);
    #endif

    // Usually force computation only requires these tables:
    memory->destroy(lj14_1);
    memory->destroy(lj14_2);
    memory->destroy(lj14_3);
    memory->destroy(lj14_4);

    #ifdef ENABLE_INTER_DIF
    memory->destroy(ljDIF_1);
    memory->destroy(ljDIF_2);
    memory->destroy(ljDIF_3);
    memory->destroy(ljDIF_4);
    #endif


    #ifdef ENABLE_SOFTCORE2
    // Arrays describing parameters for the potential used below 
    // the "core" cutoff (soft repulsive function for short distances)
    #ifdef ENABLE_INTER_DIF
    memory->destroy(cut_coreDIF_sqd);
    memory->destroy(core_coeffDIF);
    memory->destroy(core_D_coeffDIF);
    memory->destroy(core_offsetDIF);
    #endif
    // There's no way to modify 14 behavior without changing dihedral_charmm.cpp
    // Disabling the cutoff behavior for 1-4 interactions, and commenting out:
    //memory->destroy(cut_core14_sqd);
    //memory->destroy(core_coeff14);
    //memory->destroy(core_D_coeff14);
    //memory->destroy(core_offset14);
    #endif //#ifdef ENABLE_SOFTCORE2


  } //if (allocated)
} //~PairLJCharmmCoulCharmmInter()


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

void PairLJCharmmCoulCharmmInter::compute(int eflag, int vflag)
{
  int i,j,ii,jj,inum,jnum,itype,jtype;
  double q_i,q_j,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
  double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
  double   cutljsq_minus_rsq,   cutljsq_minus_rsq_sq;
  double cutcoulsq_minus_rsq, cutcoulsq_minus_rsq_sq;
  double philj, switch1_coul, switch2_coul, switch1_lj, switch2_lj;
  int *ilist,*jlist,*numneigh,**firstneigh;

  evdwl = ecoul = 0.0;
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

  double **x = atom->x;
  double **f = atom->f;
  double *q = atom->q;
  int *type = atom->type;
  int nlocal = atom->nlocal;
  double *special_coul = force->special_coul;
  double *special_lj = force->special_lj;
  int newton_pair = force->newton_pair;
  double qqrd2e = force->qqrd2e;

  int *molecule = atom->molecule;
  double _lj1, _lj2, _lj3, _lj4;
  #ifdef ENABLE_SOFTCORE2
  double _cut_core_sqd, _core_coeff, _core_D_coeff, _core_offset;
  #endif

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

  // loop over neighbors of my atoms

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    q_i = q[i];
    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_lj = special_lj[sbmask(j)];
      factor_coul = special_coul[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;

      if (rsq < cut_bothsq) {
        r2inv = 1.0/rsq;

        if (rsq < cut_coulsq) {
          //Only evaluate the expensive sqrt() operation for nonzero charges.
          //On intel i7 hardware, this if statement accelerates simulations
          //noticably for uncharged systems (approx %30).  On other hardware 
          //it might not, but the extra check should not slow things down much.
          q_j = q[j];
          if ((q_i != 0.0) && (q_j != 0.0))
            forcecoul = qqrd2e * q_i * q_j * sqrt(r2inv);
          else
            forcecoul = 0.0;

          if (rsq > cut_coul_innersq) {
            cutcoulsq_minus_rsq = (cut_coulsq-rsq);
            cutcoulsq_minus_rsq_sq = cutcoulsq_minus_rsq * cutcoulsq_minus_rsq;
            switch1_coul = cutcoulsq_minus_rsq_sq *
              (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) * denom_coul_inv;
            switch2_coul = 12.0*rsq * cutcoulsq_minus_rsq * 
              (rsq-cut_coul_innersq) * denom_coul_inv;
            forcecoul *= switch1_coul + switch2_coul;
          }
        } else forcecoul = 0.0;

        if (rsq < cut_ljsq) {
          r6inv = r2inv*r2inv*r2inv;
          jtype = type[j];

          #ifdef ENABLE_INTER_DIF
          if (molecule[i] == molecule[j]) {
          #endif
            _lj1 = lj1[itype][jtype];
            _lj2 = lj2[itype][jtype];
            _lj3 = lj3[itype][jtype];
            _lj4 = lj4[itype][jtype];
            #ifdef ENABLE_SOFTCORE2
            _cut_core_sqd = cut_core_sqd[itype][jtype];
            _core_coeff   =   core_coeff[itype][jtype];
            _core_D_coeff = core_D_coeff[itype][jtype];
            _core_offset  =  core_offset[itype][jtype];
            #endif
          #ifdef ENABLE_INTER_DIF
          }
          else {
            _lj1 = ljDIF_1[itype][jtype];
            _lj2 = ljDIF_2[itype][jtype];
            _lj3 = ljDIF_3[itype][jtype];
            _lj4 = ljDIF_4[itype][jtype];
            #ifdef ENABLE_SOFTCORE2
            _cut_core_sqd = cut_coreDIF_sqd[itype][jtype];
            _core_coeff   =   core_coeffDIF[itype][jtype];
            _core_D_coeff = core_D_coeffDIF[itype][jtype];
            _core_offset  =  core_offsetDIF[itype][jtype];
            #endif
          } //else clause for "if (molecule[i] == molecule[j])"
          #endif //#ifdef ENABLE_INTER_DIF

          #ifdef ENABLE_SOFTCORE2
          if (rsq > _cut_core_sqd) {
          #endif

            //  Standard Lennard-Jones function
            //    Ulj(r) = (_lj3/r^12- _lj4/r^6) * switch1_lj(r)
            forcelj = r6inv * (_lj1*r6inv - _lj2);

            if (rsq > cut_lj_innersq) {
              switch1_lj = switch1_coul;
              switch2_lj = switch2_coul;
              if (cut_coul_and_lj_differ) { // then recompute switch1 & switch2
                cutljsq_minus_rsq = (cut_ljsq-rsq);
                cutljsq_minus_rsq_sq = cutljsq_minus_rsq * cutljsq_minus_rsq;
                switch1_lj = cutljsq_minus_rsq_sq *
                  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
                switch2_lj = 12.0*rsq * cutljsq_minus_rsq * 
                  (rsq-cut_lj_innersq) * denom_lj_inv;
              }
              philj = r6inv * (_lj3*r6inv - _lj4);
              forcelj = forcelj*switch1_lj + philj*switch2_lj;
            }

          #ifdef ENABLE_SOFTCORE2
          }
          else {
            //  Soft-core repulsive function
            //    Ulj(r) = (_core_coeff/r^2  + _core_offset) * switch1_lj(r)
            forcelj = r2inv * _core_D_coeff;

            if (rsq > cut_lj_innersq) {
              // The next 8 lines should be identical to the lines above.  (No
              // way to move them out without creating an extra if-then check.)
              switch1_lj = switch1_coul;
              switch2_lj = switch2_coul;
              if (cut_coul_and_lj_differ) { // then recompute switch1 & switch2
                cutljsq_minus_rsq = (cut_ljsq-rsq);
                cutljsq_minus_rsq_sq = cutljsq_minus_rsq * cutljsq_minus_rsq;
                switch1_lj = cutljsq_minus_rsq_sq *
                  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
                switch2_lj = 12.0*rsq * cutljsq_minus_rsq * 
                  (rsq-cut_lj_innersq) * denom_lj_inv;
              }
              philj = r2inv * _core_coeff + _core_offset;
              //forcelj = forcelj*switch1_lj + philj*switch2_lj;
              //fprintf(stderr,"ACCESS ERROR: reached %s:%d\n",__FILE__,__LINE__);
              //exit(-1);
            }
          } //else clause for "if (rsq > _cut_core_sqd)"
          #endif //#ifdef ENABLE_SOFTCORE2

        } else forcelj = 0.0;

        fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
        //fprintf(stderr, "~~~~ fpair=%g\n", fpair);

        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) {
          if (rsq < cut_coulsq) {
            // Recompute the energy (philj) just in case it was
            // not computed above, and set evdwl equal to it.
            q_j = q[j];
            if ((q_i != 0.0) && (q_j != 0.0))
              ecoul = qqrd2e * q_i * q_j * sqrt(r2inv);
            else
              ecoul = 0.0;
            if (rsq > cut_coul_innersq)
              ecoul *= switch1_coul; //switch1_coul is calculated earlier
                                     //whenever rsq > cut_coul_innersq

            ecoul *= factor_coul;
          } else ecoul = 0.0;
          //fprintf(stderr, "~~~~ switch1_coul=%g, ecoul=%g\n", switch1_coul, ecoul);
          if (rsq < cut_ljsq) {

            #ifdef ENABLE_SOFTCORE2
            if (rsq > _cut_core_sqd)
            #endif

              philj = r6inv*(_lj3*r6inv - _lj4);

            #ifdef ENABLE_SOFTCORE2
            else {
              philj = r2inv * _core_coeff + _core_offset;
              //fprintf(stderr,"ACCESS ERROR: reached %s:%d\n",__FILE__,__LINE__);
              //exit(-1);
            }
            #endif

            if (rsq > cut_lj_innersq) {
              philj *= switch1_lj; //switch1_lj is calculated earlier
                                   //whenever rsq > cut_lj_innersq
            }
          } else philj = 0.0;

          philj *= factor_lj;
          evdwl = philj;
          //fprintf(stderr, "~~~~ switch1_lj=%g, evdwl=%g\n", switch1_lj,evdwl);
        }

        if (evflag) ev_tally(i,j,nlocal,newton_pair,
                             evdwl,ecoul,fpair,delx,dely,delz);

      } //if (rsq < cut_bothsq)
    } //for (jj = 0; jj < jnum; jj++)
  } //for (ii = 0; ii < jnum; ii++)

  if (vflag_fdotr) virial_fdotr_compute();
} //void PairLJCharmmCoulCharmmInter::compute()


/* ----------------------------------------------------------------------
   allocate all arrays
------------------------------------------------------------------------- */

void PairLJCharmmCoulCharmmInter::allocate()
{
  allocated = 1;
  int n = atom->ntypes;

  memory->create(cutsq,n+1,n+1,"pair:cutsq");

  // Tables needed for initialization
  memory->create(setflag,n+1,n+1,"pair:setflag");
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      setflag[i][j] = 0;

  PairLJRemix::allocate();     //this allocates the epsilon, sigma, 
                               //kappa, lambda, lj1, lj2, lj3, lj4 arrays

  PairLJSoftcore2::allocate(); //this allocates the cut_core_sqd, core_coeff, 
                               //core_D_coeff, core_offset arrays

  memory->create(eps14,n+1,n+1,"pair:eps14");
  memory->create(sigma14,n+1,n+1,"pair:sigma14");
  memory->create(kappa14,n+1,n+1,"pair:kappa14");
  memory->create(lambda14,n+1,n+1,"pair:lambda14");
  #ifdef DBG_LJ_REMIX_MEM
  for (int itype=0; itype<=n; itype++) {
    for (int jtype=0; jtype<=n; jtype++) {
      eps14[itype][jtype] = 0.0;
      sigma14[itype][jtype] = 0.0;
      kappa14[itype][jtype] = 0.0;
      lambda14[itype][jtype] = 0.0;
    }
  }
  #endif //#ifdef DBG_LJ_REMIX_MEM

  #ifdef ENABLE_INTER_DIF
  memory->create(epsDIF,n+1,n+1,"pair:epsDIF");
  memory->create(sigmaDIF,n+1,n+1,"pair:sigmaDIF");
  memory->create(kappaDIF,n+1,n+1,"pair:kappaDIF");
  memory->create(lambdaDIF,n+1,n+1,"pair:lambdaDIF");
  #ifdef DBG_LJ_REMIX_MEM
  for (int itype=0; itype<=n; itype++) {
    for (int jtype=0; jtype<=n; jtype++) {
      epsDIF[itype][jtype] = 0.0;
      sigmaDIF[itype][jtype] = 0.0;
      kappaDIF[itype][jtype] = 0.0;
      lambdaDIF[itype][jtype] = 0.0;
    }
  }
  #endif //#ifdef DBG_LJ_REMIX_MEM
  #endif //#ifdef ENABLE_INTER_DIF

  memory->create(lj14_1,n+1,n+1,"pair:lj14_1");
  memory->create(lj14_2,n+1,n+1,"pair:lj14_2");
  memory->create(lj14_3,n+1,n+1,"pair:lj14_3");
  memory->create(lj14_4,n+1,n+1,"pair:lj14_4");
  #ifdef DBG_LJ_REMIX_MEM
  for (int itype=0; itype<=n; itype++) {
    for (int jtype=0; jtype<=n; jtype++) {
      lj14_1[itype][jtype] = 0.0;
      lj14_2[itype][jtype] = 0.0;
      lj14_3[itype][jtype] = 0.0;
      lj14_4[itype][jtype] = 0.0;
    }
  }
  #endif //#ifdef DBG_LJ_REMIX_MEM

  #ifdef ENABLE_INTER_DIF
  memory->create(ljDIF_1,n+1,n+1,"pair:ljDIF_1");
  memory->create(ljDIF_2,n+1,n+1,"pair:ljDIF_2");
  memory->create(ljDIF_3,n+1,n+1,"pair:ljDIF_3");
  memory->create(ljDIF_4,n+1,n+1,"pair:ljDIF_4");
  #ifdef DBG_LJ_REMIX_MEM
  for (int itype=0; itype<=n; itype++) {
    for (int jtype=0; jtype<=n; jtype++) {
      ljDIF_1[itype][jtype] = 0.0;
      ljDIF_2[itype][jtype] = 0.0;
      ljDIF_3[itype][jtype] = 0.0;
      ljDIF_4[itype][jtype] = 0.0;
    }
  }
  #endif //#ifdef DBG_LJ_REMIX_MEM
  #endif //#ifdef ENABLE_INTER_DIF

  #ifdef ENABLE_SOFTCORE2
  // Parameters needed below the softcore cutoff radius:
  #ifdef ENABLE_INTER_DIF
  memory->create(cut_coreDIF_sqd,n+1,n+1,"pair:cut_coreDIF_sqd");
  memory->create(core_coeffDIF,n+1,n+1,"pair:core_coeffDIF");
  memory->create(core_D_coeffDIF,n+1,n+1,"pair:core_D_coeffDIF");
  memory->create(core_offsetDIF,n+1,n+1,"pair:core_offsetDIF");
  #ifdef DBG_LJ_REMIX_MEM
  for (int itype=0; itype<=n; itype++) {
    for (int jtype=0; jtype<=n; jtype++) {
      cut_coreDIF_sqd[itype][jtype] = 0.0;
      core_coeffDIF[itype][jtype]   = 0.0;
      core_D_coeffDIF[itype][jtype] = 0.0;
      core_offsetDIF[itype][jtype]  = 0.0;
    }
  }
  #endif //#ifdef DBG_LJ_REMIX_MEM
  #endif //#ifdef ENABLE_INTER_DIF

  // There's no way to modify 14 behavior without changing dihedral_charmm.cpp
  // Disabling the cutoff behavior for 1-4 interactions, and commenting out:
  //memory->create(cut_core14_sqd,n+1,n+1,"pair:cut_core14_sqd");
  //memory->create(core_coeff14,n+1,n+1,"pair:core_coeff14");
  //memory->create(core_D_coeff14,n+1,n+1,"pair:core_D_coeff14");
  //memory->create(core_offset14,n+1,n+1,"pair:core_offset14");
  #endif //#ifdef ENABLE_SOFTCORE2

} //PairLJCharmmCoulCharmmInter::allocate()


/* ----------------------------------------------------------------------
   global settings
   unlike other pair styles,
     there are no individual pair settings that these override
------------------------------------------------------------------------- */

void PairLJCharmmCoulCharmmInter::settings(int narg, char **arg)
{
  PairLJRemix::settings(narg, arg);
  PairLJSoftcore2::settings(narg, arg);

  // All remaining arguments should be numeric

  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
    { 
      error->all(FLERR,"Unrecognized argument to pair_style command");
    }
    i++;
  }

  // What's left over are the numerical arguments:
  if ((narg != 2) && (narg != 4))
    error->all(FLERR, "Illegal pair_style command (wrong number of arguments)");

  cut_lj_inner = force->numeric(FLERR,arg[0]);
  cut_lj = force->numeric(FLERR,arg[1]);
  cut_coul_and_lj_differ = 0;
  cut_coul_inner = cut_lj_inner;
  cut_coul = cut_lj;
  if (narg == 4) {
    cut_coul_inner = force->numeric(FLERR,arg[2]);
    cut_coul = force->numeric(FLERR,arg[3]);
    // Set cut_coul_and_lj_differ to 1 if different cutoffs for coul and lj.
    // Testing numeric equivalence is dangerous. Safer to test the strings
    if ((! strcmp(arg[0],arg[2])) || (! strcmp(arg[1],arg[3])))
      cut_coul_and_lj_differ = 1;
  }

} //void PairLJCharmmCoulCharmmInter::settings()




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

void PairLJCharmmCoulCharmmInter::coeff(int narg, char **arg)
{
  char err_msg[1024];
  if (coeff_style == COEFF_AB) {
    if ((narg != 2+2) && (narg != 4+2) && (narg != 6+2))
      error->all(FLERR,"Incorrect args for pair coefficients");
  }
  else if (! ((coeff_style == COEFF_EPS_SIG_KL) ||
              (coeff_style == COEFF_EPS_SIG_4K4L) ||
              (coeff_style == COEFF_EPS_SIG_1K2L)))
  {
    if ((narg != 2+2) && (narg != 4+2) && (narg != 6+2))
      error->all(FLERR,"Incorrect args for pair coefficients");
  }
  else
  {
    if ((narg != 4+2) && (narg != 8+2) && (narg != 12+2))
      error->all(FLERR,"Incorrect args for pair coefficients");
  }

  if (!allocated) allocate();

  int ilo,ihi,jlo,jhi;
  force->bounds(arg[0],atom->ntypes,ilo,ihi);
  force->bounds(arg[1],atom->ntypes,jlo,jhi);

  double A_one = -1.0;
  double B_one = -1.0;
  double kappa_one = 1.0;
  double lambda_one = -1.0;
  double epsilon_one = force->numeric(FLERR,arg[2]);
  double sigma_one = force->numeric(FLERR,arg[3]);

  double A14_one = -1.0;
  double B14_one = -1.0;
  double kappa14_one = 1.0;
  double lambda14_one = -1.0;
  double eps14_one = epsilon_one;
  double sigma14_one = sigma_one;

  #ifdef ENABLE_INTER_DIF
  double ADIF_one = -1.0;
  double BDIF_one = -1.0;
  double kappaDIF_one = 1.0;
  double lambdaDIF_one = -1.0;
  double epsDIF_one = epsilon_one;
  double sigmaDIF_one = sigma_one;
  #endif

  if (coeff_style == COEFF_AB) {
    A_one = force->numeric(FLERR,arg[2]);
    B_one = force->numeric(FLERR,arg[3]);
    A14_one = A_one;
    B14_one = B_one;
    #ifdef ENABLE_INTER_DIF
    ADIF_one = A_one;
    BDIF_one = B_one;
    #endif
    if (narg >= 4+2) {
      A14_one = force->numeric(FLERR,arg[2]);
      B14_one = force->numeric(FLERR,arg[3]);
      if (narg == 6+2) {
        #ifdef ENABLE_INTER_DIF
        ADIF_one = force->numeric(FLERR,arg[6]);
        BDIF_one = force->numeric(FLERR,arg[7]);
        #else
        sprintf(err_msg, "This version of LAMMPS was compiled without support for\n"
                "custom inter-vs-intra-molecular pairwise forces.  To fix this,\n"
                "the \"ENABLE_INTER_DIF\" macro must be defined in source\n"
                "file \"%s\".\n", __FILE__);
        error->all(FLERR, err_msg);
        #endif //#ifdef ENABLE_INTER_DIF
      }
    }
  }
  else if (! ((coeff_style == COEFF_EPS_SIG_KL) ||
              (coeff_style == COEFF_EPS_SIG_4K4L) ||
              (coeff_style == COEFF_EPS_SIG_1K2L)))
  {
    if (narg >= 4+2) {
      eps14_one = force->numeric(FLERR,arg[4]);
      sigma14_one = force->numeric(FLERR,arg[5]);
      if (narg == 6+2) {
        #ifdef ENABLE_INTER_DIF
        epsDIF_one = force->numeric(FLERR,arg[6]);
        sigmaDIF_one = force->numeric(FLERR,arg[7]);
        #else
        sprintf(err_msg, "This version of LAMMPS was compiled without support for\n"
                "custom inter-vs-intra-molecular pairwise forces.  To fix this,\n"
                "the \"ENABLE_INTER_DIF\" macro must be defined in source\n"
                "file \"%s\".\n", __FILE__);
        error->all(FLERR, err_msg);
        #endif //#ifdef ENABLE_INTER_DIF
      }
    }
  }
  else {
    kappa_one = force->numeric(FLERR,arg[4]);
    lambda_one = force->numeric(FLERR,arg[5]);

    kappa14_one = kappa_one;
    lambda14_one = lambda_one;

    #ifdef ENABLE_INTER_DIF
    kappaDIF_one = kappa_one;
    lambdaDIF_one = lambda_one;
    #endif

    if (narg >= 8+2) {
      eps14_one = force->numeric(FLERR,arg[6]);
      sigma14_one = force->numeric(FLERR,arg[7]);
      kappa14_one = force->numeric(FLERR,arg[8]);
      lambda14_one = force->numeric(FLERR,arg[9]);
      if (narg == 12+2) {
        #ifdef ENABLE_INTER_DIF
        epsDIF_one = force->numeric(FLERR,arg[10]);
        sigmaDIF_one = force->numeric(FLERR,arg[11]);
        kappaDIF_one = force->numeric(FLERR,arg[12]);
        lambdaDIF_one = force->numeric(FLERR,arg[13]);
        #else
        sprintf(err_msg, "This version of LAMMPS was compiled without support for\n"
                "custom inter-vs-intra-molecular pairwise forces.  To fix this,\n"
                "the \"ENABLE_INTER_DIF\" macro must be defined in source\n"
                "file \"%s\".\n", __FILE__);
        error->all(FLERR, err_msg);
        #endif //#ifdef ENABLE_INTER_DIF
      }
    }
  }

  int count = 0;
  for (int i = ilo; i <= ihi; i++) {
    for (int j = MAX(jlo,i); j <= jhi; j++) {
      if (coeff_style == COEFF_AB) {
        lj3[i][j] = A_one;
        lj4[i][j] = B_one;
        lj14_3[i][j] = A14_one;
        lj14_4[i][j] = B14_one;
        #ifdef ENABLE_INTER_DIF
        ljDIF_3[i][j] = ADIF_one;
        ljDIF_4[i][j] = BDIF_one;
        #endif
      }
      else{
        epsilon[i][j] = epsilon_one;
        sigma[i][j] = sigma_one;
        kappa[i][j] = kappa_one;
        lambda[i][j] = lambda_one;
        eps14[i][j] = eps14_one;
        sigma14[i][j] = sigma14_one;
        kappa14[i][j] = kappa14_one;
        lambda14[i][j] = lambda14_one;
        #ifdef ENABLE_INTER_DIF
        epsDIF[i][j] = epsDIF_one;
        sigmaDIF[i][j] = sigmaDIF_one;
        kappaDIF[i][j] = kappaDIF_one;
        lambdaDIF[i][j] = lambdaDIF_one;
        #endif
      }
      setflag[i][j] = 1;
      count++;
    }
  }

  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
} //PairLJCharmmCoulCharmmInter::coeff()


/* ----------------------------------------------------------------------
   init specific to this pair style
------------------------------------------------------------------------- */

void PairLJCharmmCoulCharmmInter::init_style()
{
  if (!atom->q_flag)
    error->all(FLERR,
               "Pair style lj/charmm/coul/charmm/inter requires atom attribute q");

  neighbor->request(this);

  // require cut_lj_inner < cut_lj, cut_coul_inner < cut_coul

  if (((cut_lj_inner >= cut_lj) && (cut_lj_inner > 0.0))
      || 
      ((cut_coul_inner >= cut_coul) && (cut_coul_inner > 0.0)))
  {
    error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
  }

  cut_lj_innersq = cut_lj_inner * cut_lj_inner;
  cut_ljsq = cut_lj * cut_lj;
  cut_coul_innersq = cut_coul_inner * cut_coul_inner;
  cut_coulsq = cut_coul * cut_coul;
  cut_bothsq = MAX(cut_ljsq,cut_coulsq);

  denom_lj_inv  = 1.0 / ((cut_ljsq-cut_lj_innersq)*
                         (cut_ljsq-cut_lj_innersq)*
                         (cut_ljsq-cut_lj_innersq));
  denom_coul_inv = 1.0 / ((cut_coulsq-cut_coul_innersq)*
                          (cut_coulsq-cut_coul_innersq)*
                          (cut_coulsq-cut_coul_innersq));
}




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

double PairLJCharmmCoulCharmmInter::init_one(int i, int j)
{
  // Update the intra-molecular Lennard-Jones parameter arrays:
  PairLJRemix::init_one(i, j,
                        epsilon, sigma, kappa, lambda,
                        lj1, lj2, lj3, lj4);

  // Update the intra-molecular 1-4 Lennard-Jones parameter arrays:
  PairLJRemix::init_one(i, j,
                        eps14, sigma14, kappa14, lambda14,
                        lj14_1, lj14_2, lj14_3, lj14_4);

  #ifdef ENABLE_INTER_DIF
  // Update the inter-molecular Lennard-Jones parameter arrays:
  PairLJRemix::init_one(i, j, 
                        epsDIF, sigmaDIF, kappaDIF, lambdaDIF,
                        ljDIF_1, ljDIF_2, ljDIF_3, ljDIF_4);
  #endif



  #ifdef ENABLE_SOFTCORE2

  PairLJSoftcore2::init_one(i, j, 
                            sigma,
                            cut_core_sqd,
                            core_coeff,
                            core_D_coeff,
                            core_offset,
                            lj3,
                            lj4);

  #ifdef ENABLE_INTER_DIF
  PairLJSoftcore2::init_one(i, j, 
                            sigmaDIF,
                            cut_coreDIF_sqd,
                            core_coeffDIF,
                            core_D_coeffDIF,
                            core_offsetDIF,
                            ljDIF_3,
                            ljDIF_4);
  #endif

  // There's no way to modify 14 behavior without changing dihedral_charmm.cpp
  // Disabling the cutoff behavior for 1-4 interactions, and commenting out:
  //PairLJSoftcore2::init_one(i, j, 
  //                          sigma14,
  //                          cut_core14_sqd,
  //                          core_coeff14,
  //                          core_D_coeff14,
  //                          core_offset14);
  //                          lj14_3,
  //                          lj14_4);


  if ((cut_core_sigma > 0.0) && (coeff_style != COEFF_AB))
  {
    // Here we try to insure that 1-4 interactions are explicitly 
    // turned off by the user whenever "core" cutoffs are in use.
    if ((lj14_3[i][j] > 1.0e-7 * fabs(lj3[i][j])) ||
        (lj14_4[i][j] > 1.0e-7 * fabs(lj4[i][j]))) {
      error->all(FLERR, 
                 "If you plan to use softer 1/r^2 repulsive approximation to the Lennard-Jones\n"
                 "interaction at small distances, then you must explicitly turn off all\n"
                 "1-4 interactions by setting the 1-4 epsilon coefficien to \"0\".\n"
                 "(These parameters can not be ommitted because the default is non-zero.)"
                 "\n"
                 "Reason: The dihedral code does know about this 1/r^2 approximation.)\n"
                 "  IF YOU NEED 1-4 INTERACTIONS, THIS CAN BE FIXED EASILY.\n"
                 "(Email the pair style's author or the lammps-users mailing list.)\n");
    }
  }
  #endif //#ifdef ENABLE_SOFTCORE2

  double cut = MAX(cut_lj, cut_coul);

  return cut;
} //PairLJCharmmCoulCharmmInter::init_one(int, int)



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

void PairLJCharmmCoulCharmmInter::write_restart(FILE *fp)
{
  write_restart_settings(fp);

  //Just write the raw lennard jones AB parameters.  Don't write the epsilon,
  //sigma, and lambda parameters (because these are interpreted in a way which 
  //depends on "coeff_style").

  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]) {
        PairLJRemix::write_restart_pair(fp,
                                        i, j, 
                                        epsilon, sigma, 
                                        kappa, lambda,
                                        lj3, lj4);

        PairLJRemix::write_restart_pair(fp,
                                        i, j, 
                                        eps14, sigma14, 
                                        kappa14, lambda14,
                                        lj14_3, lj14_4);

        #ifdef ENABLE_INTER_DIF
        PairLJRemix::write_restart_pair(fp,
                                        i, j, 
                                        epsDIF, sigmaDIF, 
                                        kappaDIF, lambdaDIF,
                                        ljDIF_3, ljDIF_4);
        #endif


        #ifdef ENABLE_SOFTCORE2
        PairLJSoftcore2::write_restart_pair(fp,
                                            i, j,
                                            cut_core_sqd,
                                            core_coeff,
                                            core_D_coeff,
                                            core_offset);
        #ifdef ENABLE_INTER_DIF
        PairLJSoftcore2::write_restart_pair(fp,
                                            i, j,
                                            cut_coreDIF_sqd,
                                            core_coeffDIF,
                                            core_D_coeffDIF,
                                            core_offsetDIF);
        #endif //#ifdef ENABLE_INTER_DIF
        #endif //#ifdef ENABLE_SOFTCORE2
        // There's no way to modify 14 behavior without changing 
        // dihedral_charmm.cpp.  Disabling the cutoff behavior for 1-4 
        // interactions, and commenting out:
        //PairLJSoftcore2::write_restart_pair(fp,
        //                                    i, j,
        //                                    cut_core14_sqd,
        //                                    core_coeff14,
        //                                    core_D_coeff14,
        //                                    core_offset14);

      } //if (setflag[i][j])
    } //for (j = i; j <= atom->ntypes; j++)
  } //for (i = 1; i <= atom->ntypes; i++)
} //PairLJCharmmCoulCharmmInter::write_restart(FILE *)



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


void PairLJCharmmCoulCharmmInter::read_restart(FILE *fp)
{
  read_restart_settings(fp);

  allocate();

  //Just read the raw lennard jones AB parameters.  Don't read the epsilon,
  //sigma, and lambda parameters (because these are interpreted in a way which 
  //depends on "coeff_style").

  int i,j;
  int me = comm->me;
  for (i = 1; i <= atom->ntypes; i++) {
    for (j = i; j <= atom->ntypes; j++) {
      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
      if (setflag[i][j]) {
          PairLJRemix::read_restart_pair(fp,
                                         i, j,
                                         epsilon, sigma, kappa, lambda,
                                         lj3, lj4);
          PairLJRemix::read_restart_pair(fp,
                                         i, j,
                                         eps14, sigma14, kappa14, lambda14,
                                         lj14_3, lj14_4);
          #ifdef ENABLE_INTER_DIF
          PairLJRemix::read_restart_pair(fp,
                                         i, j,
                                         epsDIF, sigmaDIF, kappaDIF, lambdaDIF,
                                         ljDIF_3, ljDIF_4);
          #endif //#ifdef ENABLE_INTER_DIF



          #ifdef ENABLE_SOFTCORE2
          PairLJSoftcore2::read_restart_pair(fp,
                                             i, j,
                                             cut_core_sqd,
                                             core_coeff,
                                             core_D_coeff,
                                             core_offset);

          #ifdef ENABLE_INTER_DIF
          PairLJSoftcore2::read_restart_pair(fp,
                                             i, j,
                                             cut_coreDIF_sqd,
                                             core_coeffDIF,
                                             core_D_coeffDIF,
                                             core_offsetDIF);
          #endif
          // There's no way to modify 14 behavior without changing 
          // dihedral_charmm.cpp.  Disabling the cutoff behavior for 1-4 
          // interactions, and commenting out:
          //PairLJSoftcore2::read_restart_pair(fp,
          //                                   i, j,
          //                           cut_coreDIF_sqd,
          //                           core_coeffDIF,
          //                           core_D_coeffDIF,
          //                           core_offsetDIF);
          #endif //#ifdef ENABLE_SOFTCORE2

      } //if (setflag[i][j])
    } //for (j = i; j <= atom->ntypes; j++)
  } //for (i = 1; i <= atom->ntypes; i++)
} //PairLJCharmmCoulCharmmInter::read_restart(FILE *)


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

void PairLJCharmmCoulCharmmInter::write_restart_settings(FILE *fp)
{
  PairLJRemix::write_restart_settings(fp);
  #ifdef ENABLE_SOFTCORE2
  PairLJSoftcore2::write_restart_settings(fp);
  #endif

  fwrite(&cut_lj_inner,sizeof(double),1,fp);
  fwrite(&cut_lj,sizeof(double),1,fp);
  fwrite(&cut_coul_inner,sizeof(double),1,fp);
  fwrite(&cut_coul,sizeof(double),1,fp);
  fwrite(&offset_flag,sizeof(int),1,fp);
  fwrite(&mix_flag,sizeof(int),1,fp);
  fwrite(&cut_coul_and_lj_differ,sizeof(int),1,fp);
}


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

void PairLJCharmmCoulCharmmInter::read_restart_settings(FILE *fp)
{
  PairLJRemix::read_restart_settings(fp);
  #ifdef ENABLE_SOFTCORE2
  PairLJSoftcore2::read_restart_settings(fp);
  #endif

  int me = comm->me;
  if (me == 0) {
    fread(&cut_lj_inner,sizeof(double),1,fp);
    fread(&cut_lj,sizeof(double),1,fp);
    fread(&cut_coul_inner,sizeof(double),1,fp);
    fread(&cut_coul,sizeof(double),1,fp);
    fread(&offset_flag,sizeof(int),1,fp);
    fread(&mix_flag,sizeof(int),1,fp);
    fread(&cut_coul_and_lj_differ,sizeof(int),1,fp);
  } //if (me == 0)
  MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&cut_coul_inner,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
  MPI_Bcast(&cut_coul_and_lj_differ,1,MPI_INT,0,world);
} //PairLJSoftcore2::read_restart_pair()



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

double PairLJCharmmCoulCharmmInter::single(int i, int j, int itype, int jtype,
                                           double rsq,
                                           double factor_coul, double factor_lj,
                                           double &fforce)
{
  double r2inv,r6inv,forcecoul,forcelj,phicoul,philj;
  double switch1_coul, switch2_coul, switch1_lj, switch2_lj;

  double _lj1, _lj2, _lj3, _lj4;
  double _cut_core_sqd, _core_coeff, _core_D_coeff, _core_offset;
  #ifdef ENABLE_INTER_DIF
  if (atom->molecule[i] == atom->molecule[j]) {
  #endif

    _lj1 = lj1[itype][jtype];
    _lj2 = lj2[itype][jtype];
    _lj3 = lj3[itype][jtype];
    _lj4 = lj4[itype][jtype];

    #ifdef ENABLE_SOFTCORE2
    _cut_core_sqd = cut_core_sqd[itype][jtype];
    _core_coeff   =   core_coeff[itype][jtype];
    _core_D_coeff = core_D_coeff[itype][jtype];
    _core_offset  =  core_offset[itype][jtype];
    #endif

  #ifdef ENABLE_INTER_DIF
  }
  else {
    _lj1 = ljDIF_1[itype][jtype];
    _lj2 = ljDIF_2[itype][jtype];
    _lj3 = ljDIF_3[itype][jtype];
    _lj4 = ljDIF_4[itype][jtype];
    #ifdef ENABLE_SOFTCORE2
    _cut_core_sqd = cut_coreDIF_sqd[itype][jtype];
    _core_coeff   =   core_coeffDIF[itype][jtype];
    _core_D_coeff = core_D_coeffDIF[itype][jtype];
    _core_offset  =  core_offsetDIF[itype][jtype];
    #endif //#ifdef ENABLE_SOFTCORE2
  }
  #endif //#ifdef ENABLE_INTER_DIF


  r2inv = 1.0/rsq;
  if (rsq < cut_coulsq) {
    forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
    phicoul   = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
    if (rsq > cut_coul_innersq) {
      double cutcoulsq_minus_rsq = (cut_coulsq-rsq);
      double cutcoulsq_minus_rsq_sq = cutcoulsq_minus_rsq * cutcoulsq_minus_rsq;
      switch1_coul = cutcoulsq_minus_rsq_sq *
          (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) * denom_coul_inv;
      switch2_coul = 12.0*rsq * cutcoulsq_minus_rsq * 
          (rsq-cut_coul_innersq) * denom_coul_inv;
      forcecoul *= switch1_coul + switch2_coul;
      phicoul   *= switch1_coul;
    }
  } else forcecoul = 0.0;

  if (rsq < cut_ljsq) {
    r6inv = r2inv*r2inv*r2inv;
    #ifdef ENABLE_SOFTCORE2
    if (rsq > _cut_core_sqd) { //--> Ulj(r)=(_lj3/r^12- _lj4/r^6)*switch1_lj
      forcelj = r6inv * (_lj1*r6inv - _lj2);
      philj   = r6inv * (_lj3*r6inv - _lj4);
    }
    else { // otherwise Ulj(r) = (_core_coeff/r^2  + _core_offset)*switch1lj
    #endif
      forcelj = r2inv * _core_D_coeff;
      philj   = r2inv * _core_coeff + _core_offset;
      //fprintf(stderr,"ACCESS ERROR: reached %s:%d\n",__FILE__,__LINE__);
      //exit(-1);
    #ifdef ENABLE_SOFTCORE2
    } //else clause for "if (rsq > _cut_core_sqd)"
    #endif
    if (rsq > cut_lj_innersq) {
      double switch1_lj = switch1_coul;
      double switch2_lj = switch2_coul;
      if (cut_coul_and_lj_differ) { // then recompute switch1 & switch2
        double cutljsq_minus_rsq = (cut_ljsq-rsq);
        double cutljsq_minus_rsq_sq = cutljsq_minus_rsq * cutljsq_minus_rsq;
        switch1_lj = cutljsq_minus_rsq_sq *
                  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
        switch2_lj = 12.0*rsq * cutljsq_minus_rsq * 
                  (rsq-cut_lj_innersq) * denom_lj_inv;
      }
      forcelj = forcelj*switch1_lj + philj*switch2_lj;
      philj  *= switch1_lj;
    }
  } else forcelj = 0.0;

  fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
  double eng = factor_coul*phicoul + factor_lj*philj;
  return eng;
} //PairLJCharmmCoulCharmmInter::single()

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

void *PairLJCharmmCoulCharmmInter::extract(const char *str, int &dim)
{
  dim = 2;
  if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
  if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
  if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
  if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
  // There's no way to modify 14 behavior without changing dihedral_charmm.cpp.
  // Disabling the cutoff behavior for 1-4 interactions, and commenting out:
  //if (strcmp(str,"cut_core14_sqd") == 0) return (void *) cut_core14_sqd;
  //if (strcmp(str,"core_coeff14") == 0) return (void *) core_coeff14;
  //if (strcmp(str,"core_D_coeff14") == 0) return (void *) core_D_coeff14;
  //if (strcmp(str,"core_offset14") == 0) return (void *) core_offset14;

  dim = 0;
  if (strcmp(str,"implicit") == 0) return (void *) &implicit;

  return NULL;
}