lammps-sys 0.6.0

/* ----------------------------------------------------------------------
   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: Ray Shan (Sandia, tnshan@sandia.gov)
   Updates and debug: Tao Liang (U Florida, liang75@ufl.edu)
                      Dundar Yilmaz (dundar.yilmaz@zirve.edu.tr)
------------------------------------------------------------------------- */

#include "pair_comb3.h"
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "my_page.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "group.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace MathConst;

#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
#define MAXNEIGH 24

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

PairComb3::PairComb3(LAMMPS *lmp) : Pair(lmp)
{
  single_enable = 0;
  restartinfo = 0;
  one_coeff = 1;
  ghostneigh = 1;

  nmax = 0;
  NCo = NULL;
  bbij = NULL;
  map = NULL;
  esm = NULL;

  nelements = 0;
  elements = NULL;
  nparams = 0;
  maxparam = 0;
  params = NULL;
  elem2param = NULL;

  intype = NULL;
  afb = NULL;
  dafb = NULL;
  fafb = NULL;
  dfafb = NULL;
  ddfafb = NULL;
  phin = NULL;
  dphin = NULL;
  erpaw = NULL;
  vvdw = NULL;
  vdvdw = NULL;
  dpl = NULL;
  xcctmp = NULL;
  xchtmp = NULL;
  xcotmp = NULL;

  sht_num = NULL;
  sht_first = NULL;

  ipage = NULL;
  pgsize = oneatom = 0;

  cflag = 0;

  // set comm size needed by this Pair
  comm_forward = 1;
  comm_reverse = 1;
}

/* ----------------------------------------------------------------------
   check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */

PairComb3::~PairComb3()
{
  memory->destroy(NCo);

  if (elements)
    for (int i = 0; i < nelements; i++) delete [] elements[i];

  delete [] elements;
  memory->sfree(params);
  memory->destroy(elem2param);

  memory->destroy(afb);
  memory->destroy(dpl);
  memory->destroy(dafb);
  memory->destroy(fafb);
  memory->destroy(phin);
  memory->destroy(bbij);
  memory->destroy(vvdw);
  memory->destroy(vdvdw);
  memory->destroy(dphin);
  memory->destroy(erpaw);
  memory->destroy(dfafb);
  memory->destroy(ddfafb);
  memory->destroy(xcctmp);
  memory->destroy(xchtmp);
  memory->destroy(xcotmp);
  memory->destroy(intype);
  memory->destroy(sht_num);
  memory->sfree(sht_first);

  delete [] ipage;

  if (allocated) {
    memory->destroy(setflag);
    memory->destroy(cutsq);
    memory->destroy(cutghost);
    delete [] map;
    delete [] esm;
  }

}

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

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

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

 map = new int[n+1];
 esm = new double[n];
}

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

void PairComb3::settings(int narg, char **arg)
{

  if (narg != 1) error->all(FLERR,"Illegal pair_style command");

  if (strcmp(arg[0],"polar_on") == 0) {
    pol_flag = 1;
    if (comm->me == 0 && screen) fprintf(screen,
                    "   PairComb3: polarization is on \n");
  } else if (strcmp(arg[0],"polar_off") == 0) {
    pol_flag = 0;
    if (comm->me == 0 && screen) fprintf(screen,
                    "   PairComb3: polarization is off \n");
  } else {
    error->all(FLERR,"Illegal pair_style command");
  }
}

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

void PairComb3::coeff(int narg, char **arg)
{
  int i,j,n;

  if (!allocated) allocate();

  if (narg != 3 + atom->ntypes)
    error->all(FLERR,"Incorrect args for pair coefficients");

  // insure I,J args are * *
  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
    error->all(FLERR,"Incorrect args for pair coefficients");

  // read args that map atom types to elements in potential file
  // map[i] = which element the Ith atom type is, -1 if NULL
  // nelements = # of unique elements
  // elements = list of element names

  if (elements) {
    for (i = 0; i < nelements; i++) delete [] elements[i];
    delete [] elements;
  }
  elements = new char*[atom->ntypes];
  for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;

  nelements = 0;
  for (i = 3; i < narg; i++) {
    if ((strcmp(arg[i],"C") == 0) && (cflag == 0)) {
      if (comm->me == 0 && screen) fprintf(screen,
      " PairComb3: Found C: reading additional library file\n");
    read_lib();
    cflag = 1;
    }
    if (strcmp(arg[i],"NULL") == 0) {
      map[i-2] = -1;
      continue;
    }
    for (j = 0; j < nelements; j++)
      if (strcmp(arg[i],elements[j]) == 0) break;
    map[i-2] = j;
    if (j == nelements) {
      n = strlen(arg[i]) + 1;
      elements[j] = new char[n];
      strcpy(elements[j],arg[i]);
      nelements++;
    }
  }

  // read potential file and initialize potential parameters

  read_file(arg[2]);
  setup_params();

  n = atom->ntypes;

  // generate Wolf 1/r energy and van der Waals look-up tables
  tables();

  // clear setflag since coeff() called once with I,J = * *
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      setflag[i][j] = 0;

  // set setflag i,j for type pairs where both are mapped to elements
  int count = 0;
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      if (map[i] >= 0 && map[j] >= 0) {
        setflag[i][j] = 1;
        count++;
      }

  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");

}

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

void PairComb3::init_style()
{
  if (atom->tag_enable == 0)
    error->all(FLERR,"Pair style COMB3 requires atom IDs");
  if (force->newton_pair == 0)
    error->all(FLERR,"Pair style COMB3 requires newton pair on");
  if (!atom->q_flag)
    error->all(FLERR,"Pair style COMB3 requires atom attribute q");

  // need a full neighbor list
  int irequest = neighbor->request(this,instance_me);
  neighbor->requests[irequest]->half = 0;
  neighbor->requests[irequest]->full = 1;
  neighbor->requests[irequest]->ghost = 1;

  // local Comb neighbor list
  // create pages if first time or if neighbor pgsize/oneatom has changed

  int create = 0;
  if (ipage == NULL) create = 1;
  if (pgsize != neighbor->pgsize) create = 1;
  if (oneatom != neighbor->oneatom) create = 1;

  if (create) {
    delete [] ipage;
    pgsize = neighbor->pgsize;
    oneatom = neighbor->oneatom;

    int nmypage = comm->nthreads;
    ipage = new MyPage<int>[nmypage];
    for (int i = 0; i < nmypage; i++)
      ipage[i].init(oneatom,pgsize);
  }
}

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

double PairComb3::init_one(int i, int j)
{
  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
  cutghost[j][i] = cutghost[i][j] = cutmax;
  return cutmax;
}

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

void PairComb3::read_lib()
{
  const unsigned int MAXLIB = 1024;
  int i,j,k,l,nwords,m;
  int ii,jj,kk,ll,mm,iii;
  char s[MAXLIB];
  char **words = new char*[80];

  // open libraray file on proc 0

  if (comm->me == 0) {
    FILE *fp = force->open_potential("lib.comb3");
    if (fp == NULL) error->one(FLERR,"Cannot open COMB3 lib.comb3 file");

    // read and store at the same time
    fgets(s,MAXLIB,fp);
    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    ccutoff[0] = atof(words[0]);
    ccutoff[1] = atof(words[1]);
    ccutoff[2] = atof(words[2]);
    ccutoff[3] = atof(words[3]);
    ccutoff[4] = atof(words[4]);
    ccutoff[5] = atof(words[5]);

    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    ch_a[0] = atof(words[0]);
    ch_a[1] = atof(words[1]);
    ch_a[2] = atof(words[2]);
    ch_a[3] = atof(words[3]);
    ch_a[4] = atof(words[4]);
    ch_a[5] = atof(words[5]);
    ch_a[6] = atof(words[6]);

    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    nsplpcn = atoi(words[0]);
    nsplrad = atoi(words[1]);
    nspltor = atoi(words[2]);

    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    maxx = atoi(words[0]);
    maxy = atoi(words[1]);
    maxz = atoi(words[2]);

    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    maxxc = atoi(words[0]);
    maxyc = atoi(words[1]);
    maxconj = atoi(words[2]);

    for (l=0; l<nsplpcn; l++) {
      fgets(s,MAXLIB,fp);
      nwords = 0;
      words[nwords++] = strtok(s," \t\n\r\f");
      while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
      maxxcn[l] = atoi(words[1]);
      vmaxxcn[l] = atof(words[2]);
      dvmaxxcn[l] = atof(words[3]);
    }

    fgets(s,MAXLIB,fp);
    nwords = 0;
    words[nwords++] = strtok(s," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
    ntab = atoi(words[0]);

    for (i=0; i<ntab+1; i++){
      fgets(s,MAXLIB,fp);
      nwords = 0;
      words[nwords++] = strtok(s," \t\n\r\f");
      while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
      pang[i]    = atof(words[1]);
      dpang[i]   = atof(words[2]);
      ddpang[i]  = atof(words[3]);
    }

    for (l=0; l<nsplpcn; l++)
      for (i=0; i<maxx+1; i++)
        for (j=0; j<maxy+1; j++)
          for (k=0; k<maxz+1; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
            ll = atoi(words[0])-1;
            ii = atoi(words[1]);
            jj = atoi(words[2]);
            kk = atoi(words[3]);
            pcn_grid[ll][ii][jj][kk]     = atof(words[4]);
            pcn_gridx[ll][ii][jj][kk]    = atof(words[5]);
            pcn_gridy[ll][ii][jj][kk]    = atof(words[6]);
            pcn_gridz[ll][ii][jj][kk]    = atof(words[7]);
          }

    for (l=0; l<nsplpcn; l++)
      for (i=0; i<maxx; i++)
        for (j=0; j<maxy; j++)
          for (k=0; k<maxz; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
           ll = atoi(words[0])-1;
           ii = atoi(words[1]);
           jj = atoi(words[2]);
           kk = atoi(words[3]);
           for(iii=0; iii<2; iii++) {
             fgets(s,MAXLIB,fp);
             nwords = 0;
             words[nwords++] = strtok(s," \t\n\r\f");
             while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
             for(m=0; m<32 ; m++) {
               mm=iii*32+m;
               pcn_cubs[ll][ii][jj][kk][mm] = atof(words[m]);
             }
           }
          }

    for (l=0; l<nsplrad; l++)
      for (i=0; i<maxxc+1; i++)
        for (j=0; j<maxyc+1; j++)
          for (k=0; k<maxconj; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
            ll = atoi(words[0])-1;
            ii = atoi(words[1]);
            jj = atoi(words[2]);
            kk = atoi(words[3])-1;
            rad_grid[ll][ii][jj][kk]     = atof(words[4]);
            rad_gridx[ll][ii][jj][kk]    = atof(words[5]);
            rad_gridy[ll][ii][jj][kk]    = atof(words[6]);
            rad_gridz[ll][ii][jj][kk]    = atof(words[7]);
          }

    for (l=0; l<nsplrad; l++)
      for (i=0; i<maxxc; i++)
        for (j=0; j<maxyc; j++)
          for (k=0; k<maxconj-1; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
            ll = atoi(words[0])-1;
            ii = atoi(words[1]);
            jj = atoi(words[2]);
            kk = atoi(words[3])-1;
            for (iii=0; iii<2; iii++) {
              fgets(s,MAXLIB,fp);
              nwords = 0;
              words[nwords++] = strtok(s," \t\n\r\f");
              while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
              for(m=0; m<32 ; m++){
                mm=iii*32+m;
                rad_spl[ll][ii][jj][kk][mm] = atof(words[m]);
              }
            }
          }

    for (l=0; l<nspltor; l++)
      for (i=0; i<maxxc+1; i++)
        for (j=0; j<maxyc+1; j++)
          for (k=0; k<maxconj; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
            ll = atoi(words[0])-1;
            ii = atoi(words[1]);
            jj = atoi(words[2]);
            kk = atoi(words[3])-1;
            tor_grid[ll][ii][jj][kk]     = atof(words[4]);
            tor_gridx[ll][ii][jj][kk]    = atof(words[5]);
            tor_gridy[ll][ii][jj][kk]    = atof(words[6]);
            tor_gridz[ll][ii][jj][kk]    = atof(words[7]);
          }

    for (l=0; l<nspltor; l++)
      for (i=0; i<maxxc; i++)
        for (j=0; j<maxyc; j++)
          for (k=0; k<maxconj-1; k++) {
            fgets(s,MAXLIB,fp);
            nwords = 0;
            words[nwords++] = strtok(s," \t\n\r\f");
            while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
            ll = atoi(words[0])-1;
            ii = atoi(words[1]);
            jj = atoi(words[2]);
            kk = atoi(words[3])-1;
            for(iii=0; iii<2; iii++) {
              fgets(s,MAXLIB,fp);
              nwords = 0;
              words[nwords++] = strtok(s," \t\n\r\f");
              while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
              for (m=0; m<32 ; m++){
                mm=iii*32+m;
                tor_spl[ll][ii][jj][kk][mm] = atof(words[m]);
              }
            }
          }

    fclose(fp);
  }

  k = 0;
  for (i=0; i<4; i++)
    for (j=0; j<4; j++) {
      iin2[k][0] = i;
      iin2[k][1] = j;
      k ++;
    }

  l = 0;
  for (i=0; i<4; i++)
    for (j=0; j<4; j++)
      for (k=0; k<4; k++) {
        iin3[l][0] = i;
        iin3[l][1] = j;
        iin3[l][2] = k;
        l ++;
      }

  MPI_Bcast(&ccutoff[0],6,MPI_DOUBLE,0,world);
  MPI_Bcast(&ch_a[0],7,MPI_DOUBLE,0,world);
  MPI_Bcast(&nsplpcn,1,MPI_INT,0,world);
  MPI_Bcast(&nsplrad,1,MPI_INT,0,world);
  MPI_Bcast(&nspltor,1,MPI_INT,0,world);
  MPI_Bcast(&maxx,1,MPI_INT,0,world);
  MPI_Bcast(&maxy,1,MPI_INT,0,world);
  MPI_Bcast(&maxz,1,MPI_INT,0,world);
  MPI_Bcast(&maxxc,1,MPI_INT,0,world);
  MPI_Bcast(&maxyc,1,MPI_INT,0,world);
  MPI_Bcast(&maxconj,1,MPI_INT,0,world);
  MPI_Bcast(&maxxcn,4,MPI_INT,0,world);
  MPI_Bcast(&vmaxxcn,4,MPI_DOUBLE,0,world);
  MPI_Bcast(&dvmaxxcn,4,MPI_DOUBLE,0,world);
  MPI_Bcast(&ntab,1,MPI_INT,0,world);
  MPI_Bcast(&pang[0],20001,MPI_DOUBLE,0,world);
  MPI_Bcast(&dpang[0],20001,MPI_DOUBLE,0,world);
  MPI_Bcast(&ddpang[0],20001,MPI_DOUBLE,0,world);
  MPI_Bcast(&pcn_grid[0][0][0][0],500,MPI_DOUBLE,0,world);
  MPI_Bcast(&pcn_gridx[0][0][0][0],500,MPI_DOUBLE,0,world);
  MPI_Bcast(&pcn_gridy[0][0][0][0],500,MPI_DOUBLE,0,world);
  MPI_Bcast(&pcn_gridz[0][0][0][0],500,MPI_DOUBLE,0,world);
  MPI_Bcast(&pcn_cubs[0][0][0][0][0],16384,MPI_DOUBLE,0,world);

  MPI_Bcast(&rad_grid[0][0][0][0],825,MPI_DOUBLE,0,world);
  MPI_Bcast(&rad_gridx[0][0][0][0],825,MPI_DOUBLE,0,world);
  MPI_Bcast(&rad_gridy[0][0][0][0],825,MPI_DOUBLE,0,world);
  MPI_Bcast(&rad_gridz[0][0][0][0],825,MPI_DOUBLE,0,world);
  MPI_Bcast(&rad_spl[0][0][0][0][0],30720,MPI_DOUBLE,0,world);

  MPI_Bcast(&tor_grid[0][0][0][0],275,MPI_DOUBLE,0,world);
  MPI_Bcast(&tor_gridx[0][0][0][0],275,MPI_DOUBLE,0,world);
  MPI_Bcast(&tor_gridy[0][0][0][0],275,MPI_DOUBLE,0,world);
  MPI_Bcast(&tor_gridz[0][0][0][0],275,MPI_DOUBLE,0,world);
  MPI_Bcast(&tor_spl[0][0][0][0][0],10240,MPI_DOUBLE,0,world);

  MPI_Bcast(&iin2[0][0],32,MPI_INT,0,world);
  MPI_Bcast(&iin3[0][0],192,MPI_INT,0,world);
  delete [] words;
}

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

void PairComb3::read_file(char *file)
{
  int params_per_line = 74;
  char **words = new char*[params_per_line+1];

  if (params) delete [] params;
  params = NULL;
  nparams = 0;

  // open file on proc 0

  FILE *fp;
  if (comm->me == 0) {
    fp = force->open_potential(file);
    if (fp == NULL) {
      char str[128];
      snprintf(str,128,"Cannot open COMB3 potential file %s",file);
      error->one(FLERR,str);
    }
  }

  // read each line out of file, skipping blank lines or leading '#'
  // store line of params if all 3 element tags are in element list

  int n,nwords,ielement,jelement,kelement;
  char line[MAXLINE],*ptr;
  int eof = 0;
  nwords=0;
  while (1) {
    if (comm->me == 0) {
      ptr = fgets(line,MAXLINE,fp);
      if (ptr == NULL) {
        eof = 1;
        fclose(fp);
      } else n = strlen(line) + 1;
    }
    MPI_Bcast(&eof,1,MPI_INT,0,world);
    if (eof) break;
    MPI_Bcast(&n,1,MPI_INT,0,world);
    MPI_Bcast(line,n,MPI_CHAR,0,world);

    // strip comment, skip line if blank

    if ((ptr = strchr(line,'#'))) *ptr = '\0';

    nwords = atom->count_words(line);
    if (nwords == 0) continue;

    // concatenate additional lines until have params_per_line words

    while (nwords < params_per_line) {
      n = strlen(line);

      if (comm->me == 0) {
        ptr = fgets(&line[n],MAXLINE-n,fp);
        if (ptr == NULL) {
          eof = 1;
          fclose(fp);
        } else n = strlen(line) + 1;
      }
      MPI_Bcast(&eof,1,MPI_INT,0,world);
      if (eof) break;
      MPI_Bcast(&n,1,MPI_INT,0,world);
      MPI_Bcast(line,n,MPI_CHAR,0,world);
      if ((ptr = strchr(line,'#'))) *ptr = '\0';
      nwords = atom->count_words(line);
    }
    if (nwords != params_per_line){
      error->all(FLERR,"Incorrect format in COMB3 potential file");
}
    // words = ptrs to all words in line

    nwords = 0;
    words[nwords++] = strtok(line," \t\n\r\f");
    while ((nwords <= params_per_line)
           && (words[nwords++] = strtok(NULL," \t\n\r\f"))) {
      continue;
    }

    // ielement,jelement,kelement = 1st args
    // if all 3 args are in element list, then parse this line
    // else skip to next line

    for (ielement = 0; ielement < nelements; ielement++)
      if (strcmp(words[0],elements[ielement]) == 0) break;
    if (ielement == nelements) continue;
    for (jelement = 0; jelement < nelements; jelement++)
      if (strcmp(words[1],elements[jelement]) == 0) break;
    if (jelement == nelements) continue;
    for (kelement = 0; kelement < nelements; kelement++)
      if (strcmp(words[2],elements[kelement]) == 0) break;
    if (kelement == nelements) continue;

    // load up parameter settings and error check their values

    if (nparams == maxparam) {
      maxparam += DELTA;
      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
                                          "pair:params");
    }

    params[nparams].ielement = ielement;
    params[nparams].jelement = jelement;
    params[nparams].kelement = kelement;
    params[nparams].ielementgp = atoi(words[3]);
    params[nparams].jelementgp = atoi(words[4]);
    params[nparams].kelementgp = atoi(words[5]);
    params[nparams].ang_flag = atoi(words[6]);
    params[nparams].pcn_flag = atoi(words[7]);
    params[nparams].rad_flag = atoi(words[8]);
    params[nparams].tor_flag = atoi(words[9]);
    params[nparams].vdwflag = atof(words[10]);
    params[nparams].powerm = atof(words[11]);
    params[nparams].veps = atof(words[12]);
    params[nparams].vsig = atof(words[13]);
    params[nparams].paaa = atof(words[14]);
    params[nparams].pbbb = atof(words[15]);
    params[nparams].lami = atof(words[16]);
    params[nparams].alfi = atof(words[17]);
    params[nparams].powern = atof(words[18]);
    params[nparams].QL = atof(words[19]);
    params[nparams].QU = atof(words[20]);
    params[nparams].DL = atof(words[21]);
    params[nparams].DU = atof(words[22]);
    params[nparams].qmin = atof(words[23]);
    params[nparams].qmax = atof(words[24]);
    params[nparams].chi = atof(words[25]);
    params[nparams].dj  = atof(words[26]);
    params[nparams].dk  = atof(words[27]);
    params[nparams].dl  = atof(words[28]);
    params[nparams].esm = atof(words[29]);
    params[nparams].cmn1 = atof(words[30]);
    params[nparams].cmn2 = atof(words[31]);
    params[nparams].pcmn1 = atof(words[32]);
    params[nparams].pcmn2 = atof(words[33]);
    params[nparams].coulcut = atof(words[34]);
    params[nparams].polz = atof(words[35]);
    params[nparams].curl = atof(words[36]);
    params[nparams].curlcut1 = atof(words[37]);
    params[nparams].curlcut2 = atof(words[38]);
    params[nparams].curl0 = atof(words[39]);
    params[nparams].alpha1 = atof(words[40]);
    params[nparams].bigB1 = atof(words[41]);
    params[nparams].alpha2 = atof(words[42]);
    params[nparams].bigB2 = atof(words[43]);
    params[nparams].alpha3 = atof(words[44]);
    params[nparams].bigB3 = atof(words[45]);
    params[nparams].lambda = atof(words[46]);
    params[nparams].bigA = atof(words[47]);
    params[nparams].beta = atof(words[48]);
    params[nparams].bigr = atof(words[49]);
    params[nparams].bigd = atof(words[50]);
    params[nparams].pcos6 = atof(words[51]);
    params[nparams].pcos5 = atof(words[52]);
    params[nparams].pcos4 = atof(words[53]);
    params[nparams].pcos3 = atof(words[54]);
    params[nparams].pcos2 = atof(words[55]);
    params[nparams].pcos1 = atof(words[56]);
    params[nparams].pcos0 = atof(words[57]);
    params[nparams].pcna = atof(words[58]);
    params[nparams].pcnb = atof(words[59]);
    params[nparams].pcnc = atof(words[60]);
    params[nparams].pcnd = atof(words[61]);
    params[nparams].p6p0 = atof(words[62]);
    params[nparams].p6p1 = atof(words[63]);
    params[nparams].p6p2 = atof(words[64]);
    params[nparams].p6p3 = atof(words[65]);
    params[nparams].p6p4 = atof(words[66]);
    params[nparams].p6p5 = atof(words[67]);
    params[nparams].p6p6 = atof(words[68]);
    params[nparams].ptork1=atof(words[69]);
    params[nparams].ptork2=atof(words[70]);
    params[nparams].addrepr=atof(words[71]);
    params[nparams].addrep=atof(words[72]);
    params[nparams].pcross = atof(words[73]);
    params[nparams].powermint = int(params[nparams].powerm);

    // parameter sanity checks

    if (params[nparams].lambda < 0.0 || params[nparams].powern < 0.0 ||
        params[nparams].beta < 0.0 || params[nparams].alpha1 < 0.0 ||
        params[nparams].bigB1< 0.0 || params[nparams].bigA< 0.0 ||
        params[nparams].bigB2< 0.0 || params[nparams].alpha2 <0.0 ||
        params[nparams].bigB3< 0.0 || params[nparams].alpha3 <0.0 ||
        params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
        params[nparams].bigd > params[nparams].bigr ||
        params[nparams].powerm - params[nparams].powermint != 0.0 ||
        params[nparams].addrepr < 0.0 || params[nparams].powermint < 1.0 ||
        params[nparams].QL > 0.0 || params[nparams].QU < 0.0 ||
        params[nparams].DL < 0.0 || params[nparams].DU > 0.0 ||
        params[nparams].pcross < 0.0 ||
        params[nparams].esm < 0.0 || params[nparams].veps < 0.0 ||
        params[nparams].vsig < 0.0 || params[nparams].vdwflag < 0.0
        )
      error->all(FLERR,"Illegal COMB3 parameter");

    nparams++;
  }

  delete [] words;
}

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

void PairComb3::setup_params()
{
  int i,j,k,m,n;

  // set elem2param for all element triplet combinations
  // must be a single exact match to lines read from file
  // do not allow for ACB in place of ABC

  memory->destroy(elem2param);
  memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param");

  for (i = 0; i < nelements; i++)
    for (j = 0; j < nelements; j++)
      for (k = 0; k < nelements; k++) {
        n = -1;
        for (m = 0; m < nparams; m++) {
          if (i == params[m].ielement && j == params[m].jelement &&
              k == params[m].kelement) {
            if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
            n = m;
          }
        }
        if (n < 0) error->all(FLERR,"Potential file is missing an entry");
        elem2param[i][j][k] = n;
      }

  // compute parameter values derived from inputs

  for (m = 0; m < nparams; m++) {
    params[m].cut = params[m].bigr + params[m].bigd;
    params[m].cutsq = params[m].cut*params[m].cut;
    params[m].c1 = pow(2.0*params[m].powern*1.0e-16,-1.0/params[m].powern);
    params[m].c2 = pow(2.0*params[m].powern*1.00e-8,-1.0/params[m].powern);
    params[m].c3 = 1.0/params[m].c2;
    params[m].c4 = 1.0/params[m].c1;

    params[m].Qo = (params[m].QU+params[m].QL)/2.0; // (A22)
    params[m].dQ = (params[m].QU-params[m].QL)/2.0; // (A21)
    params[m].aB = 1.0 /
      (1.0-pow(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
    params[m].bB = pow(fabs(params[m].aB),0.1)/params[m].dQ; // (A19)
    params[m].nD = log(params[m].DU/(params[m].DU-params[m].DL))/
                    log(params[m].QU/(params[m].QU-params[m].QL));
    params[m].bD = (pow((params[m].DL-params[m].DU),(1.0/params[m].nD)))/
                    (params[m].QU-params[m].QL);

    params[m].lcut = params[m].coulcut;
    params[m].lcutsq = params[m].lcut*params[m].lcut;
  }

  // set cutmax to max of all params

  cutmin = cutmax = 0.0;
  polar = 0;
  for (m = 0; m < nparams; m++) {
    if (params[m].cutsq > cutmin) cutmin = params[m].cutsq + 2.0;
    if (params[m].lcut > cutmax) cutmax = params[m].lcut;
  }
  chicut1 = 7.0;
  chicut2 = cutmax;
}

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

void PairComb3::Short_neigh()
{
  int nj,*neighptrj,icontrol;
  int iparam_ij,*ilist,*jlist,*numneigh,**firstneigh;
  int inum,jnum,i,j,ii,jj,itype,jtype;
  double rr1,rsq1,delrj[3];

  double **x = atom->x;
  int *type  = atom->type;

  if (atom->nmax > nmax) {
    memory->sfree(sht_first);
    nmax = atom->nmax;
    sht_first = (int **) memory->smalloc(nmax*sizeof(int *),
                                        "pair:sht_first");
    memory->grow(dpl,nmax,3,"pair:dpl");
    memory->grow(xcctmp,nmax,"pair:xcctmp");
    memory->grow(xchtmp,nmax,"pair:xchtmp");
    memory->grow(xcotmp,nmax,"pair:xcotmp");
    memory->grow(NCo,nmax,"pair:NCo");
    memory->grow(sht_num,nmax,"pair:sht_num");
    memory->grow(bbij,nmax,MAXNEIGH,"pair:bbij");
  }

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

  // create COMB neighbor list

  ipage->reset();

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    dpl[i][0] = dpl[i][1] = dpl[i][2] = 0.0;

    nj = 0;
    neighptrj = ipage->vget();

    itype = map[type[i]];
    jlist = firstneigh[i];
    jnum = numneigh[i];

    NCo[i] = 0.0;
    xcctmp[i] = 0.0;
    xchtmp[i] = 0.0;
    xcotmp[i] = 0.0;

    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj] & NEIGHMASK;

      delrj[0] = x[i][0] - x[j][0];
      delrj[1] = x[i][1] - x[j][1];
      delrj[2] = x[i][2] - x[j][2];
      rsq1 = vec3_dot(delrj,delrj);
      jtype = map[type[j]];
      iparam_ij = elem2param[itype][jtype][jtype];

      if (rsq1 > cutmin) continue;

      neighptrj[nj++] = j;
      rr1 = sqrt(rsq1);
      NCo[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;

      icontrol = params[iparam_ij].jelementgp;

      if (icontrol == 1)
          xcctmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
      if (icontrol == 2)
          xchtmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
      if (icontrol == 3)
          xcotmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;

    }

    sht_first[i] = neighptrj;
    sht_num[i] = nj;
    ipage->vgot(nj);
    if (ipage->status())
      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
  }

  // communicating coordination number to all nodes
  pack_flag = 2;
  comm->forward_comm_pair(this);

}

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

void PairComb3::compute(int eflag, int vflag)
{
  int i,ii,k,kk,j,jj,im,inum,jnum,itype,jtype,ktype;
  int iparam_i,iparam_ij,iparam_ji;
  int iparam_ijk,iparam_jik,iparam_ikj,iparam_jli,iparam_ikl;
  int sht_jnum,*sht_jlist,sht_lnum,*sht_llist;
  int sht_mnum,*sht_mlist,sht_pnum,*sht_plist;
  int *ilist,*jlist,*numneigh,**firstneigh,mr1,mr2,mr3,inty,nj;
  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
  double rsq,rsq1,rsq2,rsq3,iq,jq,yaself;
  double eng_tmp,vionij,fvionij,sr1,sr2,sr3;
  double zeta_ij,prefac_ij1,prefac_ij2,prefac_ij3,prefac_ij4,prefac_ij5;
  double zeta_ji,prefac_ji1,prefac_ji2,prefac_ji3,prefac_ji4,prefac_ji5;
  double delrj[3],delrk[3],fi[3],fj[3],fk[3],fl[3];
  double ep6p_ij,ep6p_ji,fip6p[3],fjp6p[3],fkp6p[3],flp6p[3];
  double potal,fac11,fac11e;
  tagint itag, jtag;

  int nlocal = atom->nlocal;
  int newton_pair = force->newton_pair;
  tagint *tag = atom->tag;
  int *type = atom->type;
  double **x = atom->x;
  double **f = atom->f;
  double *q = atom->q;

  // coordination terms
  double xcn, ycn;
  double kcn, lcn;
  int torindx;

  // torsion and radical variables
  int l, ll, ltype, m, mm, mtype, p, pp, ptype;
  int iparam_jil, iparam_ijl, iparam_ki, iparam_lj;
  int iparam_jl, iparam_ik, iparam_km, iparam_lp;
  double kconjug, lconjug, kradtot, lradtot;
  double delrl[3], delrm[3], delrp[3], ddprx[3], srmu;
  double zet_addi,zet_addj;

  evdwl = eng_tmp = 0.0;

  ev_init(eflag,vflag);

  // Build short range neighbor list
  Short_neigh();

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;
  yaself = vionij = fvionij = fpair = 0.0;

  // self energy correction term: potal
  potal_calc(potal,fac11,fac11e);

  // generate initial dipole tensor
  if (pol_flag )
    for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      itag = tag[i];
      itype = map[type[i]];
      xtmp = x[i][0];
      ytmp = x[i][1];
      ztmp = x[i][2];
      iq = q[i];
      jlist = firstneigh[i];
      jnum = numneigh[i];

      for (jj = 0; jj < jnum; jj++) {
        j = jlist[jj] & NEIGHMASK;
        jtag = tag[j];

        if (itag > jtag) {
          if ((itag+jtag) % 2 == 0) continue;
        } else if (itag < jtag) {
          if ((itag+jtag) % 2 == 1) continue;
        } else {
          if (x[j][2] < x[i][2]) continue;
          if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
          if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
        }

        jtype = map[type[j]];
        jq = q[j];

        delrj[0] = x[j][0] - xtmp;
        delrj[1] = x[j][1] - ytmp;
        delrj[2] = x[j][2] - ztmp;
        rsq = vec3_dot(delrj,delrj);

        iparam_ij = elem2param[itype][jtype][jtype];
        iparam_ji = elem2param[jtype][itype][itype];

        if (rsq > params[iparam_ij].lcutsq) continue;

        tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3);

        dipole_init(&params[iparam_ij],&params[iparam_ji],fac11,delrj,
                    rsq,mr1,mr2,mr3,sr1,sr2,sr3,iq,jq,i,j);
      }
    }

  // loop over full neighbor list of my atoms
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    itag = tag[i];
    itype = map[type[i]];
    xtmp = x[i][0];
    ytmp = x[i][1];
    ztmp = x[i][2];
    iq = q[i];
    nj = 0;
    iparam_i = elem2param[itype][itype][itype];

    // self energy, only on i atom
    yaself = self(&params[iparam_i],iq);

    // dipole self energy
    if (pol_flag)
      yaself += dipole_self(&params[iparam_i],i);

    if (evflag) ev_tally(i,i,nlocal,0,0.0,yaself,0.0,0.0,0.0,0.0);

    // two-body interactions (long:R + A, short: only R)

    jlist = firstneigh[i];
    jnum = numneigh[i];
    sht_jlist = sht_first[i];
    sht_jnum = sht_num[i];

    // long range interactions
    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj] & NEIGHMASK;

      jtag = tag[j];

      if (itag > jtag) {
        if ((itag+jtag) % 2 == 0) continue;
      } else if (itag < jtag) {
        if ((itag+jtag) % 2 == 1) continue;
      } else {
        if (x[j][2] < x[i][2]) continue;
        if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
        if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
      }

      jtype = map[type[j]];
      jq = q[j];

      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      rsq = delx*delx + dely*dely + delz*delz;
      iparam_ij = elem2param[itype][jtype][jtype];
      iparam_ji = elem2param[jtype][itype][itype];

      if (rsq > params[iparam_ij].lcutsq) continue;

      inty = intype[itype][jtype];

      // three-point interpolation
      tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3);

      // Q-indenpendent: van der Waals
      vdwaals(inty,mr1,mr2,mr3,rsq,sr1,sr2,sr3,eng_tmp,fpair);
      evdwl = eng_tmp;
      f[i][0] += delx*fpair;
      f[i][1] += dely*fpair;
      f[i][2] += delz*fpair;
      f[j][0] -= delx*fpair;
      f[j][1] -= dely*fpair;
      f[j][2] -= delz*fpair;

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

      // Q-dependent: Coulombic, field, polarization
      // 1/r energy and forces

      direct(&params[iparam_ij], &params[iparam_ji],
                mr1, mr2, mr3, rsq, sr1, sr2, sr3, iq, jq,
                fac11, fac11e, eng_tmp, fvionij, i, j);

      vionij = eng_tmp;

      // field correction to self energy
      field(&params[iparam_ij], &params[iparam_ji],rsq,iq,jq,
             eng_tmp,fvionij);
      vionij += eng_tmp;

      // sums up long range Q-dependent forces (excluding dipole)
      f[i][0] += delx*fvionij;
      f[i][1] += dely*fvionij;
      f[i][2] += delz*fvionij;
      f[j][0] -= delx*fvionij;
      f[j][1] -= dely*fvionij;
      f[j][2] -= delz*fvionij;

      // sums up long range Q-dependent energies (excluding dipole)
      if (evflag)
        ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);

      // polarization field
      if (pol_flag) {
        dipole_calc(&params[iparam_ij], &params[iparam_ji],fac11,
                delx,dely,delz,rsq,mr1,mr2,mr3,
                sr1,sr2,sr3,iq,jq,i,j,eng_tmp,fvionij,ddprx);
        vionij = eng_tmp;

        // sums up dipole energies
        if (evflag)
          ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);

        // sums up dipole forces
        f[i][0] += (ddprx[0] + delx*fvionij);
        f[i][1] += (ddprx[1] + dely*fvionij);
        f[i][2] += (ddprx[2] + delz*fvionij);
        f[j][0] -= (ddprx[0] + delx*fvionij);
        f[j][1] -= (ddprx[1] + dely*fvionij);
        f[j][2] -= (ddprx[2] + delz*fvionij);
      }

      if (rsq > params[iparam_ij].cutsq) continue;

      repulsive(&params[iparam_ij], &params[iparam_ji], rsq,
              fpair, eflag, eng_tmp, iq, jq);

      evdwl = eng_tmp;

      // repulsion is pure two-body, sums up pair repulsive forces
      f[i][0] += delx*fpair;
      f[i][1] += dely*fpair;
      f[i][2] += delz*fpair;
      f[j][0] -= delx*fpair;
      f[j][1] -= dely*fpair;
      f[j][2] -= delz*fpair;

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

    // many-body interactions: start of short-range
    xcn = NCo[i];
    for (jj = 0; jj < sht_jnum; jj++) {
      j = sht_jlist[jj];
      sht_llist = sht_first[j];
      sht_lnum = sht_num[j];

      jtag = tag[j];
      if ( jtag <= itag ) continue ;
      ycn = NCo[j];

      jtype = map[type[j]];
      iparam_ij = elem2param[itype][jtype][jtype];
      iparam_ji = elem2param[jtype][itype][itype];

      delrj[0] = x[j][0] - xtmp;
      delrj[1] = x[j][1] - ytmp;
      delrj[2] = x[j][2] - ztmp;
      rsq1 = vec3_dot(delrj,delrj);
      if (rsq1 > params[iparam_ij].cutsq) continue;
      nj ++;

      // this Qj for q-dependent BSi
      jq = q[j];

      // accumulate bondorder zeta for each i-j interaction via k and l loops
      zeta_ij = 0.0;
      bbtor = 0.0;
      kconjug = 0.0;

      for (kk = 0; kk < sht_jnum; kk++) {       // kk is neighbor of ii
        k = sht_jlist[kk];
        if (j == k) continue;

        ktype = map[type[k]];
        iparam_ijk = elem2param[itype][jtype][ktype];
        iparam_ikj = elem2param[itype][ktype][jtype];
        iparam_jik = elem2param[jtype][itype][ktype];
        iparam_ik  = elem2param[itype][ktype][ktype];
        delrk[0] = x[k][0] - xtmp;
        delrk[1] = x[k][1] - ytmp;
        delrk[2] = x[k][2] - ztmp;
        rsq2 = vec3_dot(delrk,delrk);

        if (rsq2 > params[iparam_ik].cutsq) continue;

        // 3-body zeta in bond order
        zeta_ij += zeta(&params[iparam_ijk], &params[iparam_ik],
                        rsq1, rsq2, delrj, delrk, i, xcn);

        // radical initialization: apply only to CC,CO,OC bonds
        if (params[iparam_ij].rad_flag > 0 &&
            params[iparam_ik].ielementgp == 1 &&
            params[iparam_ik].jelementgp == 1) {
          iparam_ki = elem2param[ktype][itype][itype];
          kcn=NCo[k];
          kconjug += rad_init(rsq2,&params[iparam_ki],i,kradtot,kcn);

        }

        // torsion: i-j-k-l: apply to all C-C bonds

        if (params[iparam_ij].tor_flag != 0) {
          srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
          srmu = sqrt(1.0-srmu*srmu);

          if(srmu > 0.1) {
            for (ll = 0; ll < sht_lnum; ll++) { // ll is neighbor of jj
              l = sht_llist[ll];

              if(l==i || l==j || l==k) continue;

              ltype = map[type[l]];

              delrl[0] = x[l][0] - x[j][0];
              delrl[1] = x[l][1] - x[j][1];
              delrl[2] = x[l][2] - x[j][2];
              rsq3 = vec3_dot(delrl,delrl);
              iparam_jl = elem2param[jtype][ltype][ltype];

              if (rsq3 > params[iparam_jl].cutsq) continue;

              iparam_ikl = elem2param[itype][ktype][ltype];
              torindx = params[iparam_ij].tor_flag;
              bbtor += bbtor1(torindx, &params[iparam_ikl],&params[iparam_jl],
                       rsq1,rsq2,rsq3,delrj,delrk,delrl,srmu);
            }
          }
        }
      }

      zeta_ji = 0.0;
      lconjug = 0.0;

      for (ll = 0; ll < sht_lnum; ll++) {
        l = sht_llist[ll];
        if (l == i) continue;

        ltype = map[type[l]];
        iparam_jil = elem2param[jtype][itype][ltype];
        iparam_ijl = elem2param[itype][jtype][ltype];
        iparam_jl  = elem2param[jtype][ltype][ltype];
        iparam_lj  = elem2param[ltype][jtype][jtype];

        delrk[0] = x[l][0] - x[j][0];
        delrk[1] = x[l][1] - x[j][1];
        delrk[2] = x[l][2] - x[j][2];
        rsq2 = vec3_dot(delrk,delrk);

        delrl[0] = x[l][0] - x[j][0];
        delrl[1] = x[l][1] - x[j][1];
        delrl[2] = x[l][2] - x[j][2];
        rsq2 = vec3_dot(delrl,delrl);

        if (rsq2 > params[iparam_jl].cutsq) continue;

        vec3_scale(-1,delrj,delrl);     // ji_hat is -(ij_hat)

        zeta_ji += zeta(&params[iparam_jil], &params[iparam_jl]
                        , rsq1, rsq2, delrl, delrk, j, ycn);

        // radical initialization: apply only to CC,CO,OC bonds
        if(params[iparam_ji].rad_flag > 0
          && params[iparam_jl].ielementgp == 1
          && params[iparam_jl].jelementgp == 1) {
          iparam_lj = elem2param[ltype][jtype][jtype];
          lcn=NCo[l];
          lconjug += rad_init(rsq2,&params[iparam_lj],j,lradtot,lcn);
        }
      }

      force_zeta(&params[iparam_ij], &params[iparam_ji],
         rsq1, xcn, ycn, zeta_ij, zeta_ji, fpair,
         prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
         prefac_ji1, prefac_ji2, prefac_ji3, prefac_ji4, prefac_ji5,
         eflag, eng_tmp, iq, jq, i, j, nj, bbtor, kconjug, lconjug);

      evdwl = eng_tmp;
      selfp6p(&params[iparam_ij],&params[iparam_ji],rsq1,eng_tmp,fpair);

      evdwl += eng_tmp;
      f[i][0] += delrj[0]*fpair;
      f[i][1] += delrj[1]*fpair;
      f[i][2] += delrj[2]*fpair;
      f[j][0] -= delrj[0]*fpair;
      f[j][1] -= delrj[1]*fpair;
      f[j][2] -= delrj[2]*fpair;

      if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,-fpair,-delrj[0],-delrj[1],-delrj[2]);

      // attractive term via loop over k (3-body forces: i-j-k)
      zet_addi=0;
      zet_addj=0;

      for (kk = 0; kk < sht_jnum; kk++) {
        k = sht_jlist[kk];
        if (j == k) continue;
        sht_mlist = sht_first[k];
        sht_mnum = sht_num[k];

        ktype = map[type[k]];
        iparam_ijk = elem2param[itype][jtype][ktype];
        iparam_ikj = elem2param[itype][ktype][jtype];
        iparam_jik = elem2param[jtype][itype][ktype];
        iparam_ik  = elem2param[itype][ktype][ktype];
        delrk[0] = x[k][0] - xtmp;
        delrk[1] = x[k][1] - ytmp;
        delrk[2] = x[k][2] - ztmp;
        rsq2 = vec3_dot(delrk,delrk);
        if (rsq2 > params[iparam_ik].cutsq) continue;

        // BO-dependent 3-body E & F
        attractive(&params[iparam_ijk], &params[iparam_jik],&params[iparam_ikj],
                prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
                rsq1,rsq2,delrj,delrk,fi,fj,fk,i,xcn);

        ep6p_ij = ep6p(&params[iparam_ijk],&params[iparam_ikj],rsq1,rsq2,delrj,delrk,zet_addi);
        fp6p(&params[iparam_ijk],&params[iparam_ikj],rsq1,rsq2,delrj,delrk,fip6p,fjp6p,fkp6p);

        // Sums up i-j-k forces: LP contribution
        for (im = 0; im < 3; im++) {
          fi[im] += fip6p[im];
          fj[im] += fjp6p[im];
          fk[im] += fkp6p[im];
        }

        // Sums up i-j-k forces: Tallies into global force vector
        for (im = 0; im < 3; im++) {
          f[i][im] += fi[im];
          f[j][im] += fj[im];
          f[k][im] += fk[im];
        }

        // torsion and radical: apply to all C-C bonds
        if (params[iparam_ijk].tor_flag != 0 && fabs(ptorr)>1.0e-8) {
          srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
          srmu = sqrt(1.0-srmu*srmu);

          if(srmu > 0.1) {
            for (ll = 0; ll < sht_lnum; ll++) { // ll is neighbor of jj
              l = sht_llist[ll];
              if (l==i||l==j||l==k) continue;

              ltype = map[type[l]];

              delrl[0] = x[l][0] - x[j][0];
              delrl[1] = x[l][1] - x[j][1];
              delrl[2] = x[l][2] - x[j][2];
              rsq3 = vec3_dot(delrl,delrl);

              iparam_jl = elem2param[jtype][ltype][ltype];
              if (rsq3 > params[iparam_jl].cutsq) continue;
              iparam_ikl = elem2param[itype][ktype][ltype];
              torindx = params[iparam_ij].tor_flag;
              tor_force(torindx, &params[iparam_ikl], &params[iparam_jl],srmu,
                          rsq1,rsq2,rsq3,delrj,delrk,delrl);

              for (im = 0; im < 3; im++) {
                f[i][im] += fi_tor[im];
                f[j][im] += fj_tor[im];
                f[k][im] += fk_tor[im];
                f[l][im] += fl_tor[im];
              }
            }
          }
        }

        if( params[iparam_ijk].rad_flag>=1 &&
          params[iparam_ijk].ielementgp==1 &&
          params[iparam_ijk].kelementgp==1) {
          iparam_ki = elem2param[ktype][itype][itype];
          kcn=NCo[k];
          double rik=sqrt(rsq2);
          kradtot = -comb_fc(rik,&params[iparam_ki])*params[iparam_ki].pcross+kcn;

          rad_forceik(&params[iparam_ki],rsq2,delrk,kconjug,kradtot);

          for (im = 0; im < 3; im++) {
            f[i][im] += fi_rad[im];
            f[k][im] += fk_rad[im];
          }

          if (fabs(radtmp) > 1.0e-12) {
            for (mm = 0; mm < sht_mnum; mm++) { // mm is neighbor of kk
               m = sht_mlist[mm];
               if (m == k) continue;

               mtype = map[type[m]];

              delrm[0] = x[m][0] - x[k][0];
              delrm[1] = x[m][1] - x[k][1];
              delrm[2] = x[m][2] - x[k][2];
              rsq3 = vec3_dot(delrm,delrm);

              iparam_km = elem2param[ktype][mtype][mtype];
              iparam_ki = elem2param[ktype][itype][itype];

              if (rsq3 > params[iparam_km].cutsq) continue;

              rad_force(&params[iparam_km],rsq3,delrm,radtmp);

              for (im = 0; im < 3; im++) {
                f[k][im] += fj_rad[im];
                f[m][im] += fk_rad[im];
              }
            }
          }
        }

        if (evflag)
          ev_tally(i,j,nlocal,newton_pair,ep6p_ij,0.0,0.0,0.0,0.0,0.0);
        if (vflag_atom)
          v_tally3(i,j,k,fj,fk,delrj,delrk);

      } // k-loop

      // attractive term via loop over l (3-body forces: j-i-l)
      for (ll = 0; ll < sht_lnum; ll++) {
        l = sht_llist[ll];
        if (l == i) continue;

        sht_plist = sht_first[l];
        sht_pnum = sht_num[l];

        ltype = map[type[l]];
        iparam_jil = elem2param[jtype][itype][ltype];
        iparam_jli = elem2param[jtype][ltype][itype];
        iparam_ijl = elem2param[itype][jtype][ltype];
        iparam_jl  = elem2param[jtype][ltype][ltype];
        delrk[0] = x[l][0] - x[j][0];
        delrk[1] = x[l][1] - x[j][1];
        delrk[2] = x[l][2] - x[j][2];

        rsq2 = vec3_dot(delrk,delrk);
        if (rsq2 > params[iparam_jl].cutsq) continue;
        vec3_scale(-1,delrj,delrl);

        attractive(&params[iparam_jil],&params[iparam_ijl],&params[iparam_jli],
                prefac_ji1,prefac_ji2,prefac_ji3,prefac_ji4,prefac_ji5,
                rsq1,rsq2,delrl,delrk,fj,fi,fl,j,ycn);

        // BO-independent 3-body j-i-l LP and BB correction and forces
        ep6p_ji = ep6p(&params[iparam_jil],&params[iparam_jli],rsq1,rsq2,delrl,delrk,zet_addj);
        fp6p(&params[iparam_jil],&params[iparam_jli],rsq1,rsq2,delrl,delrk,fjp6p,fip6p,flp6p);

        if (evflag)
          ev_tally(j,i,nlocal,newton_pair,ep6p_ji,0.0,0.0,0.0,0.0,0.0);

        // BO-dependent 3-body E & F
        for (im = 0; im < 3; im++) {
          fj[im] += fjp6p[im];
          fi[im] += fip6p[im];
          fl[im] += flp6p[im];
        }

        // Sums up j-i-l forces: Tallies into global force vector
        for (im = 0; im < 3; im++) {
          f[j][im] += fj[im];
          f[i][im] += fi[im];
          f[l][im] += fl[im];
        }

        // radical i-j-l-p: apply to all CC,CO,OC bonds
        if( params[iparam_jil].rad_flag >= 1 &&
          params[iparam_jil].ielementgp == 1 &&
          params[iparam_jil].kelementgp == 1 ) {
            iparam_lj = elem2param[ltype][jtype][jtype];
            lcn=NCo[l];
            double rjl=sqrt(rsq2);
            lradtot=-comb_fc(rjl,&params[iparam_lj])*params[iparam_lj].pcross +lcn;

            rad_forceik(&params[iparam_lj],rsq2,delrk,lconjug,lradtot);

            for (im = 0; im < 3; im++) {
              f[j][im] += fi_rad[im];
              f[l][im] += fk_rad[im];
            }

            if (fabs(radtmp)>1.0e-12) {
              for (pp = 0; pp < sht_pnum; pp++) {       // pp is neighbor of ll
                p = sht_plist[pp];
                if (p == l) continue;
                ptype = map[type[p]];

                delrp[0] = x[p][0] - x[l][0];
                delrp[1] = x[p][1] - x[l][1];
                delrp[2] = x[p][2] - x[l][2];
                rsq3 = vec3_dot(delrp,delrp);

                iparam_lp = elem2param[ltype][ptype][ptype];

                if (rsq3 > params[iparam_lp].cutsq) continue;

                vec3_scale(-1,delrj,delrj);
                rad_force(&params[iparam_lp],rsq3,delrp,radtmp);
                vec3_scale(-1,delrj,delrj);
                for (im = 0; im < 3; im++) {
                  f[l][im] += fj_rad[im];
                  f[p][im] += fk_rad[im];
                }
              }
            }
        }

        if (vflag_atom)
          v_tally3(j,i,l,fi,fl,delrl,delrk);
      }
    }
  }

  if (vflag_fdotr) virial_fdotr_compute();

}

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

void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
        double &fforce,int /*eflag*/, double &eng, double iq, double jq)
{
  double r,tmp_fc,tmp_fc_d,Di,Dj;
  double caj,vrcs,fvrcs;
  double LamDiLamDj,fcdA,rlm1,bigA;

  double romi = parami->addrep;
  double rrcs = parami->bigr + parami->bigd;
  double addr = parami->addrepr;

  r = sqrt(rsq);
  if (r > rrcs) return ;

  tmp_fc = comb_fc(r,parami);
  tmp_fc_d = comb_fc_d(r,parami);

  Di = parami->DU + pow(fabs(parami->bD*(parami->QU-iq)),parami->nD);
  Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-jq)),paramj->nD);

  bigA = parami->bigA;
  rlm1 = parami->lambda;

  fcdA = tmp_fc_d - tmp_fc * rlm1;
  LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);
  caj = bigA * LamDiLamDj;

  fforce = -caj * fcdA;

  // additional repulsion

  vrcs = 1.0; fvrcs = 0.0;
  if (romi != 0.0 && r < addr) {
    vrcs += romi * pow((1.0-r/addr),2.0);
    fvrcs = romi * 2.0 * (r/addr-1.0)/addr;
    fforce = fforce*vrcs - caj * tmp_fc * vrcs * fvrcs;
  }
  fforce /= r;

  // eng = repulsive energy
  eng = caj * tmp_fc * vrcs;
}

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

double PairComb3::zeta(Param *parami, Param *paramj, double rsqij,
        double rsqik, double *delrij, double *delrik, int /*i*/, double xcn)
{
  double rij,rik,costheta,arg,ex_delr,rlm3;

  rij = sqrt(rsqij);
  if (rij > parami->bigr+parami->bigd) return 0.0;
  rik = sqrt(rsqik);
  costheta = vec3_dot(delrij,delrik) / (rij*rik);

  rlm3 = parami->beta;
  arg = pow(rlm3*(rij-rik),int(parami->powermint));
  if (arg > 69.0776) ex_delr = 1.e30;
  else if (arg < -69.0776) ex_delr = 0.0;
  else ex_delr = exp(arg);

  return comb_fc(rik,paramj) * comb_gijk(costheta,parami,xcn) * ex_delr;
}

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

void PairComb3::selfp6p(Param *parami, Param *paramj, double rsq,
        double &eng, double &force)
{
  double r,comtti,comttj,fcj,fcj_d;

  r=sqrt(rsq);
  fcj=comb_fc(r,parami);
  fcj_d=comb_fc_d(r,parami);
  comtti = comttj = 0.0;

    double pilp0 = parami->p6p0;
    double pilp1 = parami->p6p1, pilp2 = parami->p6p2, pilp3 = parami->p6p3;
    double pilp4 = parami->p6p4, pilp5 = parami->p6p5, pilp6 = parami->p6p6;
    comtti = pilp0 + pilp1 + pilp2 + pilp3 + pilp4 + pilp5 + pilp6;

    double pjlp0 = paramj->p6p0;
    double pjlp1 = paramj->p6p1, pjlp2 = paramj->p6p2, pjlp3 = paramj->p6p3;
    double pjlp4 = paramj->p6p4, pjlp5 = paramj->p6p5, pjlp6 = paramj->p6p6;
    comttj = pjlp0 + pjlp1 + pjlp2 + pjlp3 + pjlp4 + pjlp5 + pjlp6;

  eng = 0.5 * fcj * (comtti + comttj);
  force += 0.5 * fcj_d * (comtti + comttj)/r;
}

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

double PairComb3::ep6p(Param *paramj, Param *paramk, double rsqij, double rsqik,
                     double *delrij, double *delrik , double &/*zet_add*/)
{
  double comtt;
  double pplp0 = paramj->p6p0;
  double pplp1 = paramj->p6p1, pplp2 = paramj->p6p2, pplp3 = paramj->p6p3;
  double pplp4 = paramj->p6p4, pplp5 = paramj->p6p5, pplp6 = paramj->p6p6;
  double rij,rik,costheta,lp0,lp1,lp2,lp3,lp4,lp5,lp6;
  double rmu,rmu2,rmu3,rmu4,rmu5,rmu6,fcj,fck;
    comtt=0.0;
    rij = sqrt(rsqij);
    rik = sqrt(rsqik);
    costheta = vec3_dot(delrij,delrik) / (rij*rik);
    fcj = comb_fc(rij,paramj);
    fck = comb_fc(rik,paramk);

    rmu = costheta;
    rmu2 = rmu *rmu; rmu3 = rmu2*rmu; rmu4 = rmu3*rmu;
    rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
    lp0 = pplp0;
    lp1 = pplp1*rmu;
    lp2 = pplp2*rmu2;
    lp3 = pplp3*rmu3;
    lp4 = pplp4*rmu4;
    lp5 = pplp5*rmu5;
    lp6 = pplp6*rmu6;
    comtt = lp0 + lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
    return 0.5 * fck * comtt *fcj;
  }

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

void PairComb3::fp6p(Param *paramij,Param *paramik, double rsqij, double rsqik,
                   double *delrij, double *delrik, double *drilp,
                   double *drjlp, double *drklp)
{
  double pplp0 = paramij->p6p0;
  double pplp1 = paramij->p6p1, pplp2 = paramij->p6p2, pplp3 = paramij->p6p3;
  double pplp4 = paramij->p6p4, pplp5 = paramij->p6p5, pplp6 = paramij->p6p6;
  double ffj1,ffj2,ffk1,ffk2;
  double rij,rik,costheta;
  double rmu,comtt,comtt_d,com4k,com5,com5k,fcj,fcj_d,fck,fck_d;
  double lp0,lp1,lp2,lp3,lp4,lp5,lp6;
  double lp1_d,lp2_d,lp3_d,lp4_d,lp5_d,lp6_d;
  double rmu2, rmu3, rmu4, rmu5, rmu6;

    ffj1 = 0.0, ffj2 = 0.0;
    ffk1 = 0.0, ffk2 = 0.0;

    rij = sqrt(rsqij); rik = sqrt(rsqik);
    costheta = vec3_dot(delrij,delrik) / (rij*rik);
    fcj = comb_fc(rij,paramij);
    fck = comb_fc(rik,paramik);
    fcj_d = comb_fc_d(rij,paramij);
    fck_d = comb_fc_d(rik,paramik);
    rmu = costheta;

      rmu2 = rmu *rmu; rmu3 = rmu2*rmu;
      rmu4 = rmu3*rmu; rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
      lp0 = pplp0;
      lp1 = pplp1*rmu;
      lp2 = pplp2*rmu2;
      lp3 = pplp3*rmu3;
      lp4 = pplp4*rmu4;
      lp5 = pplp5*rmu5;
      lp6 = pplp6*rmu6;
      lp1_d = pplp1;
      lp2_d = pplp2*2.0*rmu;
      lp3_d = pplp3*3.0*rmu2;
      lp4_d = pplp4*4.0*rmu3;
      lp5_d = pplp5*5.0*rmu4;
      lp6_d = pplp6*6.0*rmu5;
      comtt = lp0 + lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
      comtt_d = lp1_d + lp2_d + lp3_d + lp4_d + lp5_d + lp6_d;

      com4k = fcj * fck_d * comtt;
      com5  = fcj * fck * comtt_d;
      com5k = fck * comtt * fcj_d;

      ffj1 = 0.5*(-com5/(rij*rik));
      ffj2 = 0.5*(com5*rmu/rsqij-com5k/rij);
      ffk1 = ffj1;
      ffk2 = 0.5*(-com4k/rik+com5*rmu/rsqik);

  // j-atom
  vec3_scale(ffj1,delrik,drjlp);
  vec3_scaleadd(ffj2,delrij,drjlp,drjlp);

  // k-atom
  vec3_scale(ffk1,delrij,drklp);
  vec3_scaleadd(ffk2,delrik,drklp,drklp);

  // i-atom
  vec3_add(drjlp,drklp,drilp);
  vec3_scale(-1.0,drilp,drilp);
}

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

void PairComb3::force_zeta(Param *parami, Param *paramj, double rsq,
        double xcn, double ycn, double &zeta_ij, double &zeta_ji, double &fforce,
        double &prefac_ij1, double &prefac_ij2, double &prefac_ij3,
        double &prefac_ij4, double &prefac_ij5,
        double &prefac_ji1, double &prefac_ji2, double &prefac_ji3,
        double &prefac_ji4, double &prefac_ji5,
        int eflag, double &eng, double iq, double jq,
        int i, int j, int nj, double bbtor, double kconjug, double lconjug)
{
  double r,att_eng,att_force,bij;  // att_eng is -cbj
  double boij, dbij1, dbij2, dbij3, dbij4, dbij5;
  double boji, dbji1, dbji2, dbji3, dbji4, dbji5;
  double pradx, prady;
  r = sqrt(rsq);

  if (r > parami->bigr + parami->bigd) return;
  comb_fa(r, parami, paramj, iq, jq, att_eng, att_force);
  comb_bij_d(zeta_ij,parami,r,i,boij,dbij1,dbij2,dbij3,dbij4,dbij5,xcn);
  comb_bij_d(zeta_ji,paramj,r,j,boji,dbji1,dbji2,dbji3,dbji4,dbji5,ycn);
  bij = 0.5*(boij + boji);

  // radical energy

  if ( parami->rad_flag>0 ) {
    rad_calc( r, parami, paramj, kconjug, lconjug, i, j, xcn, ycn);
    bij +=  brad[0];
    pradx = brad[1]*att_eng;
    prady = brad[2]*att_eng;
    brad[3] = 1.0 * brad[3]*att_eng;
   }

  // torsion energy
  if ( parami->tor_flag!=0) {
     tor_calc( r, parami, paramj, kconjug, lconjug, i, j, xcn, ycn);
     bij += btor[0] * bbtor;
     ptorr =  att_eng * btor[0];
     pradx  += 1.0 *  btor[1] * bbtor * att_eng;
     prady  += 1.0 *  btor[2] * bbtor * att_eng;
     brad[3]+= 1.0 *  btor[3] * bbtor * att_eng;
  }

  fforce = 1.0*bij*att_force/r; // divide by r will done compute
  bbij[i][nj] = bij;

  prefac_ij1 = -0.5*att_eng*dbij1;      // prefac_ij1 = -pfij
  prefac_ij2 = -0.5*att_eng*dbij2;      // prefac_ij2 = -pfij1
  prefac_ij3 = -0.5*att_eng*dbij3;      // prefac_ij3 = -pfij2
  prefac_ij4 = -0.5*att_eng*dbij4;      // prefac_ij4 = -pfij3
  prefac_ij5 = -0.5*att_eng*dbij5;      // prefac_ij5 = -pfij4

  prefac_ji1 = -0.5*att_eng*dbji1;      // prefac_ji1 = -pfji
  prefac_ji2 = -0.5*att_eng*dbji2;      // prefac_ji2 = -pfji1
  prefac_ji3 = -0.5*att_eng*dbji3;      // prefac_ji3 = -pfji2
  prefac_ji4 = -0.5*att_eng*dbji4;      // prefac_ji4 = -pfji3
  prefac_ji5 = -0.5*att_eng*dbji5;      // prefac_ji5 = -pfji4

  // combines com6 & com7 below
  if ( parami->rad_flag>0 || parami->tor_flag!=0 ) {
    prefac_ij2-=pradx;
    prefac_ji2-=prady;
  }

  // eng = attraction energy
  if (eflag) eng = 1.0*bij*att_eng;
}

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

double PairComb3::comb_fc(double r, Param *param)
{
  double r_inn = param->bigr - param->bigd;
  double r_out = param->bigr + param->bigd;
  if (r <= r_inn) return 1.0;
  if (r >= r_out) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(r-r_inn)/(r_out-r_inn)));
}

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

double PairComb3::comb_fc_d(double r, Param *param)
{
  double r_inn = param->bigr - param->bigd;
  double r_out = param->bigr + param->bigd;
  if (r <= r_inn) return 0.0;
  if (r >= r_out) return 0.0;
  return -MY_PI2/(r_out-r_inn)*sin(MY_PI*(r-r_inn)/(r_out-r_inn));
}

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

double PairComb3::comb_fccc(double xcn)
{
  double cut1 = ccutoff[0];
  double cut2 = ccutoff[1];

  if (xcn <= cut1) return 1.0;
  if (xcn >= cut2) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}

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

double PairComb3::comb_fccc_d(double xcn)
{
  double cut1 = ccutoff[0];
  double cut2 = ccutoff[1];

  if (xcn <= cut1) return 0.0;
  if (xcn >= cut2) return 0.0;
  return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}

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

double PairComb3::comb_fcch(double xcn)
{
  double cut1 = ccutoff[2];
  double cut2 = ccutoff[3];

  if (xcn <= cut1) return 1.0;
  if (xcn >= cut2) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}

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

double PairComb3::comb_fcch_d(double xcn)
{
  double cut1 = ccutoff[2];
  double cut2 = ccutoff[3];

  if (xcn <= cut1) return 0.0;
  if (xcn >= cut2) return 0.0;
  return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}

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

double PairComb3::comb_fccch(double xcn)
{
  double cut1 = ccutoff[4];
  double cut2 = ccutoff[5];

  if (xcn <= cut1) return 1.0;
  if (xcn >= cut2) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}

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

double PairComb3::comb_fccch_d(double xcn)
{
  double cut1 = ccutoff[4];
  double cut2 = ccutoff[5];

  if (xcn <= cut1) return 0.0;
  if (xcn >= cut2) return 0.0;
  return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}

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

double PairComb3::comb_fcsw(double rsq)
{
  double r = sqrt(rsq);

  if (r <= chicut1) return 1.0;
  if (r >= chicut2) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(r-chicut1)/(chicut2-chicut1)));
}

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

double PairComb3::self(Param *param, double qi)
{
 double self_tmp, cmin, cmax, qmin, qmax;
 double s1=param->chi, s2=param->dj, s3=param->dk, s4=param->dl;

 self_tmp = 0.0;

 qmin = param->qmin;
 qmax = param->qmax;
 cmin = cmax = 100.0;

 self_tmp = qi*(s1+qi*(s2+qi*(s3+qi*s4)));

 if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4);
 if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4);

 return self_tmp;
}

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

void PairComb3::comb_fa(double r, Param *parami, Param *paramj, double iq,
        double jq, double &att_eng, double &att_force)
{
  double Bsi;
  double qi,qj,Di,Dj;
  double AlfDiAlfDj, YYBn, YYBj;
  double alfij1= parami->alpha1;
  double alfij2= parami->alpha2;
  double alfij3= parami->alpha3;
  double pbij1= parami->bigB1;
  double pbij2= parami->bigB2;
  double pbij3= parami->bigB3;
  if (r > parami->bigr + parami->bigd) Bsi = 0.0;

  qi = iq; qj = jq;
  Di = Dj = Bsi = 0.0;
  Di = parami->DU + pow(fabs(parami->bD*(parami->QU-qi)),parami->nD);
  Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-qj)),paramj->nD);
  YYBn = (parami->aB-fabs(pow(parami->bB*(qi-parami->Qo),10)));
  YYBj = (paramj->aB-fabs(pow(paramj->bB*(qj-paramj->Qo),10)));

  if (YYBn*YYBj > 0.0 ) {
    AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
    Bsi = (pbij1*exp(-alfij1*r)+pbij2*exp(-alfij2*r)+pbij3*exp(-alfij3*r))*
      sqrt(YYBn*YYBj)*AlfDiAlfDj;                               // Bsi is cbj

    att_eng = -Bsi * comb_fc(r,parami);
    att_force = -(Bsi*comb_fc_d(r,parami)-comb_fc(r,parami)*sqrt(YYBn*YYBj)*
        AlfDiAlfDj*(alfij1*pbij1*exp(-alfij1*r)+
        alfij2*pbij2*exp(-alfij2*r)+alfij3*pbij3*exp(-alfij3*r)));

  } else {
    att_eng = 0.0;
    att_force = 0.0;
  }

}

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

void PairComb3::comb_bij_d(double zet, Param *param, double r, int i,
        double &tbij, double &tbij1, double &tbij2,
        double &tbij3, double &tbij4, double &tbij5, double xcn)
{
  double pcorn,dpcorn,dxccij,dxchij,dxcoij;
  double zeta = zet;
  double zetang,tmp_tbij, pow_n;

  pcorn = dpcorn = dxccij = dxchij = dxcoij = 0.0;
  coord(param,r,i,pcorn,dpcorn,dxccij,dxchij,dxcoij,xcn);       // coordination term

  zetang=zeta;
  pow_n=param->powern;
  zeta = pow(zetang,pow_n)+pcorn;
  tmp_tbij=pow_n*pow(zetang,(pow_n-1.0));

  if ((1.0 + zeta) < 0.1 ){
    zeta=0.1-1.0;
    tbij = pow(1.0 + zeta, -0.5/pow_n);
    tbij1=0.0;
   }
   else if (zeta > param->c1) {
    tbij = pow(zeta,-0.5/pow_n);
    tbij1 = -0.5/pow_n*pow(zeta,(-0.5/pow_n-1.0));
   } else if (zeta > param->c2) {
    tbij = pow(zeta,-0.5/pow_n)-0.5/pow_n*pow(zeta,(-0.5/pow_n-1.0));
    tbij1 = -0.5/pow_n/zeta;
   } else if (fabs(zeta) < param->c4) {
    tbij = 1.0;
    tbij1 = 0.0;
   } else if (fabs(zeta) < param->c3) {
    tbij = 1.0 - zeta/(2.0*pow_n);
    tbij1 = -1/(2.0*pow_n);
   } else {
    tbij = pow(1.0 + zeta, -0.5/pow_n);
    tbij1 = -0.5/pow_n * pow(1.0 + zeta,(-1.0-0.5/pow_n));
   }

  tbij2 = tbij1 * dpcorn;
  tbij3 = tbij1 * dxccij;
  tbij4 = tbij1 * dxchij;
  tbij5 = tbij1 * dxcoij;
  tbij1 = tbij1 * tmp_tbij;

}

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

void PairComb3::coord(Param *param, double r, int i,
        double &pcorn, double &dpcorn, double &dxccij,
        double &dxchij, double &dxcoij, double xcn)
{
  int ixmin,iymin,izmin;
  double xcntot,xcccn,xchcn,xcocn;
  int tri_flag= param-> pcn_flag;
  int jele_gp= param->jelementgp;
  double pan = param->pcna;
  double pbn = param->pcnb;
  double pcn = param->pcnc;
  double pdn = param->pcnd;

  xcccn = xchcn = xcocn = 0.0;

  xcccn = xcctmp[i];
  xchcn = xchtmp[i];
  xcocn = xcotmp[i];
  xcntot = -comb_fc(r,param)*param->pcross + xcn;
  pcorn = dpcorn = dxccij = dxchij = dxcoij = 0.0;
  pcorn = 0.0; dpcorn = 0.0;

  if(xcntot  < 0.0) xcntot  = 0.0;

  if (tri_flag>0) {
    if(jele_gp==1) xcccn = xcccn-comb_fc(r,param)*param->pcross;
    if(jele_gp==2) xchcn = xchcn-comb_fc(r,param)*param->pcross;
    if(jele_gp==3) xcocn = xcocn-comb_fc(r,param)*param->pcross;
    if(xcccn < 0.0) xcccn = 0.0;
    if(xchcn < 0.0) xchcn = 0.0;
    if(xcocn < 0.0) xcocn = 0.0;
    if(xcccn > maxx) xcccn = maxx;
    if(xchcn > maxy) xchcn = maxy;
    if(xcocn > maxz) xcocn = maxz;

    double xcntritot=xcccn+xchcn+xcocn;

    if(xcntritot > maxxcn[tri_flag-1]) {
      pcorn  = vmaxxcn[tri_flag-1]+(xcntot-maxxcn[tri_flag-1])*dvmaxxcn[tri_flag-1];
      dxccij = dxchij = dxcoij = dvmaxxcn[tri_flag-1];
    }
    else {
      ixmin=int(xcccn+1.0e-12);
      iymin=int(xchcn+1.0e-12);
      izmin=int(xcocn+1.0e-12);
      if (fabs(float(ixmin)-xcccn)>1.0e-8 ||
          fabs(float(iymin)-xchcn)>1.0e-8 ||
          fabs(float(izmin)-xcocn)>1.0e-8) {
            cntri_int(tri_flag,xcccn,xchcn,xcocn,ixmin,iymin,izmin,
            pcorn,dxccij,dxchij,dxcoij,param);
      }
      else  {
        pcorn  = pcn_grid[tri_flag-1][ixmin][iymin][izmin];
        dxccij = pcn_gridx[tri_flag-1][ixmin][iymin][izmin];
        dxchij = pcn_gridy[tri_flag-1][ixmin][iymin][izmin];
        dxcoij = pcn_gridz[tri_flag-1][ixmin][iymin][izmin];
      }
    }
  } else {
    pcorn = pan*xcntot+pbn*exp(pcn*xcntot)+pdn;
    dpcorn = pan+pbn*pcn*exp(pcn*xcntot);
  }
}

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

void PairComb3::cntri_int(int tri_flag, double xval, double yval,
                double zval, int ixmin, int iymin, int izmin, double &vval,
                double &dvalx, double &dvaly, double &dvalz, Param * /*param*/)
{
  double x;
  vval = 0.0; dvalx = 0.0; dvaly = 0.0; dvalz = 0.0;
  if(ixmin >= maxx-1) { ixmin=maxx-1; }
  if(iymin >= maxy-1) { iymin=maxy-1; }
  if(izmin >= maxz-1) { izmin=maxz-1; }
  for (int j=0; j<64; j++) {
      x = pcn_cubs[tri_flag-1][ixmin][iymin][izmin][j]
          *pow(xval,iin3[j][0])*pow(yval,iin3[j][1])
          *pow(zval,iin3[j][2]);
    vval += x;
    if(xval>1.0e-8) {dvalx += x*iin3[j][0]/xval;}
    if(yval>1.0e-8) {dvaly += x*iin3[j][1]/yval;}
    if(zval>1.0e-8) {dvalz += x*iin3[j][2]/zval;}
  }
}

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

double PairComb3::comb_gijk(double costheta, Param *param, double nco_tmp)
{
  double rmu1 = costheta;
  double rmu2 = rmu1*rmu1;
  double rmu3 = rmu2*rmu1;
  double rmu4 = rmu3*rmu1;
  double rmu5 = rmu4*rmu1;
  double rmu6 = rmu5*rmu1;
  double co6 = param->pcos6*rmu6;
  double co5 = param->pcos5*rmu5;
  double co4 = param->pcos4*rmu4;
  double co3 = param->pcos3*rmu3;
  double co2 = param->pcos2*rmu2;
  double co1 = param->pcos1*rmu1;
  double co0 = param->pcos0;
  double pcross = param->pcross;
  double gmu;

  if (param->ang_flag==1) {
    double qtheta, gmu1, gmu2, rrmu, astep;
    int k;

    qtheta = comb_fccc(nco_tmp);
    astep = 2.0/ntab;
    rrmu = (rmu1+1.0)/astep;
    k = int(rrmu);
    gmu1 = co6+co5+co4+co3+co2+co1+co0;
    gmu2 = pang[k]+(pang[k+1]-pang[k])*(rrmu-k);
    gmu = gmu2+qtheta*(gmu1-gmu2);
    return gmu*pcross;

  } else if (param->ang_flag==2){
    double qtheta, gmu1, gmu2;
    double ch6 = ch_a[6]*rmu6;
    double ch5 = ch_a[5]*rmu5;
    double ch4 = ch_a[4]*rmu4;
    double ch3 = ch_a[3]*rmu3;
    double ch2 = ch_a[2]*rmu2;
    double ch1 = ch_a[1]*rmu1;
    double ch0 = ch_a[0];
    qtheta = comb_fccch(nco_tmp);
    gmu1 = co6+co5+co4+co3+co2+co1+co0;
    gmu2 = ch6+ch5+ch4+ch3+ch2+ch1+ch0;
    gmu = gmu2+qtheta*(gmu1-gmu2);
    return gmu*pcross;
  } else {
    gmu = co6+co5+co4+co3+co2+co1+co0;
    return gmu*pcross;
  }
}

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

void PairComb3::comb_gijk_d(double costheta, Param *param, double nco_tmp,
                double &gijk_d, double &com3jk)
{
  double rmu1 = costheta;
  double rmu2 = rmu1*rmu1;
  double rmu3 = rmu2*rmu1;
  double rmu4 = rmu3*rmu1;
  double rmu5 = rmu4*rmu1;
  double rmu6 = rmu5*rmu1;
  double co6 = param->pcos6; //*rmu5*6.0;
  double co5 = param->pcos5; //*rmu4*5.0;
  double co4 = param->pcos4; //*rmu3*4.0;
  double co3 = param->pcos3; //*rmu2*3.0;
  double co2 = param->pcos2; //*rmu1*2.0;
  double co1 = param->pcos1;
  double co0 = param->pcos0;
  double pcross = param->pcross;

  gijk_d = com3jk = 0.0;
  if (param->ang_flag==1) {
    double qtheta, dqtheta, gmu1, gmu2, dgmu1,dgmu2, rrmu, astep;
    int k;
    qtheta = comb_fccc(nco_tmp);
    dqtheta = comb_fccc_d(nco_tmp);

    astep = 2.0/ntab;
    rrmu = (rmu1+1.0)/astep;
    k = int(rrmu);

    gmu1 =rmu6*co6+rmu5*co5+rmu4*co4
         +rmu3*co3+rmu2*co2+rmu1*co1+co0;
    dgmu1 =6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
           +3.0*rmu2*co3+2.0*rmu1*co2+co1;
    gmu2 = pang[k]+(pang[k+1]-pang[k])*(rrmu-k);
    dgmu2 = dpang[k]+(dpang[k+1]-dpang[k])*(rrmu-k);
    gijk_d = pcross*(dgmu2+qtheta*(dgmu1-dgmu2));
    com3jk = dqtheta * (gmu1-gmu2);
  } else if(param->ang_flag==2) {
    double qtheta, dqtheta, gmu1, gmu2, dgmu1,dgmu2;
    double ch6 = ch_a[6];
    double ch5 = ch_a[5];
    double ch4 = ch_a[4];
    double ch3 = ch_a[3];
    double ch2 = ch_a[2];
    double ch1 = ch_a[1];
    double ch0 = ch_a[0];
    qtheta = comb_fccch(nco_tmp);
    dqtheta = comb_fccch_d(nco_tmp);

    gmu1 =rmu6*co6+rmu5*co5+rmu4*co4
         +rmu3*co3+rmu2*co2+rmu1*co1+co0;
    dgmu1 =6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
           +3.0*rmu2*co3+2.0*rmu1*co2+co1;
    gmu2 =rmu6*ch6+rmu5*ch5+rmu4*ch4
         +rmu3*ch3+rmu2*ch2+rmu1*ch1+ch0;
    dgmu2 =6.0*rmu5*ch6+5.0*rmu4*ch5+4.0*rmu3*ch4
           +3.0*rmu2*ch3+2.0*rmu1*ch2+ch1;
    gijk_d = pcross*(dgmu2+qtheta*(dgmu1-dgmu2));
    com3jk = dqtheta * (gmu1-gmu2);

  } else {
    gijk_d = pcross*(6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
                    +3.0*rmu2*co3+2.0*rmu1*co2+co1);
    com3jk = 0.0;
  }
}

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

void PairComb3::attractive(Param *parami, Param *paramj , Param *paramk, double prefac_ij1,
        double prefac_ij2, double prefac_ij3, double prefac_ij4,
        double prefac_ij5, double rsqij, double rsqik, double *delrij,
        double *delrik, double *fi, double *fj,double *fk, int /*i*/, double xcn)
{
  double rij_hat[3],rik_hat[3];
  double rij,rijinv,rik,rikinv;

  rij = sqrt(rsqij);
  rijinv = 1.0/rij;
  vec3_scale(rijinv,delrij,rij_hat);

  rik = sqrt(rsqik);
  rikinv = 1.0/rik;
  vec3_scale(rikinv,delrik,rik_hat);

  comb_zetaterm_d(prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
        rij_hat, rij,rik_hat, rik, fi, fj, fk, parami, paramj, paramk,xcn);

}

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

void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
        double prefac_ij3, double prefac_ij4, double prefac_ij5,
        double *rij_hat, double rij, double *rik_hat, double rik, double *dri,
        double *drj, double *drk, Param *parami, Param *paramj, Param *paramk, double xcn)
{
  double gijk,gijk_d,ex_delr,ex_delr_d,fc_k,cos_theta,tmp,rlm3;
  double dcosdri[3],dcosdrj[3],dcosdrk[3],dfc_i,dfc_k;
  double com6, com3j, com3k, com3jk;

  int mint = int(parami->powermint);
  double pcrossi = parami->pcross;
  double pcrossj = paramj->pcross;
  double pcrossk = paramk->pcross;
  int icontrol = parami->pcn_flag;

  dfc_i = comb_fc_d(rij,parami);
  fc_k = comb_fc(rik,paramk);
  dfc_k = comb_fc_d(rik,paramk);
  rlm3 = parami->beta;
  tmp = pow(rlm3*(rij-rik),mint);

  if (tmp > 69.0776) ex_delr = 1.e30;
  else if (tmp < -69.0776) ex_delr = 0.0;
  else ex_delr = exp(tmp);
  ex_delr *= pcrossi;

  cos_theta = vec3_dot(rij_hat,rik_hat);
  gijk = comb_gijk(cos_theta,parami,xcn);
  comb_gijk_d(cos_theta,parami,xcn,gijk_d,com3jk);
  costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);

  // com6 & com7
  if(icontrol > 0){
    if(parami->kelementgp==1) {com6 = prefac_ij3*pcrossk*dfc_k;}
    if(parami->kelementgp==2) {com6 = prefac_ij4*pcrossk*dfc_k;}
    if(parami->kelementgp==3) {com6 = prefac_ij5*pcrossk*dfc_k;}
    if(parami->rad_flag>=1 || parami->tor_flag!=0)
            {com6+=prefac_ij2*pcrossk*dfc_k;}
  } else {
    com6 = prefac_ij2*pcrossi*dfc_k;
  }

  if (parami->ang_flag==1 || parami->ang_flag==2) {
    com3j = com3jk*ex_delr*pcrossk*pcrossj*fc_k*dfc_i;
    com3k = com3jk*ex_delr*pcrossk*pcrossk*fc_k*dfc_k;
  } else {
    com3j = 0.0;
    com3k = 0.0;
  }

  ex_delr_d = mint*pow(rlm3,mint)*pow((rij-rik),(mint-1))*ex_delr; // com3
  vec3_scale(-dfc_k*gijk*ex_delr,rik_hat,dri);          // com1
  vec3_scaleadd(fc_k*gijk_d*ex_delr,dcosdri,dri,dri);   // com2
  vec3_scaleadd(fc_k*gijk*ex_delr_d,rik_hat,dri,dri);   // com3 cont'd
  vec3_scaleadd(-fc_k*gijk*ex_delr_d,rij_hat,dri,dri);  // com3 sums j
  vec3_scaleadd(-com3k,rik_hat,dri,dri);                // com3k
  vec3_scaleadd(-com3j,rij_hat,dri,dri);                // com3j
  vec3_scale(prefac_ij1,dri,dri);
  vec3_scaleadd(-com6,rik_hat,dri,dri);                 // com6

  vec3_scale(fc_k*gijk_d*ex_delr,dcosdrj,drj);          // com2
  vec3_scaleadd(fc_k*gijk*ex_delr_d,rij_hat,drj,drj);   // com3 cont'd
  vec3_scaleadd(com3j,rij_hat,drj,drj);                 // com3j
  vec3_scale(prefac_ij1,drj,drj);

  vec3_scale(dfc_k*gijk*ex_delr,rik_hat,drk);           // com1
  vec3_scaleadd(fc_k*gijk_d*ex_delr,dcosdrk,drk,drk);   // com2
  vec3_scaleadd(-fc_k*gijk*ex_delr_d,rik_hat,drk,drk);  // com3 cont'd
  vec3_scaleadd(com3k,rik_hat,drk,drk);                 // com3k
  vec3_scale(prefac_ij1,drk,drk);
  vec3_scaleadd(com6,rik_hat,drk,drk);                  // com6
}

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

void PairComb3::costheta_d(double *rij_hat, double rij, double *rik_hat,
        double rik, double *dri, double *drj, double *drk)
{
  double cos_theta = vec3_dot(rij_hat,rik_hat);

  vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
  vec3_scale(1.0/rij,drj,drj);
  vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
  vec3_scale(1.0/rik,drk,drk);
  vec3_add(drj,drk,dri);
  vec3_scale(-1.0,dri,dri);
}

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

void PairComb3::tables()

{
  int i,j,k,m, nntypes, ncoul,nnbuf, ncoul_lim, inty, itype, jtype;
  int iparam_i, iparam_ij, iparam_ji;
  double r,dra,drin,drbuf,rc,z,zr,zrc,ea,eb,ea3,eb3,alf;
  double exp2er,exp2ersh,fafash,dfafash,F1,dF1,ddF1,E1,E2,E3,E4;
  double exp2ear,exp2ebr,exp2earsh,exp2ebrsh,fafbsh,dfafbsh;
  double afbshift, dafbshift, exp2ershift;

  int n = nelements;

  dra  = 0.001;
  drin = 0.100;
  drbuf = 0.100;
  nnbuf = int(drbuf/dra) +1;
  rc = cutmax;
  alf = 0.20;
  nmax = atom->nmax;

  nntypes = int((n+1)*n/2.0)+1;
  ncoul = int((rc-drin)/dra)+ nnbuf;
  ncoul_lim = int(ncoul * 1.20);

  // allocate arrays
  memory->create(intype,n,n,"pair:intype");
  memory->create(erpaw,ncoul_lim,3,"pair:erpaw");
  memory->create(fafb,ncoul_lim,nntypes,"pair:fafb");
  memory->create(dfafb,ncoul_lim,nntypes,"pair:dfafb");
  memory->create(ddfafb,ncoul_lim,nntypes,"pair:ddfafb");
  memory->create(phin,ncoul_lim,nntypes,"pair:phin");
  memory->create(dphin,ncoul_lim,nntypes,"pair:dphin");
  memory->create(afb,ncoul_lim,nntypes,"pair:afb");
  memory->create(dafb,ncoul_lim,nntypes,"pair:dafb");
  memory->create(vvdw,ncoul,nntypes,"pair:vvdw");
  memory->create(vdvdw,ncoul,nntypes,"pair:vdvdw");
  memory->create(dpl,nmax,3,"pair:dpl");
  memory->create(bbij,nmax,MAXNEIGH,"pair:bbij");
  memory->create(xcctmp,nmax,"pair:xcctmp");
  memory->create(xchtmp,nmax,"pair:xchtmp");
  memory->create(xcotmp,nmax,"pair:xcotmp");
  memory->create(NCo,nmax,"pair:NCo");
  memory->create(sht_num,nmax,"pair:sht_num");
  sht_first = (int **) memory->smalloc(nmax*sizeof(int *),
        "pair:sht_first");

  // set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2

  m = 0; k = n;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      if (j == i) {
        intype[i][j] = m;
        m += 1;
      } else if (j != i && j > i) {
        intype[i][j] = k;
        k += 1;
      } else if (j != i && j < i) {
        intype[i][j] = intype[j][i];
      }
    }
  }

  // default arrays to zero

  for (i = 0; i < ncoul; i ++) {
    for (j = 0; j < nntypes; j ++) {
      fafb[i][j]   = 0.0;
      dfafb[i][j]  = 0.0;
      ddfafb[i][j] = 0.0;
      phin[i][j]   = 0.0;
      dphin[i][j]  = 0.0;
      afb[i][j]    = 0.0;
      dafb[i][j]   = 0.0;
    }
  }

  // direct 1/r energy with Slater 1S orbital overlap

  for (i = 0; i < n; i++) {
    if (map[i+1] < 0) continue;
    r = drin - dra;
    itype = map[i+1];
    iparam_i = elem2param[itype][itype][itype];
    z = params[iparam_i].esm;
    exp2ershift = exp(-2.0*z*rc);
    afbshift = -exp2ershift*(z+1.0/rc);
    dafbshift = exp2ershift*(2.0*z*z+2.0*z/rc+1.0/(rc*rc));

    if (comm->me == 0 && screen)
      fprintf(screen,"  element[%d] = %-2s, z = %g\n",i+1,elements[map[i+1]],z);

    for (j = 0; j < ncoul; j++) {
      exp2er = exp(-2.0 * z * r);
      phin[j][i] = 1.0 - exp2er * (1.0 + 2.0 * z * r * (1.0 + z * r));
      dphin[j][i] = (4.0 * exp2er * z * z * z * r * r);
      afb[j][i] = -exp2er*(z+1.0/r)-afbshift-(r-rc)*dafbshift;
      dafb[j][i] = -(exp2er*(2.0*z*z+2.0*z/r+1.0/(r*r))-dafbshift);
      r += dra;
    }
  }

  for (i = 0; i < n; i ++) {
    if (map[i+1] < 0) continue;
    for (j = 0; j < n; j ++) {
      if (map[j+1] < 0) continue;
      r = drin - dra;
      if (j == i) {
        itype = map[i+1];
        inty = intype[itype][itype];
        iparam_i = elem2param[itype][itype][itype];
        z = params[iparam_i].esm;
        zrc = z * rc;
        exp2ersh = exp(-2.0 * zrc);
        fafash = -exp2ersh * (1.0 / rc +
                              z * (11.0/8.0 + 3.0/4.0*zrc + zrc*zrc/6.0));
        dfafash = exp2ersh * (1.0/(rc*rc) + 2.0*z/rc +
                              z*z*(2.0 + 7.0/6.0*zrc + zrc*zrc/3.0));
        for (k = 0; k < ncoul; k ++) {
          zr = z * r;
          exp2er = exp(-2.0*zr);
          F1 = -exp2er * (1.0 / r +
                          z * (11.0/8.0 + 3.0/4.0*zr + zr*zr/6.0));
          dF1 = exp2er * (1.0/(r*r) + 2.0*z/r +
                          z*z*(2.0 + 7.0/6.0*zr + zr*zr/3.0));
          ddF1 = -exp2er * (2.0/(r*r*r) + 4.0*z/(r*r) + 4.0*z*z/r +
                            z*z*z/3.0*(17.0/2.0 + 5.0*zr + 2.0*zr*zr));
          fafb[k][inty] = F1-fafash-(r-rc)*dfafash;
          dfafb[k][inty] = -(dF1 - dfafash);
          ddfafb[k][inty] = ddF1;
                  r += dra;
        }
      } else if (j != i) {
        itype = map[i+1];
        jtype = map[j+1];
        inty = intype[itype][jtype];
        iparam_ij = elem2param[itype][jtype][jtype];
        ea = params[iparam_ij].esm;
        ea3 = ea*ea*ea;
        iparam_ji = elem2param[jtype][itype][itype];
        eb = params[iparam_ji].esm;
        eb3 = eb*eb*eb;
        E1 = ea*eb3*eb/((ea+eb)*(ea+eb)*(ea-eb)*(ea-eb));
        E2 = eb*ea3*ea/((ea+eb)*(ea+eb)*(eb-ea)*(eb-ea));
        E3 = (3.0*ea*ea*eb3*eb-eb3*eb3) /
          ((ea+eb)*(ea+eb)*(ea+eb)*(ea-eb)*(ea-eb)*(ea-eb));
        E4 = (3.0*eb*eb*ea3*ea-ea3*ea3) /
          ((ea+eb)*(ea+eb)*(ea+eb)*(eb-ea)*(eb-ea)*(eb-ea));
        exp2earsh = exp(-2.0*ea*rc);
        exp2ebrsh = exp(-2.0*eb*rc);
        fafbsh = -exp2earsh*(E1 + E3/rc)-exp2ebrsh*(E2 + E4/rc);
        dfafbsh =
          exp2earsh*(2.0*ea*(E1+E3/rc)+E3/(rc*rc)) +
          exp2ebrsh*(2.0*eb*(E2+E4/rc)+E4/(rc*rc));
        for (k = 0; k < ncoul; k ++) {
          exp2ear = exp(-2.0*ea*r);
          exp2ebr = exp(-2.0*eb*r);
          fafb[k][inty] =
            - exp2ear*(E1+E3/r) - exp2ebr*(E2+E4/r)
            - fafbsh - (r-rc) * dfafbsh;
          dfafb[k][inty] = -(exp2ear*(2.0*ea*(E1+E3/r) + E3/(r*r))
                           + exp2ebr*(2.0*eb*(E2+E4/r) + E4/(r*r))- dfafbsh);
          ddfafb[k][inty] = -exp2ear*(4.0*ea*ea*(E1+E3/r)+4.0*ea*E3/(r*r)
                                +2.0*E3/(r*r*r))
                            -exp2ebr*(4.0*eb*eb*(E2+E4/r)+4.0*eb*E4/(r*r)
                                +2.0*E4/(r*r*r));
          r += dra;
        }
      }
    }
  }

  for (i = 0; i < ncoul_lim; i ++) {
    r = dra * (i-1) + drin;
    erpaw[i][0] = erfc(r*alf);
    erpaw[i][1] = exp(-r*r*alf*alf);
  }
  // end wolf summation

  // van der Waals
  int ii,jj;
  double **rvdw, *cc2, *cc3, *vrc, *rrc;
  double r6, r7, r12, r13, rf6, rf12, drf7, drf13;
  double drcc, temp6, temp7, temp12, temp13;
  double vsigt, vepst, vdwt, dvdwt;

  vrc = new double[13];
  rrc = new double[13];
  cc2 = new double[nntypes];
  cc3 = new double[nntypes];
  memory->create(rvdw,2,nntypes,"pair:rvdw");

  vrc[0] = rc;
  for (i=1; i<13; i++) {
    vrc[i] = vrc[i-1] * vrc[0];
  }

  // generate spline coefficients for CC, CH, HH vdw
  for (ii = 0; ii < n; ii ++) {
    for (jj = ii; jj < n; jj ++) {
      itype = ii;
      jtype = jj;
      inty = intype[itype][jtype];
      iparam_ij = elem2param[itype][jtype][jtype];

      // parameter check: eps > 0
      if(params[iparam_ij].vdwflag > 0) {

        if(params[iparam_ij].vdwflag==1){
          rvdw[0][inty] = params[iparam_ij].bigr + params[iparam_ij].bigd;
        }
        else {
          rvdw[0][inty] = params[iparam_ij].bigr - params[iparam_ij].bigd;
        }

        rvdw[1][inty] = params[iparam_ij].vsig * 0.950;

        // radius check: outer radius vs. sigma
        if (rvdw[0][inty] > rvdw[1][inty])
          error->all(FLERR,"Error in vdw spline: inner radius > outer radius");

        rrc[0] = rvdw[1][inty];

        for (i=1; i<13; i++)
          rrc[i] = rrc[i-1] * rrc[0];

        drcc = rrc[0] - rvdw[0][inty];
        temp6 = 1.0/rrc[5]-1.0/vrc[5]+6.0*(rrc[0]-vrc[0])/vrc[6];
        temp7 = 6.0*(1.0/vrc[6]-1.0/rrc[6]);
        temp12 = 1.0/rrc[11]-1.0/vrc[11]+(rrc[0]-vrc[0])*12.0/vrc[12];
        temp13 = 12.0*(1.0/vrc[12]-1.0/rrc[12]);

        vsigt = params[iparam_ij].vsig;
        vepst = params[iparam_ij].veps;
        vsigt = vsigt*vsigt*vsigt*vsigt*vsigt*vsigt;

        vdwt = vepst*(vsigt*vsigt*temp12-vsigt*temp6);
        dvdwt = vepst*(vsigt*vsigt*temp13-vsigt*temp7);
        cc2[inty] = (3.0/drcc*vdwt-dvdwt)/drcc;
        cc3[inty] = (vdwt/(drcc*drcc)-cc2[inty] )/drcc;
      }
    }
  }

  // generate vdw look-up table
  for (ii = 0; ii < n; ii ++) {
    for (jj = ii; jj < n; jj ++) {
      itype = ii;
      jtype = jj;
      inty = intype[itype][jtype];
      iparam_ij = elem2param[itype][jtype][jtype];
      r = drin;
      for (k = 0; k < ncoul; k ++) {
        r6 = r*r*r*r*r*r;
        r7 = r6 * r;
        rf6 = 1.0/r6-1.0/vrc[5]+(r-vrc[0])*6.0/vrc[6];
        drf7 = 6.0*(1.0/vrc[6]-1.0/r7);
        vsigt = params[iparam_ij].vsig;
        vepst = params[iparam_ij].veps;
        vsigt = vsigt*vsigt*vsigt*vsigt*vsigt*vsigt;

        if(params[iparam_ij].vdwflag>0) {
          if(r <= rvdw[0][inty]) {
            vvdw[k][inty] = 0.0;
            vdvdw[k][inty] = 0.0;
          }
          else if ( r > rvdw[0][inty] && r <= rvdw[1][inty]) {
            drcc = r-rvdw[0][inty];
            vvdw[k][inty] = drcc*drcc*(drcc*cc3[inty]+cc2[inty]);
            vdvdw[k][inty] = drcc*(3.0*drcc*cc3[inty]+2.0*cc2[inty]);
          } else {
            r12 = r6*r6;
            r13 = r6*r7;
            rf12 = 1.0/r12-1.0/vrc[11]+(r-vrc[0])*12.0/vrc[12];
            drf13= 12.0*(1.0/vrc[12]-1.0/r13);
            vvdw[k][inty] = vepst*(vsigt*vsigt*rf12-vsigt*rf6);
            vdvdw[k][inty] = vepst*(vsigt*vsigt*drf13-vsigt*drf7);
          }
        } else {
          vvdw[k][inty]=0.0;
          vdvdw[k][inty]=0.0;
        }
          r += dra;
      }
    }
  }

  delete [] vrc;
  delete [] rrc;
  delete [] cc2;
  delete [] cc3;
  memory->destroy(rvdw);
}

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

void PairComb3::potal_calc(double &calc1, double &calc2, double &calc3)
{
  double alf,rcoul,esucon;
  int m;

  rcoul = 0.0;
  for (m = 0; m < nparams; m++)
    if (params[m].lcut > rcoul) rcoul = params[m].lcut;

  alf = 0.20;
  esucon = force->qqr2e;

  calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/MY_PIS*
           exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
  calc3 = (erfc(rcoul*alf)/rcoul)*esucon;
  calc1 = -(alf/MY_PIS*esucon+calc3*0.5);
}

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

void PairComb3::tri_point(double rsq, int &mr1, int &mr2,
                 int &mr3, double &sr1, double &sr2, double &sr3)
{
  double r, rin, dr, dd, rr1, rridr, rridr2;

  rin = 0.1000; dr = 0.0010;
  r = sqrt(rsq);
  if (r < rin + 2.0*dr)    r = rin + 2.0*dr;
  if (r > cutmax - 2.0*dr) r = cutmax - 2.0*dr;
  rridr = (r-rin)/dr;

  mr1 = int(rridr) ;
  dd = rridr - float(mr1);
  if (dd > 0.5) mr1 += 1;

  rr1 = float(mr1)*dr;
  rridr = (r - rin - rr1)/dr;
  rridr2 = rridr * rridr;

  sr1 = (rridr2 - rridr) * 0.50;
  sr2 = 1.0 - rridr2;
  sr3 = (rridr2 + rridr) * 0.50;

  mr2 = mr1 + 1;
  mr3 = mr1 + 2;
}

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

void PairComb3::vdwaals(int inty, int mr1, int mr2, int mr3, double rsq,
                      double sr1, double sr2, double sr3,
                      double &eng, double &fforce)
{
  double r = sqrt(rsq);

  eng = 1.0*(sr1*vvdw[mr1-1][inty]+sr2*vvdw[mr2-1][inty]+sr3*vvdw[mr3-1][inty]);
  fforce = -1.0/r*(sr1*vdvdw[mr1-1][inty]+sr2*vdvdw[mr2-1][inty]+sr3*vdvdw[mr3-1][inty]);
}

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

void PairComb3::direct(Param *parami, Param *paramj, int mr1,
        int mr2, int mr3, double rsq, double sr1, double sr2, double sr3,
        double iq, double jq, double fac11, double fac11e,
        double &pot_tmp, double &for_tmp, int i, int j)
{
  double r,erfcc,fafbnl,potij,esucon;
  double r3,erfcd,dfafbnl,smf2,dvdrr,alf,alfdpi;
  double afbn,afbj,sme1n,sme1j,sme1,sme2,dafbn, dafbj,smf1n,smf1j;
  double curli = parami->curl;
  double curlj = paramj->curl;
  int inti = parami->ielement;
  int intj = paramj->ielement;
  int inty = intype[inti][intj];

  double curlij0 = parami->curl0;
  double curlji0 = paramj->curl0;
  double curlij1,curlji1,dcurlij,dcurlji;
  double fcp1j,xcoij,xcoji;
  int icurl, jcurl;
  int ielegp = parami->ielementgp;
  int jelegp = paramj->ielementgp;

  r = sqrt(rsq);
  r3 = r * rsq;
  alf = 0.20;
  alfdpi = 2.0*alf/MY_PIS;
  esucon = force->qqr2e;
  pot_tmp = for_tmp = 0.0;
  icurl=jcurl=0;

  if(ielegp==2 && curli>curlij0) {
    icurl=1;
    curlij1=curli;
  }

  if(jelegp==2 && curlj>curlji0) {
    jcurl=1;
    curlji1=curlj;
  }

  if(icurl==1 || jcurl ==1) {
    xcoij = xcotmp[i];
    xcoji = xcotmp[j];
    fcp1j = comb_fc_d(r,parami);

    if(icurl==1) {
      curli=curlij1+(curlij0-curlij1)*comb_fc_curl(xcoij,parami);
      dcurlij=fcp1j*(curlij0-curlij1)*comb_fc_curl_d(xcoij,parami);
    }

    if(jcurl==1) {
      curlj=curlji1+(curlji0-curlji1)*comb_fc_curl(xcoji,paramj);
      dcurlji=fcp1j*(curlji0-curlji1)*comb_fc_curl_d(xcoji,paramj);
    }
  }

  erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
  afbn = sr1*afb[mr1][inti] + sr2*afb[mr2][inti] + sr3*afb[mr3][inti];
  afbj = sr1*afb[mr1][intj] + sr2*afb[mr2][intj] + sr3*afb[mr3][intj];
  fafbnl= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
  potij = (erfcc/r * esucon - fac11e);

  sme1n = iq*curlj*(afbn-fafbnl)*esucon;
  sme1j = jq*curli*(afbj-fafbnl)*esucon;
  sme1 = sme1n + sme1j;
  sme2 = (potij + fafbnl * esucon) * iq * jq;
  pot_tmp = 1.0 * (sme1+sme2);

  // 1/r force (wrt r)

  erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
  dafbn = sr1*dafb[mr1][inti] + sr2*dafb[mr2][inti] + sr3*dafb[mr3][inti];
  dafbj = sr1*dafb[mr1][intj] + sr2*dafb[mr2][intj] + sr3*dafb[mr3][intj];
  dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];

  dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
  smf1n = iq * curlj * (dafbn-dfafbnl)*esucon/r;
  smf1j = jq * curli * (dafbj-dfafbnl)*esucon/r;

  if(jcurl==1 && ielegp == 3 && dcurlji != 0.0){
   smf1n += dcurlji*iq*(afbn-fafbnl)*esucon/r;
  }
  if(icurl==1 && jelegp == 3 && dcurlij != 0.0){
   smf1j += dcurlij*jq*(afbj-fafbnl)*esucon/r;
  }

  smf2 = dvdrr + dfafbnl * esucon/r;
  for_tmp =  1.0 * iq * jq * smf2 + smf1n + smf1j;
}

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

void PairComb3::field(Param *parami, Param *paramj, double rsq, double iq,
                double jq, double &eng_tmp,double &for_tmp)
{
  double r,r3,r4,r5,rc,rc2,rc3,rc4,rc5;
  double cmi1,cmi2,cmj1,cmj2,pcmi1,pcmi2;
  double rf3i,rcf3i,rf5i,rcf5i;
  double drf3i,drcf3i,drf5i,drcf5i;
  double rf3,rf5,drf4,drf6;
  double smpn,smpl,rfx1,rfx2;

  r = sqrt(rsq);
  r3 = r * r * r;
  r4 = r3 * r;
  r5 = r4 * r;
  rc = parami->lcut;
  rc2 = rc * rc;
  rc3 = rc*rc*rc;
  rc4 = rc3 * rc;
  rc5 = rc4 * rc;
  cmi1 = parami->cmn1;
  cmi2 = parami->cmn2;
  cmj1 = paramj->cmn1;
  cmj2 = paramj->cmn2;
  pcmi1 = parami->pcmn1;
  pcmi2 = parami->pcmn2;

  rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
  rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
  rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
  rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));

  drf3i = 3/r*rf3i-6*rsq*rf3i*rf3i;
  drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
  drf5i = 5/r*rf5i-10*r4*rf5i*rf5i;
  drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;

  rf3 = rf3i-rcf3i-(r-rc)*drcf3i;
  rf5 = rf5i-rcf5i-(r-rc)*drcf5i;
  drf4 = drf3i - drcf3i;
  drf6 = drf5i - drcf5i;

 // field correction energy
  smpn = jq*(cmi1*rf3+jq*cmi2*rf5);
  smpl = iq*(cmj1*rf3+iq*cmj2*rf5);
  eng_tmp = 1.0 * (smpn + smpl);

 // field correction force
  rfx1 = jq*(cmi1*drf4+jq*cmi2*drf6)/r;
  rfx2 = iq*(cmj1*drf4+iq*cmj2*drf6)/r;
  for_tmp -= 1.0 * (rfx1 + rfx2);
}

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

double PairComb3::rad_init(double rsq2,Param *param,int /*i*/,
                double &radtot, double cnconj)
{
  double r, fc1k, radcut;

  r = sqrt(rsq2);
  fc1k = comb_fc(r,param);
  radtot = -fc1k * param->pcross + cnconj;
  radcut = comb_fcch(radtot);
  return fc1k * param->pcross * radcut;
}

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

void PairComb3::rad_calc(double r, Param *parami, Param *paramj,
        double kconjug, double lconjug, int /*i*/, int /*j*/, double xcn, double ycn)
{
  int ixmin, iymin, izmin;
  int radindx;
  double xrad, yrad, zcon, vrad, pradx, prady, pradz;

  vrad = pradx = prady = pradz = 0.0;
  xrad = -comb_fc(r,parami)*parami->pcross + xcn;
  yrad = -comb_fc(r,paramj)*paramj->pcross + ycn;
  zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);

  if(xrad < 0.0) xrad = 0.0;
  if(yrad < 0.0) yrad = 0.0;
  if(zcon < 1.0) zcon = 1.0;
  if(xrad > maxxc) xrad = maxxc;
  if(yrad > maxyc) yrad = maxyc;
  if(zcon > maxconj) zcon = maxconj;
  ixmin = int(xrad+1.0e-12);
  iymin = int(yrad+1.0e-12);
  izmin = int(zcon+1.0e-12);
  radindx=parami->rad_flag-1;
  if (fabs(float(ixmin)-xrad)>1.0e-8 ||
      fabs(float(iymin)-yrad)>1.0e-8 ||
      fabs(float(izmin)-zcon)>1.0e-8) {
    rad_int(radindx,xrad,yrad,zcon,ixmin,iymin,izmin,
              vrad,pradx,prady,pradz);
  } else {
    vrad  = rad_grid[radindx][ixmin][iymin][izmin-1];
    pradx = rad_gridx[radindx][ixmin][iymin][izmin-1];
    prady = rad_gridy[radindx][ixmin][iymin][izmin-1];
    pradz = rad_gridz[radindx][ixmin][iymin][izmin-1];
  }

  brad[0] = vrad;
  brad[1] = pradx;
  brad[2] = prady;
  brad[3] = pradz;
}

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

void PairComb3::rad_int(int radindx,double xrad, double yrad, double zcon, int l,
                int m, int n, double &vrad, double &pradx, double &prady,
                double &pradz)
{
  int j;
  double x;
  vrad = pradx = prady = pradz = 0.0;
  if(l >= maxxc-1) { l=maxxc-1;}
  if(m >= maxyc-1) { m=maxyc-1; }
  if(n >= maxconj-1) { n=maxconj-1;}

  for (j=0; j<64; j++) {
    x = rad_spl[radindx][l][m][n-1][j] * pow(xrad,iin3[j][0])
          * pow(yrad,iin3[j][1]) * pow(zcon,iin3[j][2]);
    vrad  += x;
    if(xrad > 1.0e-8) pradx += x*iin3[j][0]/xrad;
    if(yrad > 1.0e-8) prady += x*iin3[j][1]/yrad;
    if(zcon > 1.0e-8) pradz += x*iin3[j][2]/zcon;
  }
}


/* ---------------------------------------------------------------------- */
void PairComb3::rad_forceik(Param *paramk, double rsq2, double *delrk,
        double conjug, double radtot)
{
  int nm;
  double  rik, fc1k, fcp1k;
  double pradk, ffkk2, fktmp[3];
  double radcut = comb_fcch(radtot);
  double dradcut = comb_fcch_d(radtot);

  for (nm=0; nm<3; nm++) {
    fi_rad[nm] =  fk_rad[nm] = 0.0;
  }
    radtmp =0.0;

  rik = sqrt(rsq2);

  fc1k = comb_fc(rik, paramk);
  fcp1k = comb_fc_d(rik,paramk);

  pradk = brad[3]*fcp1k*radcut*paramk->pcross*2.0*conjug;
  radtmp= brad[3]*fc1k*dradcut*paramk->pcross*2.0*conjug;

  ffkk2 = -pradk/rik;

  for (nm=0; nm<3; nm++) {
    fktmp[nm] = - ffkk2 * delrk[nm];
  }

  for (nm=0; nm<3; nm++) {
    fi_rad[nm] =  fktmp[nm];
    fk_rad[nm] = -fktmp[nm];
  }
}

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

void PairComb3::rad_force(Param *paramm, double rsq3,
        double *delrm, double dpradk)
{
  int nm;
  double rkm, fcp1m;
  double comkm, ffmm2, fkm[3];

  for (nm=0; nm<3; nm++) {
    fj_rad[nm] = fk_rad[nm] = 0.0;
    fkm[nm]=0.0;
  }

  rkm = sqrt(rsq3);

  fcp1m = comb_fc_d(rkm, paramm);

  comkm = dpradk * fcp1m *  paramm->pcross;
  ffmm2 = -comkm/rkm;

  for (nm=0; nm<3; nm++) {
    fkm[nm] = -ffmm2 * delrm[nm];
  }

  for (nm=0; nm<3; nm++) {
    fj_rad[nm] =  fkm[nm];
    fk_rad[nm] = -fkm[nm];
  }
}

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

double PairComb3::bbtor1(int torindx, Param *paramk, Param *paraml,
        double rsq1, double rsq2, double rsq3, double *delrj,
        double *delrk, double *delrl, double srmu)
{
  double rmul, rij, rik, rjl;

  rij = sqrt(rsq1);
  rik = sqrt(rsq2);
  rjl = sqrt(rsq3);

  vec3_scale(-1.0,delrl,delrl);
  rmul = vec3_dot(delrj,delrl)/(rij*rjl);
  vec3_scale(-1.0,delrl,delrl);
  rmul = sqrt(1.0-rmul*rmul);

  if(rmul > 0.1 ) {
    double fc1k, fc1l, TT1, TT2, rmut, btt, tork[3], torl[3];

    fc1k = comb_fc(rik,paramk);
    fc1l = comb_fc(rjl,paraml);

    TT1 = rik*rjl*rij*rij*srmu*rmul;
    tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
    torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
    tork[1] = delrk[2]*delrj[0] - delrk[0]*delrj[2];
    torl[1] = delrj[2]*delrl[0] - delrj[0]*delrl[2];
    tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
    torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
    TT2 = vec3_dot(tork,torl);
    rmut = pow((TT2/TT1),2);
    if(torindx>=1) {
      btt = 1.0 - rmut;
      return btt * fc1k * fc1l;
    }
    else {
      btt=paramk->ptork1-TT2/TT1;
      btt=paramk->ptork2*pow(btt,2);
      return btt * fc1k * fc1l;
    }

  } else {
    return 0.0;
  }
}

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

void PairComb3::tor_calc(double r, Param *parami, Param *paramj,
        double kconjug, double lconjug, int /*i*/, int /*j*/, double xcn, double ycn)
{
  int ixmin, iymin, izmin;
  double vtor, dtorx, dtory, dtorz;
  double xtor, ytor, zcon;
  int torindx;

  vtor = dtorx = dtory = dtorz = 0.0;
  torindx=parami->tor_flag;

  if(torindx<0){
    vtor=1.0;
    dtorx=0.0;
    dtory=0.0;
    dtorz=0.0;
  } else {
    xtor = -comb_fc(r, parami) * parami->pcross + xcn;
    ytor = -comb_fc(r, paramj) * paramj->pcross + ycn;
    zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
    if (xtor < 0.0) xtor = 0.0;
    if (ytor < 0.0) ytor = 0.0;
    if (zcon < 1.0) zcon = 1.0;
    if (xtor > maxxc) xtor = maxxc;
    if (ytor > maxyc) ytor = maxyc;
    if (zcon > maxconj) zcon = maxconj;

    ixmin = int(xtor+1.0e-12);
    iymin = int(ytor+1.0e-12);
    izmin = int(zcon+1.0e-12);

    torindx=torindx-1;

    if (fabs(float(ixmin)-xtor)>1.0e-8 ||
      fabs(float(iymin)-ytor)>1.0e-8 ||
      fabs(float(izmin)-zcon)>1.0e-8) {
      tor_int(torindx,xtor,ytor,zcon,ixmin,iymin,izmin,
              vtor,dtorx,dtory,dtorz);
    } else {
      vtor  = tor_grid[torindx][ixmin][iymin][izmin-1];
      dtorx = tor_gridx[torindx][ixmin][iymin][izmin-1];
      dtory = tor_gridy[torindx][ixmin][iymin][izmin-1];
      dtorz = tor_gridz[torindx][ixmin][iymin][izmin-1];
    }
  }

  btor[0] = vtor;
  btor[1] = dtorx;
  btor[2] = dtory;
  btor[3] = dtorz;
}

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

void PairComb3::tor_int(int torindx,double xtor, double ytor, double zcon, int l,
        int m, int n, double &vtor, double &dtorx, double &dtory, double &dtorz)
{
  int j;
  double x;

  vtor = dtorx = dtory = dtorz = 0.0;
  if(l >= maxxc-1) { l=maxxc-1; }  //boundary condition changed
  if(m >= maxyc-1) { m=maxyc-1; }
  if(n >= maxconj-1) { n=maxconj-1; }

  for (j=0; j<64; j++) {
    x = tor_spl[torindx][l][m][n-1][j] * pow(xtor,iin3[j][0])
          * pow(ytor,iin3[j][1]) * pow(zcon,iin3[j][2]);
    vtor += x;

  if(xtor > 1.0e-8 ) dtorx += x*iin3[j][0]/xtor;
  if(ytor > 1.0e-8 ) dtory += x*iin3[j][1]/ytor;
  if(zcon > 1.0e-8 ) dtorz += x*iin3[j][2]/zcon;
  }
}

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

void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
        double srmu, double rsq1,double rsq2, double rsq3,
        double *delrj, double *delrk, double *delrl)
{
  int nm;
  double rmu, rmul, srmul, rij, rik, rjl;

  for (nm=0; nm<3; nm++) {
    fi_tor[nm] = fj_tor[nm] = fk_tor[nm] = fl_tor[nm] = 0.0;
  }

  rij = sqrt(rsq1);
  rik = sqrt(rsq2);
  rjl = sqrt(rsq3);

  rmu = vec3_dot(delrj,delrk)/(rij*rik);
  vec3_scale(-1.0,delrl,delrl);
  rmul = vec3_dot(delrj,delrl)/(rij*rjl);
  vec3_scale(-1.0,delrl,delrl);
  srmul = sqrt(1.0-rmul*rmul);
  if(acos(rmul) > MY_PI) srmul = -srmul;

  if(srmul > 0.1 ) {
    double fc1k, fcp1k, fc1l, fcp1l, srmul2, dt1dik, dt1djl;
    double TT1, TT2, rmut, btt, tork[3], torl[3];
    double dt2dik[3], dt2djl[3], dt2dij[3], AA, AA2;
    double tfij[4], tfik[2], tfjl[2], tjx[3], tjy[3], tjz[3];
    double tkx[2], tky[2], tkz[2], tlx[2], tly[2], tlz[2];

    fc1k  = comb_fc(rik,paramk);
    fcp1k = comb_fc_d(rik,paramk);
    fc1l  = comb_fc(rjl,paraml);
    fcp1l = comb_fc_d(rjl,paraml);
    srmul2 = pow(srmul,2);

    TT1 = rik*rjl*rij*rij*srmu*srmul;
    dt1dik = -rmu/pow(srmu,2);
    dt1djl = -rmul/srmul2;
    tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
    torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
    tork[1] = delrk[2]*delrj[0] - delrk[0]*delrj[2];
    torl[1] = delrj[2]*delrl[0] - delrj[0]*delrl[2];
    tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
    torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
    TT2 = vec3_dot(tork,torl);

    dt2dik[0] = -delrj[1]*torl[2] + delrj[2]*torl[1];
    dt2dik[1] = -delrj[2]*torl[0] + delrj[0]*torl[2];
    dt2dik[2] = -delrj[0]*torl[1] + delrj[1]*torl[0];
    dt2djl[0] =  delrj[1]*tork[2] - delrj[2]*tork[1];
    dt2djl[1] =  delrj[2]*tork[0] - delrj[0]*tork[2];
    dt2djl[2] =  delrj[0]*tork[1] - delrj[1]*tork[0];
    dt2dij[0] = -delrk[2]*torl[1] + delrl[2]*tork[1]
               + delrk[1]*torl[2] - delrl[1]*tork[2];
    dt2dij[1] = -delrk[0]*torl[2] + delrl[0]*tork[2]
               + delrk[2]*torl[0] - delrl[2]*tork[0];
    dt2dij[2] = -delrk[1]*torl[0] + delrl[1]*tork[0]
               + delrk[0]*torl[1] - delrl[0]*tork[1];

    rmut = TT2/TT1;

    if(torindx>=1) {
        btt = 1.0 - pow(rmut,2);
        AA = -2.0 * ptorr * rmut * fc1k * fc1l / TT1;
    }
    else {
        btt=paramk->ptork1-rmut;
        btt=paramk->ptork2*pow(btt,2);
        AA = -2.0 * ptorr * paramk->ptork2 *
          (paramk->ptork1-rmut) * fc1k * fc1l /TT1;
   }

    AA2 = AA * TT2;
    tfij[0] = -(dt1dik*AA2)/rij/rik;
    tfij[1] = AA2/rij/rij - dt1dik*AA2*rmu/rij/rij;
    tfij[2] = -dt1djl*AA2/rij/rjl;
    tfij[3] = AA2/rij/rij - dt1djl*AA2*rmul/rij/rij;
    tfik[0] = tfij[0];
    tfik[1] = (AA2/rik - btt*ptorr*fc1l*fcp1k)/rik -
            dt1dik*AA2*rmu/rik/rik;
    tfjl[0] = tfij[2];
    tfjl[1] = (AA2/rjl - btt*ptorr*fc1k*fcp1l)/rjl -
            dt1djl*AA2*rmul/rjl/rjl;

    tjx[0] = tfij[0]*delrk[0] - tfij[1]*delrj[0];
    tjy[0] = tfij[0]*delrk[1] - tfij[1]*delrj[1];
    tjz[0] = tfij[0]*delrk[2] - tfij[1]*delrj[2];
    tjx[1] = -tfij[2]*delrl[0] - tfij[3]*delrj[0];
    tjy[1] = -tfij[2]*delrl[1] - tfij[3]*delrj[1];
    tjz[1] = -tfij[2]*delrl[2] - tfij[3]*delrj[2];
    tjx[2] = -dt2dij[0] * AA;
    tjy[2] = -dt2dij[1] * AA;
    tjz[2] = -dt2dij[2] * AA;

    tkx[0] = tfik[0]*delrj[0] - tfik[1]*delrk[0];
    tky[0] = tfik[0]*delrj[1] - tfik[1]*delrk[1];
    tkz[0] = tfik[0]*delrj[2] - tfik[1]*delrk[2];
    tkx[1] = -dt2dik[0] * AA;
    tky[1] = -dt2dik[1] * AA;
    tkz[1] = -dt2dik[2] * AA;

    tlx[0] = -tfjl[0]*delrj[0] - tfjl[1]*delrl[0];
    tly[0] = -tfjl[0]*delrj[1] - tfjl[1]*delrl[1];
    tlz[0] = -tfjl[0]*delrj[2] - tfjl[1]*delrl[2];
    tlx[1] = -dt2djl[0] * AA;
    tly[1] = -dt2djl[1] * AA;
    tlz[1] = -dt2djl[2] * AA;

    fi_tor[0] = tjx[0]+tjx[1]+tjx[2]+tkx[0]+tkx[1];
    fi_tor[1] = tjy[0]+tjy[1]+tjy[2]+tky[0]+tky[1];
    fi_tor[2] = tjz[0]+tjz[1]+tjz[2]+tkz[0]+tkz[1];

    fj_tor[0] = -tjx[0]-tjx[1]-tjx[2]+tlx[0]+tlx[1];
    fj_tor[1] = -tjy[0]-tjy[1]-tjy[2]+tly[0]+tly[1];
    fj_tor[2] = -tjz[0]-tjz[1]-tjz[2]+tlz[0]+tlz[1];

    fk_tor[0] = -tkx[0]-tkx[1];
    fk_tor[1] = -tky[0]-tky[1];
    fk_tor[2] = -tkz[0]-tkz[1];

    fl_tor[0] = -tlx[0]-tlx[1];
    fl_tor[1] = -tly[0]-tly[1];
    fl_tor[2] = -tlz[0]-tlz[1];

  }
}

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

double PairComb3::combqeq(double *qf_fix, int &igroup)
{
  int i,j,ii, jj,itype,jtype,jnum;
  int iparam_i,iparam_ji,iparam_ij;
  int *ilist,*jlist,*numneigh,**firstneigh;
  int mr1,mr2,mr3,inty,nj;
  double xtmp,ytmp,ztmp,rsq1,delrj[3];
  double iq,jq,fqi,fqij,fqji,sr1,sr2,sr3;
  double potal,fac11,fac11e;
  int sht_jnum,*sht_jlist;
  tagint itag, jtag;

  double **x = atom->x;
  double *q = atom->q;
  tagint *tag = atom->tag;
  int *type = atom->type;
  int inum = list->inum;
  int *mask = atom->mask;
  int groupbit = group->bitmask[igroup];

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

  qf = qf_fix;
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    if (mask[i] & groupbit) {
      qf[i] = 0.0;
      dpl[i][0] = dpl[i][1] = dpl[i][2] = 0.0;
    }
  }
  // communicating charge force to all nodes, first forward then reverse

  pack_flag = 1;
  comm->forward_comm_pair(this);

  // self energy correction term: potal

  potal_calc(potal,fac11,fac11e);

  // loop over full neighbor list of my atoms

  fqi = fqij = fqji = 0.0;

  for (ii = 0; ii < inum; ii ++) {
    i = ilist[ii];
    itag = tag[i];
    nj = 0;
    if (mask[i] & groupbit) {
      itype = map[type[i]];
      xtmp = x[i][0];
      ytmp = x[i][1];
      ztmp = x[i][2];
      iq = q[i];
      iparam_i = elem2param[itype][itype][itype];

      // charge force from self energy
      fqi = qfo_self(&params[iparam_i],iq);

      jlist = firstneigh[i];
      jnum = numneigh[i];

      sht_jlist = sht_first[i];
      sht_jnum = sht_num[i];

      // two-body interactions

      for (jj = 0; jj < jnum; jj++) {
        j = jlist[jj] & NEIGHMASK;

        jtag = tag[j];
        if (itag >= jtag) continue;

        jtype = map[type[j]];
        inty = intype[itype][jtype];
        jq = q[j];

        delrj[0] = xtmp - x[j][0];
        delrj[1] = ytmp - x[j][1];
        delrj[2] = ztmp - x[j][2];
        rsq1 = vec3_dot(delrj,delrj);

        iparam_ij = elem2param[itype][jtype][jtype];
        iparam_ji = elem2param[jtype][itype][itype];

        // long range q-dependent

        if (rsq1 > params[iparam_ij].lcutsq) continue;

        // polynomial three-point interpolation
        tri_point(rsq1,mr1,mr2,mr3,sr1,sr2,sr3);

        // 1/r charge forces
        qfo_direct(&params[iparam_ij],&params[iparam_ji],
                mr1,mr2,mr3,rsq1,sr1,sr2,sr3,fac11e,fqij,fqji,
                iq,jq,i,j);

        fqi += fqij;  qf[j] += fqji;

       // field correction to self energy and charge force
        qfo_field(&params[iparam_ij],&params[iparam_ji],rsq1,
                iq,jq,fqij,fqji);

        fqi += fqij;  qf[j] += fqji;

        // polarization field charge force
        if (pol_flag) {
          qfo_dipole(fac11,mr1,mr2,mr3,inty,rsq1,delrj,sr1,sr2,sr3,
                fqij,fqji,i,j);

          fqi += fqij;  qf[j] += fqji;
        }
      }

      for (jj = 0; jj < sht_jnum; jj++) {
        j = sht_jlist[jj];

        jtag = tag[j];
        if (itag >= jtag) continue;

        jtype = map[type[j]];
        inty = intype[itype][jtype];
        jq = q[j];

        delrj[0] = xtmp - x[j][0];
        delrj[1] = ytmp - x[j][1];
        delrj[2] = ztmp - x[j][2];
        rsq1 = vec3_dot(delrj,delrj);

        iparam_ij = elem2param[itype][jtype][jtype];
        iparam_ji = elem2param[jtype][itype][itype];

        if (rsq1 >= params[iparam_ij].cutsq) continue;
        nj ++;

        // charge force in Aij and Bij
         qfo_short(&params[iparam_ij],&params[iparam_ji],
                rsq1,iq,jq,fqij,fqji,i,j,nj);

        fqi += fqij;  qf[j] += fqji;
      }
      qf[i] += fqi;
    }
  }

  comm->reverse_comm_pair(this);

  // sum charge force on each node and return it

  double eneg = 0.0;
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    if (mask[i] & groupbit){
      eneg += qf[i];
          itag=tag[i];
    }
  }

  MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
  MPI_Bcast(&enegtot,1,MPI_DOUBLE,0,world);
  return enegtot;

}

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

double PairComb3::qfo_self(Param *param, double qi)
{
  double self_d,cmin,cmax,qmin,qmax;
  double s1 = param->chi;
  double s2 = param->dj;
  double s3 = param->dk;
  double s4 = param->dl;

  self_d = 0.0;

  qmin = param->qmin;
  qmax = param->qmax;
  cmin = cmax = 100.0;
  self_d = s1+qi*(2.0*s2+qi*(3.0*s3+qi*4.0*s4));

  if (qi < qmin) self_d += 4.0 * cmin * pow((qi-qmin),3);
  if (qi > qmax) self_d += 4.0 * cmax * pow((qi-qmax),3);

  return self_d;
}

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

void PairComb3::qfo_direct(Param *parami, Param *paramj, int mr1,
          int mr2, int mr3, double rsq, double sr1, double sr2,
          double sr3, double fac11e, double &fqij, double &fqji,
          double iq, double jq, int i, int j)
{
  double r, erfcc, fafbnl, vm, vmfafb, esucon;
  double afbn, afbj, sme1n, sme1j;
  double curli = parami->curl;
  double curlj = paramj->curl;
  int inti = parami->ielement;
  int intj = paramj->ielement;
  int inty = intype[inti][intj];

  double curlij0 = parami->curl0;
  double curlji0 = paramj->curl0;
  double curlij1,curlji1;
  int icurl, jcurl;
  int ielegp = parami->ielementgp;
  int jelegp = paramj->ielementgp;

  r = sqrt(rsq);
  esucon=force->qqr2e;

  icurl = jcurl = 0;
  if(ielegp==2 && curli>curlij0) {
    icurl=1;
    curlij1=curli;
  }

  if(jelegp==2 && curlj>curlji0) {
    jcurl=1;
    curlji1=curlj;
  }
  if(icurl==1 || jcurl ==1) {
    double xcoij= xcotmp[i];
    double xcoji= xcotmp[j];

    if(icurl==1) {
      curli=curlij1+(curlij0-curlij1)*comb_fc_curl(xcoij,parami);
    }
    if(jcurl==1) {
      curlj=curlji1+(curlji0-curlji1)*comb_fc_curl(xcoji,paramj);
    }
  }

  // 1/r force (wrt q)

  erfcc = sr1*erpaw[mr1][0]   + sr2*erpaw[mr2][0]   + sr3*erpaw[mr3][0];
  fafbnl= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
  afbn  = sr1*afb[mr1][inti]  + sr2*afb[mr2][inti]  + sr3*afb[mr3][inti];
  afbj  = sr1*afb[mr1][intj]  + sr2*afb[mr2][intj]  + sr3*afb[mr3][intj];
  vm = (erfcc/r * esucon - fac11e);
  vmfafb = vm + esucon * fafbnl;
  sme1n = curlj * (afbn - fafbnl) * esucon;
  sme1j = curli * (afbj - fafbnl) * esucon;
  fqij = 1.0 * (jq * vmfafb + sme1n);
  fqji = 1.0 * (iq * vmfafb + sme1j);
}

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

void PairComb3::qfo_field(Param *parami, Param *paramj, double rsq,
        double iq,double jq, double &fqij, double &fqji)
{
  double r,r3,r5,rc,rc2,rc3,rc4,rc5;
  double cmi1,cmi2,cmj1,cmj2,pcmi1,pcmi2;
  double rf3i,rcf3i,rf5i,rcf5i;
  double drcf3i,drcf5i,rf3,rf5;

  r  = sqrt(rsq);
  r3 = r * rsq;
  r5 = r3 * rsq;
  rc = parami->lcut;
  rc2=  rc*rc;
  rc3 = rc*rc*rc;
  rc4 = rc3 * rc;
  rc5 = rc4 * rc;
  cmi1 = parami->cmn1;
  cmi2 = parami->cmn2;
  cmj1 = paramj->cmn1;
  cmj2 = paramj->cmn2;
  pcmi1 = parami->pcmn1;
  pcmi2 = parami->pcmn2;

  rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
  rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
  rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
  rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));

  drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
  drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;

  rf3 = rf3i-rcf3i-(r-rc)*drcf3i;
  rf5 = rf5i-rcf5i-(r-rc)*drcf5i;

  // field correction charge force
  fqij = 1.0 * cmj1*rf3+2.0*iq*cmj2*rf5;
  fqji = 1.0 * cmi1*rf3+2.0*jq*cmi2*rf5;

}

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

void PairComb3::qfo_dipole(double fac11, int mr1, int mr2, int mr3,
        int inty, double rsq, double *delrj, double sr1, double sr2,
        double sr3, double &fqij, double &fqji, int i, int j)
{
  double erfcc, erfcd, dvdrr, dfafbnl, smf2;
  double r, r3, alfdpi, esucon;

  r = sqrt(rsq);
  r3 = r * rsq;
  alfdpi = 0.4/MY_PIS;
  esucon = force->qqr2e;

  erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
  erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
  dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
  dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
  smf2 = (dvdrr + dfafbnl*esucon)/r;

  fqij = dpl[i][0]*delrj[0] + dpl[i][1]*delrj[1] +dpl[i][2]*delrj[2];
  fqji = dpl[j][0]*delrj[0] + dpl[j][1]*delrj[1] +dpl[j][2]*delrj[2];
  fqij *= smf2;
  fqji *= smf2;

}

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

void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
        double iq, double jq, double &fqij, double &fqji,
        int i, int /*j*/, int nj)
{
  double r, tmp_fc;
  double Di, Dj, dDi, dDj, Bsi, Bsj, dBsi, dBsj;
  double QUchi, QOchi, QUchj, QOchj;
  double bij, caj, cbj, caqpn, caqpj, cbqpn, cbqpj;
  double LamDiLamDj, AlfDiAlfDj;
  double rlm1 = parami->lambda;
  double alfij1= parami->alpha1;
  double alfij2= parami->alpha2;
  double alfij3= parami->alpha3;
  double pbij1= parami->bigB1;
  double pbij2= parami->bigB2;
  double pbij3= parami->bigB3;

  caj = cbj = caqpn = caqpj = cbqpn = cbqpj = 0.0;
  r = sqrt(rsq);
  tmp_fc = comb_fc(r,parami);
  bij = bbij[i][nj];

  QUchi = (parami->QU - iq) * parami->bD;
  QUchj = (paramj->QU - jq) * paramj->bD;
  QOchi = (iq - parami->Qo) * parami->bB;
  QOchj = (jq - paramj->Qo) * paramj->bB;

  if (iq < parami->QL-0.2) {
    iq = parami->QL-0.2;
    Di = parami->DL;
    dDi = Bsi = dBsi = 0.0;
  } else if (iq > parami->QU+0.2) {
    iq = parami->QU+0.2;
    Di = parami->DU;
    dDi = Bsi = dBsi = 0.0;
  } else {
    Di = parami->DU + pow(QUchi,parami->nD);                            // YYDin
    dDi = -parami->nD * parami->bD * pow(QUchi,(parami->nD-1.0));       // YYDiqp
    Bsi = parami->aB - pow(QOchi,10);                                   // YYBsin
    dBsi = -parami->bB * 10.0 * pow(QOchi,9.0);                         // YYBsiqp
  }

  if (jq < paramj->QL-0.2) {
    jq = paramj->QL-0.2;
    Dj = paramj->DL;
    dDj = Bsj = dBsj = 0.0;
  } else if (jq > paramj->QU+0.2) {
    jq = paramj->QU+0.2;
    Dj = paramj->DU;
    dDj = Bsj = dBsj = 0.0;
  } else {
    Dj = paramj->DU + pow(QUchj,paramj->nD);                            // YYDij
    dDj = -paramj->nD * paramj->bD * pow(QUchj,(paramj->nD-1.0));       // YYDiqpj
    Bsj = paramj->aB - pow(QOchj,10);                                   // YYBsij
    dBsj = -paramj->bB * 10.0 * pow(QOchj,9.0);                         // YYBsiqpj
  }

  LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);
  caj = 0.5 * tmp_fc * parami->bigA * LamDiLamDj;

  if (Bsi*Bsj > 0.0) {
    AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
    cbj=-0.5*tmp_fc*bij*sqrt(Bsi*Bsj)*AlfDiAlfDj*
                (pbij1*exp(-alfij1*r)+pbij2*exp(-alfij2*r)+pbij3*exp(-alfij3*r));
    cbqpn = cbj * (parami->alfi * dDi + dBsi/Bsi);
    cbqpj = cbj * (paramj->alfi * dDj + dBsj/Bsj);
  } else {
    cbj = cbqpn = cbqpj = 0.0;
  }

  caqpn = caj * parami->lami * dDi;
  caqpj = caj * paramj->lami * dDj;

  fqij = 1.0 * (caqpn + cbqpn);
  fqji = 1.0 * (caqpj + cbqpj);

}

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

void PairComb3::dipole_init(Param *parami, Param *paramj, double fac11,
        double *delrj, double rsq, int mr1, int mr2, int mr3, double sr1,
        double sr2, double sr3, double iq, double jq, int i, int j)
{
  double erfcc, erfcd, dvdrr, dfafbnl, smf2, phinn, phinj, efn, efj;
  double r, r3, alfdpi, esucon;
  double rcd, rct, tmurn, tmurj, poln[3], polj[3], Qext[3];
  int nm;
  int inti = parami->ielement;
  int intj = paramj->ielement;
  int inty = intype[inti][intj];

  for(nm=0; nm<3; nm++) Qext[nm] = 0.0;

  r = sqrt(rsq);
  r3 = r * rsq;
  rcd = 1.0/(r3);
  rct = 3.0*rcd/rsq;
  alfdpi = 0.4/MY_PIS;
  esucon = force->qqr2e;

  erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
  erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
  dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
  dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
  smf2 = dvdrr/esucon + dfafbnl/r;
  phinn = sr1*phin[mr1][inti] + sr2*phin[mr2][inti] + sr3*phin[mr3][inti];
  phinj = sr1*phin[mr1][intj] + sr2*phin[mr2][intj] + sr3*phin[mr3][intj];
  efn = jq * smf2;
  efj = iq * smf2;

  tmurn = dpl[i][0]*delrj[0] + dpl[i][1]*delrj[1] + dpl[i][2]*delrj[2];
  tmurj = dpl[j][0]*delrj[0] + dpl[j][1]*delrj[1] + dpl[j][2]*delrj[2];

  for (nm=0; nm<3; nm++) {
    poln[nm] = (tmurj*delrj[nm]*rct - dpl[j][nm]*rcd)*phinj;
    polj[nm] = (tmurn*delrj[nm]*rct - dpl[i][nm]*rcd)*phinn;
  }

  for (nm=0; nm<3; nm++) {
    dpl[i][nm] += (Qext[nm]/esucon + delrj[nm]*efn + poln[nm])*parami->polz*0.50;
    dpl[j][nm] += (Qext[nm]/esucon - delrj[nm]*efj + polj[nm])*paramj->polz*0.50;
  }

}

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

double PairComb3::dipole_self(Param *parami, int i)
{
  double esucon = force->qqr2e;
  double apn = parami->polz;
  double selfdpV = 0.0;

  if (apn != 0.0) {
      selfdpV= (dpl[i][0]*dpl[i][0]+dpl[i][1]*dpl[i][1]+dpl[i][2]*dpl[i][2])
                *esucon/(2.0*apn); }
  return selfdpV;
}

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

void PairComb3::dipole_calc(Param *parami, Param *paramj, double fac11,
        double delx, double dely, double delz, double rsq,
        int mr1, int mr2, int mr3, double sr1, double sr2, double sr3,
        double iq, double jq, int i, int j, double &vionij,
        double &fvionij, double *ddprx)
{
  double erfcc, erfcd, dvdrr, dfafbnl, ef, phinn, phinj, efn, efj;
  double r, r3, alf, alfdpi, esucon, dphinn, dphinj, ddfafbnl;
  double def, defn, defj, tmun, tmuj, emuTmu, edqn, edqj, ddvdrr;
  double rcd, rct, tmurn, tmurj, tmumu, poln[3], polj[3], delr1[3];
  double demuTmu, ddpr, dcoef;
  int nm;
  int inti = parami->ielement;
  int intj = paramj->ielement;
  int inty = intype[inti][intj];

  r = sqrt(rsq);
  r3 = r * rsq;
  esucon = force->qqr2e;
  rcd = esucon/r3;
  rct = 3.0*rcd/rsq;
  alf = 0.2;
  alfdpi = 2.0*alf/MY_PIS;
  delr1[0] = delx;
  delr1[1] = dely;
  delr1[2] = delz;

  // generate energy & force information from tables
  erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
  erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
  dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
  ddvdrr = (2.0*erfcc/r3 + 2.0*alfdpi*erfcd*(1.0/rsq+alf*alf))*esucon;
  dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
  phinn = sr1*phin[mr1][inti] + sr2*phin[mr2][inti] + sr3*phin[mr3][inti];
  phinj = sr1*phin[mr1][intj] + sr2*phin[mr2][intj] + sr3*phin[mr3][intj];
  dphinn = sr1*dphin[mr1][inti] + sr2*dphin[mr2][inti] + sr3*dphin[mr3][inti];
  dphinj = sr1*dphin[mr1][intj] + sr2*dphin[mr2][intj] + sr3*dphin[mr3][intj];
  ddfafbnl= sr1*ddfafb[mr1][inty] + sr2*ddfafb[mr2][inty] + sr3*ddfafb[mr3][inty];
  ef = (dvdrr + dfafbnl * esucon)/r;
  efn =  jq * ef;
  efj = -iq * ef;
  def = (ddvdrr + ddfafbnl * esucon)/r;
  defn =  jq * def;
  defj = -iq * def;

  // dipole - dipole field tensor (Tij)
  tmurn = dpl[i][0]*delr1[0] + dpl[i][1]*delr1[1] + dpl[i][2]*delr1[2];
  tmurj = dpl[j][0]*delr1[0] + dpl[j][1]*delr1[1] + dpl[j][2]*delr1[2];
  tmumu = dpl[i][0]*dpl[j][0] + dpl[i][1]*dpl[j][1] + dpl[i][2]*dpl[j][2];

  for (nm=0; nm<3; nm++) {
    poln[nm] = (tmurj*delr1[nm]*rct - dpl[j][nm]*rcd);
    polj[nm] = (tmurn*delr1[nm]*rct - dpl[i][nm]*rcd);
  }
  tmun = dpl[j][0]*polj[0] + dpl[j][1]*polj[1] + dpl[j][2]*polj[2];
  tmuj = dpl[i][0]*poln[0] + dpl[i][1]*poln[1] + dpl[i][2]*poln[2];

  // dipole - dipole energy
  emuTmu = -0.5*(tmun*phinn+tmuj*phinj);

  // dipole - charge energy
  edqn = -0.5 * (tmurn * efn);
  edqj = -0.5 * (tmurj * efj);

  // overall dipole energy
  vionij = emuTmu + edqn + edqj;

  // dipole - dipole force
  demuTmu = (tmun*dphinn + tmuj*dphinj)/r;
  ddpr =  5.0*tmurn*tmurj/rsq - tmumu;
  dcoef = rct * (phinn+phinj);

  for (nm = 0; nm < 3; nm ++) {
    ddprx[nm] = dcoef * (ddpr*delr1[nm] - tmurn*dpl[j][nm] - tmurj*dpl[i][nm])
          + demuTmu * delr1[nm];
  }

  // dipole - charge force
  fvionij = -tmurn*defn - tmurj*defj;
}

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

double PairComb3::comb_fc_curl(double rocn, Param *param)
{
  double r_inn = param->curlcut1;
  double r_out = param->curlcut2;
  if (rocn <= r_inn) return 1.0;
  if (rocn >= r_out) return 0.0;
  return 0.5*(1.0 + cos(MY_PI*(rocn-r_inn)/(r_out-r_inn)));
}

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

double PairComb3::comb_fc_curl_d(double rocn, Param *param)
{
  double r_inn = param->curlcut1;
  double r_out = param->curlcut2;
  if (rocn <= r_inn) return 0.0;
  if (rocn >= r_out) return 0.0;
  return -MY_PI2/(r_out-r_inn)*sin(MY_PI*(rocn-r_inn)/(r_out-r_inn));
}

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

int PairComb3::heaviside(double rr)
{
  if (rr <= 0.0) return 0;
  else return 1;
}

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

double PairComb3::switching(double rr)
{
  if (rr <= 0.0) return 1.0;
  else if (rr >= 1.0) return 0.0;
  else return heaviside(-rr)+heaviside(rr)*heaviside(1.0-rr)
          * (1.0-(3.0-2.0*rr)*rr*rr);
}

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

double PairComb3::switching_d(double rr)
{
  if (rr <= 0.0) return 0.0;
  else if (rr >= 1.0) return 0.0;
  else return heaviside(rr)*heaviside(1.0-rr)
          * 6.0*rr*(rr-1.0);
}

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

int PairComb3::pack_forward_comm(int n, int *list, double *buf,
                                 int /*pbc_flag*/, int * /*pbc*/)
{
  int i,j,m;

  m = 0;
  if (pack_flag == 1) {
    for (i = 0; i < n; i ++) {
      j = list[i];
      buf[m++] = qf[j];
    }
  } else if (pack_flag == 2) {
    for (i = 0; i < n; i ++) {
      j = list[i];
      buf[m++] = NCo[j];
    }
  }
  return m;
}

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

void PairComb3::unpack_forward_comm(int n, int first, double *buf)
{
  int i,m,last;

  m = 0;
  last = first + n ;
  if (pack_flag == 1) {
    for (i = first; i < last; i++)
      qf[i] = buf[m++];
  } else if (pack_flag == 2) {
    for (i = first; i < last; i++)
      NCo[i] = buf[m++];
  }
}

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

int PairComb3::pack_reverse_comm(int n, int first, double *buf)
{
  int i,m,last;

  m = 0;
  last = first + n;
  if (pack_flag == 1) {
    for (i = first; i < last; i++)
      buf[m++] = qf[i];
  } else if (pack_flag == 2) {
    for (i = first; i < last; i++)
      buf[m++] = NCo[i];
  }
  return m;
}

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

void PairComb3::unpack_reverse_comm(int n, int *list, double *buf)
{
  int i,j,m;

  m = 0;
  if (pack_flag == 1) {
    for (i = 0; i < n; i++) {
      j = list[i];
      qf[j] += buf[m++];
    }
  } else if (pack_flag == 2) {
    for (i = 0; i < n; i++) {
      j = list[i];
      NCo[j] += buf[m++];
    }
  }
}

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

double PairComb3::memory_usage()
{
  double bytes = maxeatom * sizeof(double);
  bytes += maxvatom*6 * sizeof(double);
  bytes += nmax * sizeof(int);
  bytes += nmax * 8.0 * sizeof(double);
  bytes += 25000*2*sizeof(double);

  for (int i = 0; i < comm->nthreads; i++)
    bytes += ipage[i].size();

  return bytes;
}