diff --git a/src/OPT/pair_ilp_graphene_hbn_opt.cpp b/src/OPT/pair_ilp_graphene_hbn_opt.cpp
index cceeec79ca..a35ee6a95d 100644
--- a/src/OPT/pair_ilp_graphene_hbn_opt.cpp
+++ b/src/OPT/pair_ilp_graphene_hbn_opt.cpp
@@ -41,10 +41,9 @@
 #include "neigh_request.h"
 #include "neighbor.h"
 #include "potential_file_reader.h"
-
+#include "tokenizer.h"
 #include <cmath>
 #include <cstring>
-#include <map>
 
 using namespace LAMMPS_NS;
 using namespace InterLayer;
@@ -52,6 +51,7 @@ using namespace InterLayer;
 #define MAXLINE 1024
 #define DELTA 4
 #define PGDELTA 1
+
 static const char cite_ilp[] =
     "ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848\n"
     "@Article{Ouyang2018\n"
@@ -82,7 +82,6 @@ static const char cite_ilp_cur[] =
     " location = {St. Louis, Missouri},\n"
     " series = {SC'21},\n"
     "}\n\n";
-
 /* ---------------------------------------------------------------------- */
 
 PairILPGrapheneHBNOpt::PairILPGrapheneHBNOpt(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
@@ -105,13 +104,16 @@ PairILPGrapheneHBNOpt::PairILPGrapheneHBNOpt(LAMMPS *lmp) : PairILPGrapheneHBN(l
   params = nullptr;
   cutILPsq = nullptr;
 
+  layered_neigh = nullptr;
+  first_layered_neigh = nullptr;
+  num_intra = nullptr;
+  num_inter = nullptr;
+  num_vdw = nullptr;
+  inum_max = 0;
+  jnum_max = 0;
   nmax = 0;
   maxlocal = 0;
 
-  normal = nullptr;
-  dnormal = nullptr;
-  dnormdri = nullptr;
-
   // always compute energy offset
   offset_flag = 1;
 
@@ -123,9 +125,7 @@ PairILPGrapheneHBNOpt::PairILPGrapheneHBNOpt(LAMMPS *lmp) : PairILPGrapheneHBN(l
 
 PairILPGrapheneHBNOpt::~PairILPGrapheneHBNOpt()
 {
-  memory->destroy(normal);
-  memory->destroy(dnormal);
-  memory->destroy(dnormdri);
+  delete[] pvector;
 
   if (allocated) {
     memory->destroy(setflag);
@@ -133,6 +133,12 @@ PairILPGrapheneHBNOpt::~PairILPGrapheneHBNOpt()
     memory->destroy(offset);
   }
 
+  memory->destroy(layered_neigh);
+  memory->sfree(first_layered_neigh);
+  memory->destroy(num_intra);
+  memory->destroy(num_inter);
+  memory->destroy(num_vdw);
+
   memory->destroy(elem2param);
   memory->destroy(cutILPsq);
   memory->sfree(params);
@@ -150,248 +156,172 @@ void PairILPGrapheneHBNOpt::init_style()
     error->all(FLERR, "Pair style ilp/graphene/hbn requires atom attribute molecule");
 
   // need a full neighbor list, including neighbors of ghosts
+
   neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST);
 }
 
 /* ---------------------------------------------------------------------- */
 void PairILPGrapheneHBNOpt::compute(int eflag, int vflag)
 {
-  ev_init(eflag, vflag);
-
-  //refactor and optimize for this pair style
-  computeILP(eflag, vflag);
-  return;
-}
-
-/* ----------------------------------------------------------------------
-  refactor and optimize for this pair style
-------------------------------------------------------------------------- */
-
-void PairILPGrapheneHBNOpt::computeILP(int eflag, int vflag)
-{
-  pvector[0] = pvector[1] = 0.0;
-  int i, j, ii, jj, inum, jnum, itype, jtype, k, l, kk, ll;
-  tagint itag, jtag;
-  double xtmp, ytmp, ztmp, delx, dely, delz, erep, fpair;
-  double rsq, r, Rcut, r2inv, r6inv, r8inv, Tap, dTap, Vilp, TSvdw, TSvdw2inv, fsum;
+  int i, j, ii, jj, inum, jnum, itype, itype_map, jtype, k, kk;
+  double prodnorm1, fkcx, fkcy, fkcz;
+  double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
+  double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp;
+  double frho1, Erep, fsum, rdsq1;
   int *ilist, *jlist, *numneigh, **firstneigh;
+  int *ILP_neighs_i;
+  tagint itag, jtag;
+  evdwl = 0.0;
 
-  erep = 0.0;
   double **x = atom->x;
   double **f = atom->f;
-  int *type = atom->type;
   tagint *tag = atom->tag;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
-
-  //Vars for calc_normal
-  int id, ip, m;
-  double nn, nn2;
-  double pv12[3], pv31[3], pv23[3], n1[3], dni[3], dnn[3][3], vet[3][3], dpvdri[3][3];
-  double dn1[3][3][3], dpv12[3][3][3], dpv23[3][3][3], dpv31[3][3][3];
-  double normal[3], dnormal[3][3][3], dnormdri[3][3];
-  //Vars for calc_Frep
-  double prodnorm1, fkcx, fkcy, fkcz;
-  double rhosq1, exp0, exp1;
-  double frho1, Erep, rdsq1, fpair1;
   double dprodnorm1[3] = {0.0, 0.0, 0.0};
   double fp1[3] = {0.0, 0.0, 0.0};
   double fprod1[3] = {0.0, 0.0, 0.0};
-  double delkj[3] = {0.0, 0.0, 0.0};
+  double delki[3] = {0.0, 0.0, 0.0};
   double fk[3] = {0.0, 0.0, 0.0};
-  double fkk[3][3], ekk[3], vkk[3][6];
-
-  //more for overflow
-  int ilp_neigh[6];
+  double dproddni[3] = {0.0, 0.0, 0.0};
+  double cij;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  ev_init(eflag, vflag);
+  pvector[0] = pvector[1] = 0.0;
+
+  if (neighbor->ago == 0) update_internal_list();
+
+  if (eflag_global || eflag_atom) {
+    if (vflag_either) {
+      if (tap_flag) {
+        eval<1, 1, 1>();
+      } else {
+        eval<1, 1, 0>();
+      }
+    } else {
+      if (tap_flag) {
+        eval<1, 0, 1>();
+      } else {
+        eval<1, 0, 0>();
+      }
+    }
+  } else {
+    if (vflag_either) {
+      if (tap_flag) {
+        eval<0, 1, 1>();
+      } else {
+        eval<0, 1, 0>();
+      }
+    } else {
+      if (tap_flag) {
+        eval<0, 0, 1>();
+      } else {
+        eval<0, 0, 0>();
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+template <int EFLAG, int VFLAG_EITHER, int TAP_FLAG> void PairILPGrapheneHBNOpt::eval()
+{
+  constexpr int EVFLAG = EFLAG || VFLAG_EITHER;
+  int i, j, ii, jj, inum, jnum, itype, itype_map, jtype, k, kk;
+  double prodnorm1, fkcx, fkcy, fkcz;
+  double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
+  double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp;
+  double frho1, Erep, fsum, rdsq1;
+  int *ilist, *jlist, *numneigh, **firstneigh;
+  int *ILP_neighs_i;
+  tagint itag, jtag;
+  evdwl = 0.0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+  double dprodnorm1[3] = {0.0, 0.0, 0.0};
+  double fp1[3] = {0.0, 0.0, 0.0};
+  double fprod1[3] = {0.0, 0.0, 0.0};
+  double delki[3] = {0.0, 0.0, 0.0};
+  double fk[3] = {0.0, 0.0, 0.0};
+  double dproddni[3] = {0.0, 0.0, 0.0};
+  double cij;
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  // loop over neighbors of my atoms
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
     itype = type[i];
-    int itypem = map[itype];
+    itype_map = map[type[i]];
     itag = tag[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    int nilp = 0;
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
-      jtype = type[j];
-      int jtypem = map[jtype];
-      jtag = tag[j];
+    int *jlist_intra = first_layered_neigh[i];
+    int *jlist_inter = first_layered_neigh[i] + num_intra[i];
+    int jnum_intra = num_intra[i];
+    int jnum_inter = num_inter[i];
+    int jnum_vdw = num_vdw[i];
+    int ILP_neigh[3];
+    int ILP_nneigh = 0;
+    for (jj = 0; jj < jnum_intra; jj++) {
+      j = jlist_intra[jj];
 
-      // two-body interactions from full neighbor list, skip half of them
+      jtype = map[type[j]];
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx * delx + dely * dely + delz * delz;
+
+      if (rsq < cutILPsq[itype_map][jtype]) {
+        if (ILP_nneigh >= 3)
+          error->one(FLERR, "There are too many neighbors for calculating normals");
+        ILP_neigh[ILP_nneigh++] = j;
+      }
+    }    // loop over jj
+
+    dproddni[0] = 0.0;
+    dproddni[1] = 0.0;
+    dproddni[2] = 0.0;
+
+    double norm[3], dnormdxi[3][3], dnormdxk[3][3][3];
+    calc_single_normal(i, ILP_neigh, ILP_nneigh, norm, dnormdxi, dnormdxk);
+
+    for (jj = 0; jj < jnum_inter; jj++) {
+      j = jlist_inter[jj];
+      jtype = type[j];
+      jtag = tag[j];
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
       rsq = delx * delx + dely * dely + delz * delz;
 
-      if (rsq != 0 && rsq < cutILPsq[itypem][jtypem] && atom->molecule[i] == atom->molecule[j]) {
-        ilp_neigh[nilp] = j;
-        vet[nilp][0] = -delx;
-        vet[nilp][1] = -dely;
-        vet[nilp][2] = -delz;
-        nilp++;
-      }
-    }
-    //nilp;
-    for (id = 0; id < 3; id++) {
-      pv12[id] = 0.0;
-      pv31[id] = 0.0;
-      pv23[id] = 0.0;
-      n1[id] = 0.0;
-      dni[id] = 0.0;
-      normal[id] = 0.0;
-      for (ip = 0; ip < 3; ip++) {
-        //vet[ip][id] = 0.0;
-        dnn[ip][id] = 0.0;
-        dpvdri[ip][id] = 0.0;
-        dnormdri[ip][id] = 0.0;
-        for (m = 0; m < 3; m++) {
-          dpv12[m][ip][id] = 0.0;
-          dpv31[m][ip][id] = 0.0;
-          dpv23[m][ip][id] = 0.0;
-          dn1[m][ip][id] = 0.0;
-          dnormal[m][ip][id] = 0.0;
-        }
-      }
-    }
-
-    if (nilp <= 1) {
-      normal[0] = 0.0;
-      normal[1] = 0.0;
-      normal[2] = 1.0;
-      for (id = 0; id < 3; id++) {
-        for (ip = 0; ip < 3; ip++) {
-          dnormdri[id][ip] = 0.0;
-          for (m = 0; m < 3; m++) { dnormal[m][id][ip] = 0.0; }
-        }
-      }
-    } else if (nilp == 2) {
-      cross_deriv(pv12, dpv12, vet, 0, 1, 2);
-      // derivatives of pv12[0] to ri
-      dpvdri[0][0] = 0.0;
-      dpvdri[0][1] = vet[0][2] - vet[1][2];
-      dpvdri[0][2] = vet[1][1] - vet[0][1];
-      // derivatives of pv12[1] to ri
-      dpvdri[1][0] = vet[1][2] - vet[0][2];
-      dpvdri[1][1] = 0.0;
-      dpvdri[1][2] = vet[0][0] - vet[1][0];
-      // derivatives of pv12[2] to ri
-      dpvdri[2][0] = vet[0][1] - vet[1][1];
-      dpvdri[2][1] = vet[1][0] - vet[0][0];
-      dpvdri[2][2] = 0.0;
-
-      n1[0] = pv12[0];
-      n1[1] = pv12[1];
-      n1[2] = pv12[2];
-      // the magnitude of the normal vector
-      nn2 = n1[0] * n1[0] + n1[1] * n1[1] + n1[2] * n1[2];
-      nn = sqrt(nn2);
-      if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
-      // the unit normal vector
-      normal[0] = n1[0] / nn;
-      normal[1] = n1[1] / nn;
-      normal[2] = n1[2] / nn;
-      // derivatives of nn, dnn:3x1 vector
-      dni[0] = (n1[0] * dpvdri[0][0] + n1[1] * dpvdri[1][0] + n1[2] * dpvdri[2][0]) / nn;
-      dni[1] = (n1[0] * dpvdri[0][1] + n1[1] * dpvdri[1][1] + n1[2] * dpvdri[2][1]) / nn;
-      dni[2] = (n1[0] * dpvdri[0][2] + n1[1] * dpvdri[1][2] + n1[2] * dpvdri[2][2]) / nn;
-      // derivatives of unit vector ni respect to ri, the result is 3x3 matrix
-      for (id = 0; id < 3; id++) {
-        for (ip = 0; ip < 3; ip++) {
-          dnormdri[id][ip] = dpvdri[id][ip] / nn - n1[id] * dni[ip] / nn2;
-        }
-      }
-      // derivatives of non-normalized normal vector, dn1:3x3x3 array
-      for (m = 0; m < 3; m++) {
-        for (id = 0; id < 3; id++) {
-          for (ip = 0; ip < 3; ip++) { dn1[m][id][ip] = dpv12[m][id][ip]; }
-        }
-      }
-      calc_dnormal(dnormal, dn1, n1, nn, nn2);
-
-    } else if (nilp == 3) {
-      cross_deriv(pv12, dpv12, vet, 0, 1, 2);
-      cross_deriv(pv31, dpv31, vet, 2, 0, 1);
-      cross_deriv(pv23, dpv23, vet, 1, 2, 0);
-
-      n1[0] = (pv12[0] + pv31[0] + pv23[0]) / jnum;
-      n1[1] = (pv12[1] + pv31[1] + pv23[1]) / jnum;
-      n1[2] = (pv12[2] + pv31[2] + pv23[2]) / jnum;
-      // the magnitude of the normal vector
-      nn2 = n1[0] * n1[0] + n1[1] * n1[1] + n1[2] * n1[2];
-      nn = sqrt(nn2);
-      if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
-      // the unit normal vector
-      normal[0] = n1[0] / nn;
-      normal[1] = n1[1] / nn;
-      normal[2] = n1[2] / nn;
-
-      // for the central atoms, dnormdri is always zero
-      for (id = 0; id < 3; id++) {
-        for (ip = 0; ip < 3; ip++) { dnormdri[id][ip] = 0.0; }
-      }
-
-      // derivatives of non-normalized normal vector, dn1:3x3x3 array
-      for (m = 0; m < 3; m++) {
-        for (id = 0; id < 3; id++) {
-          for (ip = 0; ip < 3; ip++) {
-            dn1[m][id][ip] = (dpv12[m][id][ip] + dpv23[m][id][ip] + dpv31[m][id][ip]) / jnum;
-          }
-        }
-      }
-      calc_dnormal(dnormal, dn1, n1, nn, nn2);
-    } else {
-      error->one(FLERR, "There are too many neighbors for calculating normals");
-    }
-    for (kk = 0; kk < nilp; kk++) {
-      fkk[kk][0] = 0.0;
-      fkk[kk][1] = 0.0;
-      fkk[kk][2] = 0.0;
-      vkk[kk][0] = 0.0;
-      vkk[kk][1] = 0.0;
-      vkk[kk][2] = 0.0;
-      vkk[kk][3] = 0.0;
-      vkk[kk][4] = 0.0;
-      vkk[kk][5] = 0.0;
-      ekk[kk] = 0.0;
-    }
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
-
-      jtype = type[j];
-      jtag = tag[j];
-      delx = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
-      rsq = delx * delx + dely * dely + delz * delz;
-
-      // only include the interation between different layers
-      if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
+      // only include the interaction between different layers
+      if (rsq < cutsq[itype][jtype]) {
 
         int iparam_ij = elem2param[map[itype]][map[jtype]];
-
-        //Param& p = params[iparam_ij];
-        Param *p = &params[iparam_ij];
+        Param &p = params[iparam_ij];
 
         r = sqrt(rsq);
-        double rinv = 1.0 / r;
-
+        double r2inv = 1.0 / rsq;
+        double rinv = r * r2inv;
         // turn on/off taper function
-        if (tap_flag) {
+        if (TAP_FLAG) {
           Rcut = sqrt(cutsq[itype][jtype]);
           Tap = calc_Tap(r, Rcut);
           dTap = calc_dTap(r, Rcut);
@@ -400,290 +330,499 @@ void PairILPGrapheneHBNOpt::computeILP(int eflag, int vflag)
           dTap = 0.0;
         }
 
-        int vdwflag = 1;
-        if (itag > jtag) {
-          if ((itag + jtag) % 2 == 0) vdwflag = 0;
-        } else if (itag < jtag) {
-          if ((itag + jtag) % 2 == 1) vdwflag = 0;
-        } else {
-          if (x[j][2] < ztmp) vdwflag = 0;
-          if (x[j][2] == ztmp && x[j][1] < ytmp) vdwflag = 0;
-          if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) vdwflag = 0;
-        }
-
-        double fvdw = 0, evdw = 0;
-
-        if (vdwflag) {
-          r2inv = rinv * rinv;
-          r6inv = r2inv * r2inv * r2inv;
-          r8inv = r6inv * r2inv;
-
-          TSvdw = 1.0 + exp(-p->d * (r / p->seff - 1.0));
-          TSvdw2inv = 1 / (TSvdw * TSvdw);    //pow(TSvdw,-2.0);
-          double vvdw = -p->C6 * r6inv / TSvdw;
-
-          // derivatives
-          fpair = -6.0 * p->C6 * r8inv / TSvdw +
-              p->C6 * p->d / p->seff * (TSvdw - 1.0) * TSvdw2inv * r8inv * r;
-          fsum = fpair * Tap - vvdw * dTap * rinv;
-
-          fvdw = fsum;
-          evdw = vvdw * Tap;
-          pvector[0] += evdw;
-        }
-
         // Calculate the transverse distance
-        prodnorm1 = normal[0] * delx + normal[1] * dely + normal[2] * delz;
+        prodnorm1 = norm[0] * delx + norm[1] * dely + norm[2] * delz;
         rhosq1 = rsq - prodnorm1 * prodnorm1;    // rho_ij
-        rdsq1 = rhosq1 * p->delta2inv;           // (rho_ij/delta)^2
+        rdsq1 = rhosq1 * p.delta2inv;            // (rho_ij/delta)^2
 
         // store exponents
-        exp0 = exp(-p->lambda * (r - p->z0));
+        exp0 = exp(-p.lambda * (r - p.z0));
         exp1 = exp(-rdsq1);
 
-        frho1 = exp1 * p->C;
-        Erep = 0.5 * p->epsilon + frho1;
+        frho1 = exp1 * p.C;
+        Erep = 0.5 * p.epsilon + frho1;
+
         Vilp = exp0 * Erep;
 
         // derivatives
-        fpair = p->lambda * exp0 * rinv * Erep;
-        fpair1 = 2.0 * exp0 * frho1 * p->delta2inv;
+        fpair = p.lambda * exp0 * rinv * Erep;
+        fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
         fsum = fpair + fpair1;
         // derivatives of the product of rij and ni, the result is a vector
-        dprodnorm1[0] = dnormdri[0][0] * delx + dnormdri[1][0] * dely + dnormdri[2][0] * delz;
-        dprodnorm1[1] = dnormdri[0][1] * delx + dnormdri[1][1] * dely + dnormdri[2][1] * delz;
-        dprodnorm1[2] = dnormdri[0][2] * delx + dnormdri[1][2] * dely + dnormdri[2][2] * delz;
-        fp1[0] = prodnorm1 * normal[0] * fpair1;
-        fp1[1] = prodnorm1 * normal[1] * fpair1;
-        fp1[2] = prodnorm1 * normal[2] * fpair1;
-        fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
-        fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
-        fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
+
+        fp1[0] = prodnorm1 * norm[0] * fpair1;
+        fp1[1] = prodnorm1 * norm[1] * fpair1;
+        fp1[2] = prodnorm1 * norm[2] * fpair1;
 
         fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx * rinv;
         fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely * rinv;
         fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz * rinv;
 
-        //This should be no use because fkcx need a lot of variables
-        //fi + fj + sum(fk) = 0
-        //-sum(fk) = fi + fj = -fprod*Tap
-        //sum(fk) = fprod * Tap
-        //fj = -fi - fprod*Tap
-        double ftotx = fvdw * delx + fkcx;
-        double ftoty = fvdw * dely + fkcy;
-        double ftotz = fvdw * delz + fkcz;
-        f[i][0] += ftotx - fprod1[0] * Tap;
-        f[i][1] += ftoty - fprod1[1] * Tap;
-        f[i][2] += ftotz - fprod1[2] * Tap;
-        f[j][0] -= ftotx;
-        f[j][1] -= ftoty;
-        f[j][2] -= ftotz;
+        f[i][0] += fkcx;
+        f[i][1] += fkcy;
+        f[i][2] += fkcz;
+        f[j][0] -= fkcx;
+        f[j][1] -= fkcy;
+        f[j][2] -= fkcz;
 
-        // calculate the forces acted on the neighbors of atom i from atom j
-        //ILP_neighs_i = ilp_neigh;
-        for (kk = 0; kk < nilp; kk++) {
-          k = ilp_neigh[kk];
-          //if (k == i) continue;
-          // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
-          dprodnorm1[0] =
-              dnormal[kk][0][0] * delx + dnormal[kk][1][0] * dely + dnormal[kk][2][0] * delz;
-          dprodnorm1[1] =
-              dnormal[kk][0][1] * delx + dnormal[kk][1][1] * dely + dnormal[kk][2][1] * delz;
-          dprodnorm1[2] =
-              dnormal[kk][0][2] * delx + dnormal[kk][1][2] * dely + dnormal[kk][2][2] * delz;
-          fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
-          fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
-          fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
-          fkk[kk][0] += fk[0];
-          fkk[kk][1] += fk[1];
-          fkk[kk][2] += fk[2];
-          delkj[0] = x[i][0] + vet[kk][0] - x[j][0];
-          delkj[1] = x[i][1] + vet[kk][1] - x[j][1];
-          delkj[2] = x[i][2] + vet[kk][2] - x[j][2];
-          //if (evflag) ev_tally_xyz(k,j,nlocal,newton_pair,0.0,0.0,fk[0],fk[1],fk[2],delkj[0],delkj[1],delkj[2]);
-          if (vflag_atom) {
-            if (evflag)
-              ev_tally_buffer(k, j, ekk + kk, vkk[kk], eatom + j, vatom[j], nlocal, newton_pair,
-                              0.0, 0.0, fk[0], fk[1], fk[2], delkj[0], delkj[1], delkj[2]);
-          } else {
-            if (evflag)
-              ev_tally_buffer(k, j, ekk + kk, vkk[kk], eatom + j, NULL, nlocal, newton_pair, 0.0,
-                              0.0, fk[0], fk[1], fk[2], delkj[0], delkj[1], delkj[2]);
-          }
-        }
-        erep = Tap * Vilp;
-        if (eflag) pvector[1] += erep;    // = Tap*Vilp;
-        //if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,fkcx,fkcy,fkcz,delx,dely,delz);
-        if (evflag)
-          ev_tally_xyz(i, j, nlocal, newton_pair, erep + evdw, 0.0, ftotx, ftoty, ftotz, delx, dely,
+        cij = -prodnorm1 * fpair1 * Tap;
+        dproddni[0] += cij * delx;
+        dproddni[1] += cij * dely;
+        dproddni[2] += cij * delz;
+
+        if (EFLAG) pvector[1] += evdwl = Tap * Vilp;
+        if (EVFLAG)
+          ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
+
+        /* ----------------------------------------------------------------------
+	   van der Waals forces and energy
+	   ------------------------------------------------------------------------- */
+        if (jj >= jnum_vdw) continue;
+        double r6inv = r2inv * r2inv * r2inv;
+        double r8inv = r6inv * r2inv;
+
+        //double TSvdw = 1.0 + exp(-p.d * (r * p.seffinv - 1.0));
+        double TSvdw = 1.0 + exp(-p.d * (r * 1.0 / p.seff - 1.0));
+        double TSvdwinv = 1.0 / TSvdw;
+        double TSvdw2inv = TSvdwinv * TSvdwinv;    //pow(TSvdw, -2.0);
+        Vilp = -p.C6 * r6inv * TSvdwinv;
+
+        fpair = -6.0 * p.C6 * r8inv * TSvdwinv +
+            p.C6 * p.d * 1.0 / p.seff * (TSvdw - 1.0) * TSvdw2inv * r8inv * r;
+        //p.C6 * p.d * p.seffinv * (TSvdw - 1.0) * TSvdw2inv * r8inv * r;
+        fsum = fpair * Tap - Vilp * dTap * rinv;
+
+        double fvdwx = fsum * delx;
+        double fvdwy = fsum * dely;
+        double fvdwz = fsum * delz;
+
+        f[i][0] += fvdwx;
+        f[i][1] += fvdwy;
+        f[i][2] += fvdwz;
+        f[j][0] -= fvdwx;
+        f[j][1] -= fvdwy;
+        f[j][2] -= fvdwz;
+
+        if (EFLAG) pvector[0] += evdwl = Tap * Vilp;
+        if (EVFLAG)
+          ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fvdwx, fvdwy, fvdwz, delx, dely,
                        delz);
       }
     }    // loop over jj
-    for (kk = 0; kk < nilp; kk++) {
-      int k = ilp_neigh[kk];
-      f[k][0] += fkk[kk][0];
-      f[k][1] += fkk[kk][1];
-      f[k][2] += fkk[kk][2];
-      if (eflag_atom) eatom[k] += ekk[kk];
-      if (vflag_atom) {
-        vatom[k][0] += vkk[kk][0];
-        vatom[k][1] += vkk[kk][1];
-        vatom[k][2] += vkk[kk][2];
-        vatom[k][3] += vkk[kk][3];
-        vatom[k][4] += vkk[kk][4];
-        vatom[k][5] += vkk[kk][5];
+
+    for (kk = 0; kk < ILP_nneigh; kk++) {
+      k = ILP_neigh[kk];
+      if (k == i) continue;
+      // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
+      fk[0] = dnormdxk[0][0][kk] * dproddni[0] + dnormdxk[1][0][kk] * dproddni[1] +
+          dnormdxk[2][0][kk] * dproddni[2];
+      fk[1] = dnormdxk[0][1][kk] * dproddni[0] + dnormdxk[1][1][kk] * dproddni[1] +
+          dnormdxk[2][1][kk] * dproddni[2];
+      fk[2] = dnormdxk[0][2][kk] * dproddni[0] + dnormdxk[1][2][kk] * dproddni[1] +
+          dnormdxk[2][2][kk] * dproddni[2];
+
+      f[k][0] += fk[0];
+      f[k][1] += fk[1];
+      f[k][2] += fk[2];
+
+      delki[0] = x[k][0] - x[i][0];
+      delki[1] = x[k][1] - x[i][1];
+      delki[2] = x[k][2] - x[i][2];
+
+      if (VFLAG_EITHER) {
+        ev_tally_xyz(k, i, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0], delki[1],
+                     delki[2]);
+      }
+    }
+    f[i][0] +=
+        dnormdxi[0][0] * dproddni[0] + dnormdxi[1][0] * dproddni[1] + dnormdxi[2][0] * dproddni[2];
+    f[i][1] +=
+        dnormdxi[0][1] * dproddni[0] + dnormdxi[1][1] * dproddni[1] + dnormdxi[2][1] * dproddni[2];
+    f[i][2] +=
+        dnormdxi[0][2] * dproddni[0] + dnormdxi[1][2] * dproddni[1] + dnormdxi[2][2] * dproddni[2];
+  }    // loop over ii
+}
+
+/* ----------------------------------------------------------------------
+   Calculate the normals for one atom
+------------------------------------------------------------------------- */
+void PairILPGrapheneHBNOpt::calc_single_normal(int i, int *ILP_neigh, int nneigh, double *normal,
+                                               double (*dnormdri)[3], double (*dnormdrk)[3][3])
+{
+  int j, ii, jj, inum, jnum;
+  int cont, id, ip, m;
+  double nn, xtp, ytp, ztp, delx, dely, delz, nn2;
+  int *ilist, *jlist;
+  double pv12[3], pv31[3], pv23[3], n1[3], dni[3], dnn[3][3], vet[3][3], dpvdri[3][3];
+  double dn1[3][3][3], dpv12[3][3][3], dpv23[3][3][3], dpv31[3][3][3];
+
+  double **x = atom->x;
+
+  for (id = 0; id < 3; id++) {
+    pv12[id] = 0.0;
+    pv31[id] = 0.0;
+    pv23[id] = 0.0;
+    n1[id] = 0.0;
+    dni[id] = 0.0;
+    normal[id] = 0.0;
+    for (ip = 0; ip < 3; ip++) {
+      vet[ip][id] = 0.0;
+      dnn[ip][id] = 0.0;
+      dpvdri[ip][id] = 0.0;
+      dnormdri[ip][id] = 0.0;
+      for (m = 0; m < 3; m++) {
+        dpv12[ip][id][m] = 0.0;
+        dpv31[ip][id][m] = 0.0;
+        dpv23[ip][id][m] = 0.0;
+        dn1[ip][id][m] = 0.0;
+        dnormdrk[ip][id][m] = 0.0;
       }
     }
   }
-  if (vflag_fdotr) virial_fdotr_compute();
-}
 
-inline void PairILPGrapheneHBNOpt::cross_deriv(double *pv, double (*dpv)[3][3], double (*vet)[3],
-                                               int j, int k, int l)
-{
-  pv[0] = vet[j][1] * vet[k][2] - vet[k][1] * vet[j][2];
-  pv[1] = vet[j][2] * vet[k][0] - vet[k][2] * vet[j][0];
-  pv[2] = vet[j][0] * vet[k][1] - vet[k][0] * vet[j][1];
+  xtp = x[i][0];
+  ytp = x[i][1];
+  ztp = x[i][2];
 
-  dpv[j][0][0] = 0.0;
-  dpv[j][0][1] = vet[k][2];
-  dpv[j][0][2] = -vet[k][1];
-  dpv[j][1][0] = -vet[k][2];
-  dpv[j][1][1] = 0.0;
-  dpv[j][1][2] = vet[k][0];
-  dpv[j][2][0] = vet[k][1];
-  dpv[j][2][1] = -vet[k][0];
-  dpv[j][2][2] = 0.0;
+  cont = 0;
+  jlist = ILP_neigh;
+  jnum = nneigh;
+  for (jj = 0; jj < jnum; jj++) {
+    j = jlist[jj];
+    j &= NEIGHMASK;
 
-  dpv[k][0][0] = 0.0;
-  dpv[k][0][1] = -vet[j][2];
-  dpv[k][0][2] = vet[j][1];
-  dpv[k][1][0] = vet[j][2];
-  dpv[k][1][1] = 0.0;
-  dpv[k][1][2] = -vet[j][0];
-  dpv[k][2][0] = -vet[j][1];
-  dpv[k][2][1] = vet[j][0];
-  dpv[k][2][2] = 0.0;
-
-  dpv[l][0][0] = 0.0;
-  dpv[l][0][1] = 0.0;
-  dpv[l][0][2] = 0.0;
-  dpv[l][1][0] = 0.0;
-  dpv[l][1][1] = 0.0;
-  dpv[l][1][2] = 0.0;
-  dpv[l][2][0] = 0.0;
-  dpv[l][2][1] = 0.0;
-  dpv[l][2][2] = 0.0;
-}
-
-inline void PairILPGrapheneHBNOpt::calc_dnormal(double (*dnormal)[3][3], double (*dn1)[3][3],
-                                                double *n1, double nn, double nn2)
-{
-  double dnn[3][3];
-  int m, id, ip;
-  // derivatives of nn, dnn:3x3 vector
-  // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
-  // r[id][m]: the id's component of atom m
-  for (m = 0; m < 3; m++) {
-    for (id = 0; id < 3; id++) {
-      dnn[id][m] = (n1[0] * dn1[m][0][id] + n1[1] * dn1[m][1][id] + n1[2] * dn1[m][2][id]) / nn;
-    }
+    delx = x[j][0] - xtp;
+    dely = x[j][1] - ytp;
+    delz = x[j][2] - ztp;
+    vet[cont][0] = delx;
+    vet[cont][1] = dely;
+    vet[cont][2] = delz;
+    cont++;
   }
-  // dnormal[m][id][ip]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
-  // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
-  for (m = 0; m < 3; m++) {
+
+  if (cont <= 1) {
+    normal[0] = 0.0;
+    normal[1] = 0.0;
+    normal[2] = 1.0;
     for (id = 0; id < 3; id++) {
       for (ip = 0; ip < 3; ip++) {
-        dnormal[m][id][ip] = dn1[m][id][ip] / nn - n1[id] * dnn[ip][m] / nn2;
+        dnormdri[id][ip] = 0.0;
+        for (m = 0; m < 3; m++) { dnormdrk[id][ip][m] = 0.0; }
+      }
+    }
+  } else if (cont == 2) {
+    pv12[0] = vet[0][1] * vet[1][2] - vet[1][1] * vet[0][2];
+    pv12[1] = vet[0][2] * vet[1][0] - vet[1][2] * vet[0][0];
+    pv12[2] = vet[0][0] * vet[1][1] - vet[1][0] * vet[0][1];
+    // derivatives of pv12[0] to ri
+    dpvdri[0][0] = 0.0;
+    dpvdri[0][1] = vet[0][2] - vet[1][2];
+    dpvdri[0][2] = vet[1][1] - vet[0][1];
+    // derivatives of pv12[1] to ri
+    dpvdri[1][0] = vet[1][2] - vet[0][2];
+    dpvdri[1][1] = 0.0;
+    dpvdri[1][2] = vet[0][0] - vet[1][0];
+    // derivatives of pv12[2] to ri
+    dpvdri[2][0] = vet[0][1] - vet[1][1];
+    dpvdri[2][1] = vet[1][0] - vet[0][0];
+    dpvdri[2][2] = 0.0;
+
+    dpv12[0][0][0] = 0.0;
+    dpv12[0][1][0] = vet[1][2];
+    dpv12[0][2][0] = -vet[1][1];
+    dpv12[1][0][0] = -vet[1][2];
+    dpv12[1][1][0] = 0.0;
+    dpv12[1][2][0] = vet[1][0];
+    dpv12[2][0][0] = vet[1][1];
+    dpv12[2][1][0] = -vet[1][0];
+    dpv12[2][2][0] = 0.0;
+
+    // derivatives respect to the second neighbor, atom l
+    dpv12[0][0][1] = 0.0;
+    dpv12[0][1][1] = -vet[0][2];
+    dpv12[0][2][1] = vet[0][1];
+    dpv12[1][0][1] = vet[0][2];
+    dpv12[1][1][1] = 0.0;
+    dpv12[1][2][1] = -vet[0][0];
+    dpv12[2][0][1] = -vet[0][1];
+    dpv12[2][1][1] = vet[0][0];
+    dpv12[2][2][1] = 0.0;
+
+    // derivatives respect to the third neighbor, atom n
+    // derivatives of pv12 to rn is zero
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) { dpv12[id][ip][2] = 0.0; }
+    }
+
+    n1[0] = pv12[0];
+    n1[1] = pv12[1];
+    n1[2] = pv12[2];
+    // the magnitude of the normal vector
+    nn2 = n1[0] * n1[0] + n1[1] * n1[1] + n1[2] * n1[2];
+    nn = sqrt(nn2);
+    if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
+    // the unit normal vector
+    normal[0] = n1[0] / nn;
+    normal[1] = n1[1] / nn;
+    normal[2] = n1[2] / nn;
+    // derivatives of nn, dnn:3x1 vector
+    dni[0] = (n1[0] * dpvdri[0][0] + n1[1] * dpvdri[1][0] + n1[2] * dpvdri[2][0]) / nn;
+    dni[1] = (n1[0] * dpvdri[0][1] + n1[1] * dpvdri[1][1] + n1[2] * dpvdri[2][1]) / nn;
+    dni[2] = (n1[0] * dpvdri[0][2] + n1[1] * dpvdri[1][2] + n1[2] * dpvdri[2][2]) / nn;
+    // derivatives of unit vector ni respect to ri, the result is 3x3 matrix
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) {
+        dnormdri[id][ip] = dpvdri[id][ip] / nn - n1[id] * dni[ip] / nn2;
+      }
+    }
+    // derivatives of non-normalized normal vector, dn1:3x3x3 array
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) {
+        for (m = 0; m < 3; m++) { dn1[id][ip][m] = dpv12[id][ip][m]; }
+      }
+    }
+    // derivatives of nn, dnn:3x3 vector
+    // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+    // r[id][m]: the id's component of atom m
+    for (m = 0; m < 3; m++) {
+      for (id = 0; id < 3; id++) {
+        dnn[id][m] = (n1[0] * dn1[0][id][m] + n1[1] * dn1[1][id][m] + n1[2] * dn1[2][id][m]) / nn;
+      }
+    }
+    // dnormdrk[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+    // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+    for (m = 0; m < 3; m++) {
+      for (id = 0; id < 3; id++) {
+        for (ip = 0; ip < 3; ip++) {
+          dnormdrk[id][ip][m] = dn1[id][ip][m] / nn - n1[id] * dnn[ip][m] / nn2;
+        }
       }
     }
   }
+  //##############################################################################################
+
+  else if (cont == 3) {
+    pv12[0] = vet[0][1] * vet[1][2] - vet[1][1] * vet[0][2];
+    pv12[1] = vet[0][2] * vet[1][0] - vet[1][2] * vet[0][0];
+    pv12[2] = vet[0][0] * vet[1][1] - vet[1][0] * vet[0][1];
+    // derivatives respect to the first neighbor, atom k
+    dpv12[0][0][0] = 0.0;
+    dpv12[0][1][0] = vet[1][2];
+    dpv12[0][2][0] = -vet[1][1];
+    dpv12[1][0][0] = -vet[1][2];
+    dpv12[1][1][0] = 0.0;
+    dpv12[1][2][0] = vet[1][0];
+    dpv12[2][0][0] = vet[1][1];
+    dpv12[2][1][0] = -vet[1][0];
+    dpv12[2][2][0] = 0.0;
+    // derivatives respect to the second neighbor, atom l
+    dpv12[0][0][1] = 0.0;
+    dpv12[0][1][1] = -vet[0][2];
+    dpv12[0][2][1] = vet[0][1];
+    dpv12[1][0][1] = vet[0][2];
+    dpv12[1][1][1] = 0.0;
+    dpv12[1][2][1] = -vet[0][0];
+    dpv12[2][0][1] = -vet[0][1];
+    dpv12[2][1][1] = vet[0][0];
+    dpv12[2][2][1] = 0.0;
+
+    // derivatives respect to the third neighbor, atom n
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) { dpv12[id][ip][2] = 0.0; }
+    }
+
+    pv31[0] = vet[2][1] * vet[0][2] - vet[0][1] * vet[2][2];
+    pv31[1] = vet[2][2] * vet[0][0] - vet[0][2] * vet[2][0];
+    pv31[2] = vet[2][0] * vet[0][1] - vet[0][0] * vet[2][1];
+    // derivatives respect to the first neighbor, atom k
+    dpv31[0][0][0] = 0.0;
+    dpv31[0][1][0] = -vet[2][2];
+    dpv31[0][2][0] = vet[2][1];
+    dpv31[1][0][0] = vet[2][2];
+    dpv31[1][1][0] = 0.0;
+    dpv31[1][2][0] = -vet[2][0];
+    dpv31[2][0][0] = -vet[2][1];
+    dpv31[2][1][0] = vet[2][0];
+    dpv31[2][2][0] = 0.0;
+    // derivatives respect to the third neighbor, atom n
+    dpv31[0][0][2] = 0.0;
+    dpv31[0][1][2] = vet[0][2];
+    dpv31[0][2][2] = -vet[0][1];
+    dpv31[1][0][2] = -vet[0][2];
+    dpv31[1][1][2] = 0.0;
+    dpv31[1][2][2] = vet[0][0];
+    dpv31[2][0][2] = vet[0][1];
+    dpv31[2][1][2] = -vet[0][0];
+    dpv31[2][2][2] = 0.0;
+    // derivatives respect to the second neighbor, atom l
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) { dpv31[id][ip][1] = 0.0; }
+    }
+
+    pv23[0] = vet[1][1] * vet[2][2] - vet[2][1] * vet[1][2];
+    pv23[1] = vet[1][2] * vet[2][0] - vet[2][2] * vet[1][0];
+    pv23[2] = vet[1][0] * vet[2][1] - vet[2][0] * vet[1][1];
+    // derivatives respect to the second neighbor, atom k
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) { dpv23[id][ip][0] = 0.0; }
+    }
+    // derivatives respect to the second neighbor, atom l
+    dpv23[0][0][1] = 0.0;
+    dpv23[0][1][1] = vet[2][2];
+    dpv23[0][2][1] = -vet[2][1];
+    dpv23[1][0][1] = -vet[2][2];
+    dpv23[1][1][1] = 0.0;
+    dpv23[1][2][1] = vet[2][0];
+    dpv23[2][0][1] = vet[2][1];
+    dpv23[2][1][1] = -vet[2][0];
+    dpv23[2][2][1] = 0.0;
+    // derivatives respect to the third neighbor, atom n
+    dpv23[0][0][2] = 0.0;
+    dpv23[0][1][2] = -vet[1][2];
+    dpv23[0][2][2] = vet[1][1];
+    dpv23[1][0][2] = vet[1][2];
+    dpv23[1][1][2] = 0.0;
+    dpv23[1][2][2] = -vet[1][0];
+    dpv23[2][0][2] = -vet[1][1];
+    dpv23[2][1][2] = vet[1][0];
+    dpv23[2][2][2] = 0.0;
+
+    //############################################################################################
+    // average the normal vectors by using the 3 neighboring planes
+    n1[0] = (pv12[0] + pv31[0] + pv23[0]) / cont;
+    n1[1] = (pv12[1] + pv31[1] + pv23[1]) / cont;
+    n1[2] = (pv12[2] + pv31[2] + pv23[2]) / cont;
+    // the magnitude of the normal vector
+    nn2 = n1[0] * n1[0] + n1[1] * n1[1] + n1[2] * n1[2];
+    nn = sqrt(nn2);
+    if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
+    // the unit normal vector
+    normal[0] = n1[0] / nn;
+    normal[1] = n1[1] / nn;
+    normal[2] = n1[2] / nn;
+
+    // for the central atoms, dnormdri is always zero
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) { dnormdri[id][ip] = 0.0; }
+    }
+
+    // derivatives of non-normalized normal vector, dn1:3x3x3 array
+    for (id = 0; id < 3; id++) {
+      for (ip = 0; ip < 3; ip++) {
+        for (m = 0; m < 3; m++) {
+          dn1[id][ip][m] = (dpv12[id][ip][m] + dpv23[id][ip][m] + dpv31[id][ip][m]) / cont;
+        }
+      }
+    }
+    // derivatives of nn, dnn:3x3 vector
+    // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+    // r[id][m]: the id's component of atom m
+    for (m = 0; m < 3; m++) {
+      for (id = 0; id < 3; id++) {
+        dnn[id][m] = (n1[0] * dn1[0][id][m] + n1[1] * dn1[1][id][m] + n1[2] * dn1[2][id][m]) / nn;
+      }
+    }
+    // dnormdrk[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+    // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+    for (m = 0; m < 3; m++) {
+      for (id = 0; id < 3; id++) {
+        for (ip = 0; ip < 3; ip++) {
+          dnormdrk[id][ip][m] = dn1[id][ip][m] / nn - n1[id] * dnn[ip][m] / nn2;
+        }
+      }
+    }
+  } else {
+    error->one(FLERR, "There are too many neighbors for calculating normals");
+  }
+
+  //##############################################################################################
 }
 
-inline void PairILPGrapheneHBNOpt::ev_tally_buffer(int i, int j, double *ei, double *vi, double *ej,
-                                                   double *vj, int nlocal, int newton_pair,
-                                                   double evdwl, double ecoul, double fx, double fy,
-                                                   double fz, double delx, double dely, double delz)
+static bool check_vdw(tagint itag, tagint jtag, double *xi, double *xj)
 {
-  double evdwlhalf, ecoulhalf, epairhalf, v[6];
-  if (eflag_either) {
-    if (eflag_global) {
-      if (newton_pair) {
-        eng_vdwl += evdwl;
-        eng_coul += ecoul;
-      } else {
-        evdwlhalf = 0.5 * evdwl;
-        ecoulhalf = 0.5 * ecoul;
-        if (i < nlocal) {
-          eng_vdwl += evdwlhalf;
-          eng_coul += ecoulhalf;
-        }
-        if (j < nlocal) {
-          eng_vdwl += evdwlhalf;
-          eng_coul += ecoulhalf;
-        }
-      }
-    }
-    if (eflag_atom) {
-      epairhalf = 0.5 * (evdwl + ecoul);
-      if (newton_pair || i < nlocal) *ei += epairhalf;
-      if (newton_pair || j < nlocal) *ej += epairhalf;
-    }
+  if (itag > jtag) {
+    if ((itag + jtag) % 2 == 0) return false;
+  } else if (itag < jtag) {
+    if ((itag + jtag) % 2 == 1) return false;
+  } else {
+    if (xj[2] < xi[2]) return false;
+    if (xj[2] == xi[2] && xj[1] < xi[1]) return false;
+    if (xj[2] == xi[2] && xj[1] == xi[1] && xj[0] < xi[0]) return false;
+  }
+  return true;
+}
+void PairILPGrapheneHBNOpt::update_internal_list()
+{
+  int jnum_sum = 0;
+  int inum = list->inum;
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  int *jlist, jnum;
+  tagint *tag = atom->tag;
+  double **x = atom->x;
+  for (int ii = 0; ii < inum; ii++) { jnum_sum += numneigh[ilist[ii]]; }
+  if (inum > inum_max) {
+    memory->destroy(num_intra);
+    memory->destroy(num_inter);
+    memory->destroy(num_vdw);
+    memory->sfree(first_layered_neigh);
+    //golden ratio grow
+    inum_max = (int) ceil(inum / 0.618);
+    memory->create(num_intra, inum_max, "PairILPGrapheneHBN:intra_layer_count");
+    memory->create(num_inter, inum_max, "PairILPGrapheneHBN:inter_layer_count");
+    memory->create(num_vdw, inum_max, "PairILPGrapheneHBN:vdw_count");
+    first_layered_neigh = (int **) memory->smalloc(inum_max * sizeof(int *),
+                                                   "PairILPGrapheneHBN:first_layered_neigh");
+  }
+  if (jnum_sum > jnum_max) {
+    memory->destroy(layered_neigh);
+    jnum_max = (int) ceil(jnum_sum / 0.618);
+    memory->create(layered_neigh, jnum_max, "PairILPGrapheneHBN:layered_neigh");
   }
 
-  if (vflag_either) {
-    v[0] = delx * fx;
-    v[1] = dely * fy;
-    v[2] = delz * fz;
-    v[3] = delx * fy;
-    v[4] = delx * fz;
-    v[5] = dely * fz;
+  double cut_intra = 0;
+  for (int i = 0; i < nparams; i++) {
+    if (params[i].rcut > cut_intra) { cut_intra = params[i].rcut; }
+  }
+  double cut_intra_listsq = (cut_intra + neighbor->skin) * (cut_intra + neighbor->skin);
 
-    if (vflag_global) {
-      if (newton_pair) {
-        virial[0] += v[0];
-        virial[1] += v[1];
-        virial[2] += v[2];
-        virial[3] += v[3];
-        virial[4] += v[4];
-        virial[5] += v[5];
-      } else {
-        if (i < nlocal) {
-          virial[0] += 0.5 * v[0];
-          virial[1] += 0.5 * v[1];
-          virial[2] += 0.5 * v[2];
-          virial[3] += 0.5 * v[3];
-          virial[4] += 0.5 * v[4];
-          virial[5] += 0.5 * v[5];
-        }
-        if (j < nlocal) {
-          virial[0] += 0.5 * v[0];
-          virial[1] += 0.5 * v[1];
-          virial[2] += 0.5 * v[2];
-          virial[3] += 0.5 * v[3];
-          virial[4] += 0.5 * v[4];
-          virial[5] += 0.5 * v[5];
-        }
+  int total_neigh = 0;
+  for (int ii = 0; ii < inum; ii++) {
+    int i = ilist[ii];
+    tagint itag = tag[i];
+    int jnum = numneigh[i];
+    int *jlist = firstneigh[i];
+    int *jlist_layered = first_layered_neigh[i] = layered_neigh + total_neigh;
+    int ninter = 0, nintra = 0;
+
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = jlist[jj] & NEIGHMASK;
+      if (atom->molecule[j] == atom->molecule[i]) {
+        double delx = x[i][0] - x[j][0];
+        double dely = x[i][1] - x[j][1];
+        double delz = x[i][2] - x[j][2];
+        double rsq = delx * delx + dely * dely + delz * delz;
+        if (rsq < cut_intra_listsq) jlist_layered[nintra++] = j;
       }
     }
-
-    if (vflag_atom) {
-      if (newton_pair || i < nlocal) {
-        vi[0] += 0.5 * v[0];
-        vi[1] += 0.5 * v[1];
-        vi[2] += 0.5 * v[2];
-        vi[3] += 0.5 * v[3];
-        vi[4] += 0.5 * v[4];
-        vi[5] += 0.5 * v[5];
-      }
-      if (newton_pair || j < nlocal) {
-        vj[0] += 0.5 * v[0];
-        vj[1] += 0.5 * v[1];
-        vj[2] += 0.5 * v[2];
-        vj[3] += 0.5 * v[3];
-        vj[4] += 0.5 * v[4];
-        vj[5] += 0.5 * v[5];
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = jlist[jj] & NEIGHMASK;
+      tagint jtag = tag[j];
+      if (atom->molecule[j] != atom->molecule[i]) {
+        if (check_vdw(itag, jtag, x[i], x[j])) jlist_layered[nintra + ninter++] = j;
       }
     }
+    num_vdw[i] = ninter;
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = jlist[jj] & NEIGHMASK;
+      tagint jtag = tag[j];
+      if (atom->molecule[j] != atom->molecule[i]) {
+        if (!check_vdw(itag, jtag, x[i], x[j])) jlist_layered[nintra + ninter++] = j;
+      }
+    }
+    num_intra[i] = nintra;
+    num_inter[i] = ninter;
+    total_neigh += nintra + ninter;
   }
 }
diff --git a/src/OPT/pair_ilp_graphene_hbn_opt.h b/src/OPT/pair_ilp_graphene_hbn_opt.h
index b39d335bf5..9d15b0e600 100644
--- a/src/OPT/pair_ilp_graphene_hbn_opt.h
+++ b/src/OPT/pair_ilp_graphene_hbn_opt.h
@@ -24,7 +24,7 @@ PairStyle(ilp/graphene/hbn/opt,PairILPGrapheneHBNOpt);
 
 namespace LAMMPS_NS {
 
-class PairILPGrapheneHBNOpt : public PairILPGrapheneHBN {
+class PairILPGrapheneHBNOpt : virtual public PairILPGrapheneHBN {
  public:
   PairILPGrapheneHBNOpt(class LAMMPS *);
   ~PairILPGrapheneHBNOpt() override;
@@ -33,13 +33,14 @@ class PairILPGrapheneHBNOpt : public PairILPGrapheneHBN {
   void init_style() override;
 
  protected:
-  void computeILP(int, int);
-  inline void cross_deriv(double *pv, double (*dpv)[3][3], double (*vet)[3], int j, int k, int l);
-  inline void calc_dnormal(double (*dnormal)[3][3], double (*dn1)[3][3], double *n1, double nn,
-                           double nn2);
-  inline void ev_tally_buffer(int i, int j, double *ei, double *vi, double *ej, double *vj,
-                              int nlocal, int newton_pair, double evdwl, double ecoul, double fx,
-                              double fy, double fz, double delx, double dely, double delz);
+  void calc_single_normal(int i, int *ILP_neigh, int nneigh, double *normal, double (*dnormdri)[3],
+                          double (*dnormdrk)[3][3]);
+  void update_internal_list();
+  template <int, int, int> void eval();
+  int *layered_neigh;
+  int **first_layered_neigh;
+  int *num_intra, *num_inter, *num_vdw;
+  int inum_max, jnum_max;
 };
 
 }    // namespace LAMMPS_NS