diff --git a/src/MESONT/pair_mesocnt.cpp b/src/MESONT/pair_mesocnt.cpp
index ecd46e1f17..6751083438 100644
--- a/src/MESONT/pair_mesocnt.cpp
+++ b/src/MESONT/pair_mesocnt.cpp
@@ -107,158 +107,6 @@ void PairMesoCNT::compute(int eflag, int vflag)
 {
   ev_init(eflag, vflag);
 
-  mesolj();
-
-  if (vflag_fdotr) virial_fdotr_compute();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMesoCNT::allocate()
-{
-  allocated = 1;
-  int ntypes = atom->ntypes;
-  int np1 = ntypes + 1;
-  int init_size = 1;
-
-  memory->create(cutsq, np1, np1, "pair:cutsq");
-  memory->create(setflag, np1, np1, "pair:setflag");
-  for (int i = 1; i <= ntypes; i++)
-    for (int j = i; j <= ntypes; j++) setflag[i][j] = 0;
-
-  memory->create(end_types, nend_types, "pair:end_types");
-
-  memory->create(uinf_coeff, uinf_points, 4, "pair:uinf_coeff");
-  memory->create(gamma_coeff, gamma_points, 4, "pair:gamma_coeff");
-  memory->create(phi_coeff, phi_points, phi_points, 4, 4, "pair:phi_coeff");
-  memory->create(usemi_coeff, usemi_points, usemi_points, 4, 4, "pair:usemi_coeff");
-
-  memory->create(numchainlist, init_size, "pair:numchainlist");
-  memory->create(nchainlist, init_size, init_size, "pair:nchainlist");
-  memory->create(endlist, init_size, init_size, "pair:endlist");
-  memory->create(modelist, init_size, init_size, "pair:modelist");
-  memory->create(chainlist, init_size, init_size, init_size, "pair:chainlist");
-
-  memory->create(w, init_size, "pair:w");
-  memory->create(wnode, init_size, "pair:wnode");
-  memory->create(dq_w, init_size, 3, "pair:dq_w");
-  memory->create(q1_dq_w, init_size, 3, 3, "pair:q1_dq_w");
-  memory->create(q2_dq_w, init_size, 3, 3, "pair:q2_dq_w");
-
-  memory->create(param, 7, "pair:param");
-
-  memory->create(flocal, 2, 3, "pair:flocal");
-  memory->create(fglobal, 4, 3, "pair:fglobal");
-  memory->create(basis, 3, 3, "pair:basis");
-}
-
-/* ----------------------------------------------------------------------
-   global settings
-------------------------------------------------------------------------- */
-
-void PairMesoCNT::settings(int narg, char ** /* arg */)
-{
-  if (narg != 0) error->all(FLERR, "Illegal pair_style command");
-}
-
-/* ----------------------------------------------------------------------
-   set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairMesoCNT::coeff(int narg, char **arg)
-{
-  if (narg < 4) error->all(FLERR, "Incorrect args for pair coefficients");
-  read_file(arg[2]);
-
-  nend_types = narg - 3;
-
-  if (!allocated) allocate();
-
-  // end atom types
-  for (int i = 3; i < narg; i++) end_types[i - 3] = utils::inumeric(FLERR, arg[i], false, lmp);
-
-  // units, eV to energy unit conversion
-  ang = force->angstrom;
-  ang_inv = 1.0 / ang;
-  if (strcmp(update->unit_style, "lj") == 0)
-    error->all(FLERR, "Pair style mesocnt does not support lj units");
-  else if (strcmp(update->unit_style, "real") == 0)
-    eunit = 23.06054966;
-  else if (strcmp(update->unit_style, "metal") == 0)
-    eunit = 1.0;
-  else if (strcmp(update->unit_style, "si") == 0)
-    eunit = 1.6021765e-19;
-  else if (strcmp(update->unit_style, "cgs") == 0)
-    eunit = 1.6021765e-12;
-  else if (strcmp(update->unit_style, "electron") == 0)
-    eunit = 3.674932248e-2;
-  else if (strcmp(update->unit_style, "micro") == 0)
-    eunit = 1.6021765e-4;
-  else if (strcmp(update->unit_style, "nano") == 0)
-    eunit = 1.6021765e2;
-  else
-    error->all(FLERR, "Pair style mesocnt does not recognize this units style");
-  funit = eunit * ang_inv;
-
-  // potential variables
-  sig = sig_ang * ang;
-  r = r_ang * ang;
-  rsq = r * r;
-  d = 2.0 * r;
-  d_ang = 2.0 * r_ang;
-  rc = 3.0 * sig;
-  cutoff = rc + d;
-  cutoffsq = cutoff * cutoff;
-  cutoff_ang = cutoff * ang_inv;
-  cutoffsq_ang = cutoff_ang * cutoff_ang;
-  comega = 0.275 * (1.0 - 1.0 / (1.0 + 0.59 * r_ang));
-  ctheta = 0.35 + 0.0226 * (r_ang - 6.785);
-
-  // compute spline coefficients
-  spline_coeff(uinf_data, uinf_coeff, delh_uinf, uinf_points);
-  spline_coeff(gamma_data, gamma_coeff, delh_gamma, gamma_points);
-  spline_coeff(phi_data, phi_coeff, delh_phi, delpsi_phi, phi_points);
-  spline_coeff(usemi_data, usemi_coeff, delh_usemi, delxi_usemi, usemi_points);
-
-  memory->destroy(uinf_data);
-  memory->destroy(gamma_data);
-  memory->destroy(phi_data);
-  memory->destroy(usemi_data);
-
-  int ntypes = atom->ntypes;
-  for (int i = 1; i <= ntypes; i++)
-    for (int j = i; j <= ntypes; j++) setflag[i][j] = 1;
-}
-
-/* ----------------------------------------------------------------------
-   init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairMesoCNT::init_style()
-{
-  if (atom->tag_enable == 0) error->all(FLERR, "Pair style mesocnt requires atom IDs");
-  if (force->newton_pair == 0) error->all(FLERR, "Pair style mesocnt requires newton pair on");
-
-  // need a full neighbor list
-
-  neighbor->add_request(this, NeighConst::REQ_FULL);
-}
-
-/* ----------------------------------------------------------------------
-   init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairMesoCNT::init_one(int /* i */, int /* j */)
-{
-  return cutoff;
-}
-
-/* ----------------------------------------------------------------------
-   calculate energy and forces due to mesoscopic LJ potential
-------------------------------------------------------------------------- */
-
-void PairMesoCNT::mesolj()
-{
   int i, j, k, i1, i2, j1, j2;
   int endflag, endindex;
   int clen, numchain;
@@ -743,6 +591,149 @@ void PairMesoCNT::mesolj()
       }
     }
   }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMesoCNT::allocate()
+{
+  allocated = 1;
+  int ntypes = atom->ntypes;
+  int np1 = ntypes + 1;
+  int init_size = 1;
+
+  memory->create(cutsq, np1, np1, "pair:cutsq");
+  memory->create(setflag, np1, np1, "pair:setflag");
+  for (int i = 1; i <= ntypes; i++)
+    for (int j = i; j <= ntypes; j++) setflag[i][j] = 0;
+
+  memory->create(end_types, nend_types, "pair:end_types");
+
+  memory->create(uinf_coeff, uinf_points, 4, "pair:uinf_coeff");
+  memory->create(gamma_coeff, gamma_points, 4, "pair:gamma_coeff");
+  memory->create(phi_coeff, phi_points, phi_points, 4, 4, "pair:phi_coeff");
+  memory->create(usemi_coeff, usemi_points, usemi_points, 4, 4, "pair:usemi_coeff");
+
+  memory->create(numchainlist, init_size, "pair:numchainlist");
+  memory->create(nchainlist, init_size, init_size, "pair:nchainlist");
+  memory->create(endlist, init_size, init_size, "pair:endlist");
+  memory->create(modelist, init_size, init_size, "pair:modelist");
+  memory->create(chainlist, init_size, init_size, init_size, "pair:chainlist");
+
+  memory->create(w, init_size, "pair:w");
+  memory->create(wnode, init_size, "pair:wnode");
+  memory->create(dq_w, init_size, 3, "pair:dq_w");
+  memory->create(q1_dq_w, init_size, 3, 3, "pair:q1_dq_w");
+  memory->create(q2_dq_w, init_size, 3, 3, "pair:q2_dq_w");
+
+  memory->create(param, 7, "pair:param");
+
+  memory->create(flocal, 2, 3, "pair:flocal");
+  memory->create(fglobal, 4, 3, "pair:fglobal");
+  memory->create(basis, 3, 3, "pair:basis");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::settings(int narg, char ** /* arg */)
+{
+  if (narg != 0) error->all(FLERR, "Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::coeff(int narg, char **arg)
+{
+  if (narg < 4) error->all(FLERR, "Incorrect args for pair coefficients");
+  read_file(arg[2]);
+
+  nend_types = narg - 3;
+
+  if (!allocated) allocate();
+
+  // end atom types
+  for (int i = 3; i < narg; i++) end_types[i - 3] = utils::inumeric(FLERR, arg[i], false, lmp);
+
+  // units, eV to energy unit conversion
+  ang = force->angstrom;
+  ang_inv = 1.0 / ang;
+  if (strcmp(update->unit_style, "lj") == 0)
+    error->all(FLERR, "Pair style mesocnt does not support lj units");
+  else if (strcmp(update->unit_style, "real") == 0)
+    eunit = 23.06054966;
+  else if (strcmp(update->unit_style, "metal") == 0)
+    eunit = 1.0;
+  else if (strcmp(update->unit_style, "si") == 0)
+    eunit = 1.6021765e-19;
+  else if (strcmp(update->unit_style, "cgs") == 0)
+    eunit = 1.6021765e-12;
+  else if (strcmp(update->unit_style, "electron") == 0)
+    eunit = 3.674932248e-2;
+  else if (strcmp(update->unit_style, "micro") == 0)
+    eunit = 1.6021765e-4;
+  else if (strcmp(update->unit_style, "nano") == 0)
+    eunit = 1.6021765e2;
+  else
+    error->all(FLERR, "Pair style mesocnt does not recognize this units style");
+  funit = eunit * ang_inv;
+
+  // potential variables
+  sig = sig_ang * ang;
+  r = r_ang * ang;
+  rsq = r * r;
+  d = 2.0 * r;
+  d_ang = 2.0 * r_ang;
+  rc = 3.0 * sig;
+  cutoff = rc + d;
+  cutoffsq = cutoff * cutoff;
+  cutoff_ang = cutoff * ang_inv;
+  cutoffsq_ang = cutoff_ang * cutoff_ang;
+  comega = 0.275 * (1.0 - 1.0 / (1.0 + 0.59 * r_ang));
+  ctheta = 0.35 + 0.0226 * (r_ang - 6.785);
+
+  // compute spline coefficients
+  spline_coeff(uinf_data, uinf_coeff, delh_uinf, uinf_points);
+  spline_coeff(gamma_data, gamma_coeff, delh_gamma, gamma_points);
+  spline_coeff(phi_data, phi_coeff, delh_phi, delpsi_phi, phi_points);
+  spline_coeff(usemi_data, usemi_coeff, delh_usemi, delxi_usemi, usemi_points);
+
+  memory->destroy(uinf_data);
+  memory->destroy(gamma_data);
+  memory->destroy(phi_data);
+  memory->destroy(usemi_data);
+
+  int ntypes = atom->ntypes;
+  for (int i = 1; i <= ntypes; i++)
+    for (int j = i; j <= ntypes; j++) setflag[i][j] = 1;
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::init_style()
+{
+  if (atom->tag_enable == 0) error->all(FLERR, "Pair style mesocnt requires atom IDs");
+  if (force->newton_pair == 0) error->all(FLERR, "Pair style mesocnt requires newton pair on");
+
+  // need a full neighbor list
+
+  neighbor->add_request(this, NeighConst::REQ_FULL);
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairMesoCNT::init_one(int /* i */, int /* j */)
+{
+  return cutoff;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MESONT/pair_mesocnt.h b/src/MESONT/pair_mesocnt.h
index c58e7439ff..c1100755c7 100644
--- a/src/MESONT/pair_mesocnt.h
+++ b/src/MESONT/pair_mesocnt.h
@@ -76,8 +76,6 @@ class PairMesoCNT : public Pair {
   void finf(const double *, double &, double **);
   void fsemi(const double *, double &, double &, double **);
 
-  void mesolj();
-
   // inlined functions for efficiency
 
   inline void weight(const double *, const double *, const double *, const double *, double &,
diff --git a/src/MESONT/pair_mesocnt_viscous.cpp b/src/MESONT/pair_mesocnt_viscous.cpp
index 9d8e72b443..a288bb8933 100644
--- a/src/MESONT/pair_mesocnt_viscous.cpp
+++ b/src/MESONT/pair_mesocnt_viscous.cpp
@@ -40,7 +40,7 @@ using namespace MathExtra;
 void PairMesoCNTViscous::compute(int eflag, int vflag)
 {
   ev_init(eflag, vflag);
-  
+
   int i, j, k, i1, i2, j1, j2;
   int endflag, endindex;
   int clen, numchain;
@@ -48,7 +48,6 @@ void PairMesoCNTViscous::compute(int eflag, int vflag)
   int **chain;
   double fend, lp, scale, sumw, sumw_inv;
   double evdwl, evdwl_chain;
-  double vtot, fvisc_tot;
   double *r1, *r2, *q1, *q2, *qe;
   double *vr1, *vr2, *vq1, *vq2;
   double vr[3], vp1[3], vp2[3], vp[3], vrel[3], fvisc[3];
@@ -83,7 +82,7 @@ void PairMesoCNTViscous::compute(int eflag, int vflag)
 
     r1 = x[i1];
     r2 = x[i2];
-
+    
     vr1 = v[i1];
     vr2 = v[i2];
     add3(vr1, vr2, vr);
@@ -100,274 +99,464 @@ void PairMesoCNTViscous::compute(int eflag, int vflag)
       clen = nchain[j];
       if (clen < 2) continue;
 
-      // assign end position
+      if (modelist[i][j]) {
+        for (k = 0; k < clen - 1; k++) {
+          j1 = chain[j][k];
+          j2 = chain[j][k + 1];
+          j1 &= NEIGHMASK;
+          j2 &= NEIGHMASK;
+          q1 = x[j1];
+          q2 = x[j2];
 
-      endflag = end[j];
-      if (endflag == 1) {
-        endindex = chain[j][0];
-        qe = x[endindex];
-      } else if (endflag == 2) {
-        endindex = chain[j][clen - 1];
-        qe = x[endindex];
+          geometry(r1, r2, q1, q2, q1, p, m, param, basis);
+          if (param[0] > cutoff) continue;
+          
+          double calpha = cos(param[1]);
+          double salpha = sin(param[1]);
+          double hsq = param[0] * param[0];
+
+          double ceta = calpha * param[4];
+          double seta = salpha * param[4];
+          double dsq1 = hsq + seta * seta;
+          if (ceta < param[2]) dsq1 += pow(ceta - param[2], 2);
+          else if (ceta > param[3]) dsq1 += pow(ceta - param[3], 2);
+
+          ceta = calpha * param[5];
+          seta = salpha * param[5];
+          
+          double dsq2 = hsq + seta * seta;
+          if (ceta < param[2]) dsq2 += pow(ceta - param[2], 2);
+          else if (ceta > param[3]) dsq2 += pow(ceta - param[3], 2);
+
+          if (dsq1 > cutoffsq && dsq2 > cutoffsq) continue;
+
+          int jj1, jj2;
+          
+          if (dsq1 < dsq2) {
+            jj1 = j1;
+            jj2 = j2;
+          }
+          else {
+            if (param[1] > MY_PI) param[1] -= MY_PI;
+            else param[1] += MY_PI;
+
+            double temp = -param[5];
+            param[5] = -param[4];
+            param[4] = temp;
+            param[6] = temp;
+
+            negate3(m);
+          
+            jj1 = j2;
+            jj2 = j1;
+          }
+
+          // first force contribution
+
+          fsemi(param, evdwl, fend, flocal);          
+
+          if (evdwl == 0.0) continue;
+          
+          // transform to global coordinate system
+
+          matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
+          matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
+
+          // forces acting on segment 2
+
+          add3(fglobal[0], fglobal[1], ftotal);
+          scaleadd3(fend, m, ftotal, ftotal);
+          scale3(-0.5, ftotal);
+
+          sub3(r1, p, delr1);
+          sub3(r2, p, delr2);
+          cross3(delr1, fglobal[0], t1);
+          cross3(delr2, fglobal[1], t2);
+          add3(t1, t2, torque);
+        
+          cross3(torque, m, ftorque);
+          lp = param[5] - param[4];
+          scale3(1.0 / lp, ftorque);
+          
+          add3(ftotal, ftorque, fglobal[2]);
+          sub3(ftotal, ftorque, fglobal[3]);
+
+          // add forces to nodes
+
+          scaleadd3(0.5, fglobal[0], f[i1], f[i1]);
+          scaleadd3(0.5, fglobal[1], f[i2], f[i2]);
+          scaleadd3(0.5, fglobal[2], f[jj1], f[jj1]);
+          scaleadd3(0.5, fglobal[3], f[jj2], f[jj2]);
+          scaleadd3(0.5 * fend, m, f[jj1], f[jj1]);
+          
+          // add energy
+
+          if (eflag_either) {
+            evdwl = 0.5 * evdwl;
+            double evdwl_atom = 0.25 * evdwl;
+            if (eflag_global) eng_vdwl += evdwl;
+            if (eflag_atom) {
+              eatom[i1] += evdwl_atom;
+              eatom[i2] += evdwl_atom;
+              eatom[jj1] += evdwl_atom;
+              eatom[jj2] += evdwl_atom;
+            }
+          }
+          
+          // second force contribution
+
+          param[6] += lp;
+          
+          fsemi(param, evdwl, fend, flocal);
+          
+          if (evdwl == 0.0) continue;
+
+          // transform to global coordinate system
+
+          matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
+          matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
+
+          // forces acting on segment 2
+
+          add3(fglobal[0], fglobal[1], ftotal);
+          scaleadd3(fend, m, ftotal, ftotal);
+          scale3(-0.5, ftotal);
+
+          scaleadd3(lp, m, p, p);
+
+          sub3(r1, p, delr1);
+          sub3(r2, p, delr2);
+          cross3(delr1, fglobal[0], t1);
+          cross3(delr2, fglobal[1], t2);
+          add3(t1, t2, torque);
+        
+          cross3(torque, m, ftorque);
+          scale3(1.0 / lp, ftorque);
+          
+          add3(ftotal, ftorque, fglobal[2]);
+          sub3(ftotal, ftorque, fglobal[3]);
+
+          // add forces to nodes
+
+          scaleadd3(-0.5, fglobal[0], f[i1], f[i1]);
+          scaleadd3(-0.5, fglobal[1], f[i2], f[i2]);
+          scaleadd3(-0.5, fglobal[2], f[jj2], f[jj2]);
+          scaleadd3(0.5, fglobal[3], f[jj1], f[jj1]);
+          sub3(f[jj2], fglobal[3], f[jj2]);
+          scaleadd3(-0.5 * fend, m, f[jj2], f[jj2]);
+
+          // add energy
+
+          if (eflag_either) {
+            evdwl = -0.5 * evdwl;
+            double evdwl_atom = 0.25 * evdwl;
+            if (eflag_global) eng_vdwl += evdwl;
+            if (eflag_atom) {
+              eatom[i1] += evdwl_atom;
+              eatom[i2] += evdwl_atom;
+              eatom[jj1] += evdwl_atom;
+              eatom[jj2] += evdwl_atom;
+            }
+          }
+
+        }
       }
+      else {
+        // assign end position
 
-      // compute substitute straight (semi-)infinite CNT
-
-      zero3(p1);
-      zero3(p2);
-      zero3(dr1_sumw);
-      zero3(dr2_sumw);
-      zero3(vp1);
-      zero3(vp2);
-      zeromat3(q1_dr1_w);
-      zeromat3(q2_dr1_w);
-      zeromat3(q1_dr2_w);
-      zeromat3(q2_dr2_w);
-      for (k = 0; k < clen; k++) {
-        wnode[k] = 0.0;
-        zero3(dq_w[k]);
-        zeromat3(q1_dq_w[k]);
-        zeromat3(q2_dq_w[k]);
-      }
-      sumw = 0.0;
-
-      for (k = 0; k < clen - 1; k++) {
-        j1 = chain[j][k];
-        j2 = chain[j][k + 1];
-        j1 &= NEIGHMASK;
-        j2 &= NEIGHMASK;
-        q1 = x[j1];
-        q2 = x[j2];
-        vq1 = v[j1];
-        vq2 = v[j2];
-
-        weight(r1, r2, q1, q2, w[k], dr1_w, dr2_w, dq1_w, dq2_w);
-
-        if (w[k] == 0.0) {
-          if (endflag == 1 && k == 0)
-            endflag = 0;
-          else if (endflag == 2 && k == clen - 2)
-            endflag = 0;
-          continue;
+        endflag = end[j];
+        if (endflag == 1) {
+          endindex = chain[j][0];
+          qe = x[endindex];
+        } else if (endflag == 2) {
+          endindex = chain[j][clen - 1];
+          qe = x[endindex];
         }
 
-        sumw += w[k];
-        wnode[k] += w[k];
-        wnode[k + 1] += w[k];
+        // compute substitute straight (semi-)infinite CNT
 
-        scaleadd3(w[k], q1, p1, p1);
-        scaleadd3(w[k], q2, p2, p2);
+        zero3(p1);
+        zero3(p2);
+        zero3(dr1_sumw);
+        zero3(dr2_sumw);
+        zeromat3(q1_dr1_w);
+        zeromat3(q2_dr1_w);
+        zeromat3(q1_dr2_w);
+        zeromat3(q2_dr2_w);
+        for (k = 0; k < clen; k++) {
+          wnode[k] = 0.0;
+          zero3(dq_w[k]);
+          zeromat3(q1_dq_w[k]);
+          zeromat3(q2_dq_w[k]);
+        }
+        sumw = 0.0;
+
+        for (k = 0; k < clen - 1; k++) {
+          j1 = chain[j][k];
+          j2 = chain[j][k + 1];
+          j1 &= NEIGHMASK;
+          j2 &= NEIGHMASK;
+          q1 = x[j1];
+          q2 = x[j2];
+
+          weight(r1, r2, q1, q2, w[k], dr1_w, dr2_w, dq1_w, dq2_w);
+
+          if (w[k] == 0.0) {
+            if (endflag == 1 && k == 0)
+              endflag = 0;
+            else if (endflag == 2 && k == clen - 2)
+              endflag = 0;
+            continue;
+          }
+
+          sumw += w[k];
+          wnode[k] += w[k];
+          wnode[k + 1] += w[k];
+
+          scaleadd3(w[k], q1, p1, p1);
+          scaleadd3(w[k], q2, p2, p2);
         
-        // weighted velocity for friction
+          // weighted velocity for friction
 
-        scaleadd3(w[k], vq1, vp1, vp1);
-        scaleadd3(w[k], vq2, vp2, vp2);
+          scaleadd3(w[k], vq1, vp1, vp1);
+          scaleadd3(w[k], vq2, vp2, vp2);
 
-        // weight gradient terms
+          // weight gradient terms
 
-        add3(dr1_w, dr1_sumw, dr1_sumw);
-        add3(dr2_w, dr2_sumw, dr2_sumw);
+          add3(dr1_w, dr1_sumw, dr1_sumw);
+          add3(dr2_w, dr2_sumw, dr2_sumw);
 
-        outer3(q1, dr1_w, temp);
-        plus3(temp, q1_dr1_w, q1_dr1_w);
-        outer3(q2, dr1_w, temp);
-        plus3(temp, q2_dr1_w, q2_dr1_w);
-        outer3(q1, dr2_w, temp);
-        plus3(temp, q1_dr2_w, q1_dr2_w);
-        outer3(q2, dr2_w, temp);
-        plus3(temp, q2_dr2_w, q2_dr2_w);
+          outer3(q1, dr1_w, temp);
+          plus3(temp, q1_dr1_w, q1_dr1_w);
+          outer3(q2, dr1_w, temp);
+          plus3(temp, q2_dr1_w, q2_dr1_w);
+          outer3(q1, dr2_w, temp);
+          plus3(temp, q1_dr2_w, q1_dr2_w);
+          outer3(q2, dr2_w, temp);
+          plus3(temp, q2_dr2_w, q2_dr2_w);
 
-        add3(dq1_w, dq_w[k], dq_w[k]);
-        add3(dq2_w, dq_w[k + 1], dq_w[k + 1]);
+          add3(dq1_w, dq_w[k], dq_w[k]);
+          add3(dq2_w, dq_w[k + 1], dq_w[k + 1]);
 
-        outer3(q1, dq1_w, temp);
-        plus3(temp, q1_dq_w[k], q1_dq_w[k]);
-        outer3(q1, dq2_w, temp);
-        plus3(temp, q1_dq_w[k + 1], q1_dq_w[k + 1]);
-        outer3(q2, dq1_w, temp);
-        plus3(temp, q2_dq_w[k], q2_dq_w[k]);
-        outer3(q2, dq2_w, temp);
-        plus3(temp, q2_dq_w[k + 1], q2_dq_w[k + 1]); 
-      }
+          outer3(q1, dq1_w, temp);
+          plus3(temp, q1_dq_w[k], q1_dq_w[k]);
+          outer3(q1, dq2_w, temp);
+          plus3(temp, q1_dq_w[k + 1], q1_dq_w[k + 1]);
+          outer3(q2, dq1_w, temp);
+          plus3(temp, q2_dq_w[k], q2_dq_w[k]);
+          outer3(q2, dq2_w, temp);
+          plus3(temp, q2_dq_w[k + 1], q2_dq_w[k + 1]);
+        }
 
-      if (sumw == 0.0) continue;
+        if (sumw == 0.0) continue;
 
-      sumw_inv = 1.0 / sumw;
-      scale3(sumw_inv, p1);
-      scale3(sumw_inv, p2);
+        sumw_inv = 1.0 / sumw;
+        scale3(sumw_inv, p1);
+        scale3(sumw_inv, p2);
 
-      // compute geometry and forces
+        // compute geometry and forces
 
-      // infinite CNT case
+        if (endflag == 0) {
+          
+          // infinite CNT case
+          
+          geometry(r1, r2, p1, p2, nullptr, p, m, param, basis);
+          
+          if (param[0] > cutoff) continue;
+          double salpha = sin(param[1]);
+          double sxi1 = salpha * param[2];
+          double sxi2 = salpha * param[3];
+          double hsq = param[0] * param[0];
+          if (sxi1 * sxi1 + hsq > cutoffsq && sxi2 * sxi2 + hsq > cutoffsq)
+            continue;
 
-      if (endflag == 0) {
-        geometry(r1, r2, p1, p2, nullptr, p, m, param, basis);
-        if (param[0] > cutoff) continue;
-        finf(param, evdwl, flocal);
+          finf(param, evdwl, flocal);
 
-        // semi-infinite CNT case with end at start of chain
+        } else {
+          
+          // semi-infinite CNT case
+          // endflag == 1: CNT end at start of chain
+          // endflag == 2: CNT end at end of chain
+        
+          if (endflag == 1) geometry(r1, r2, p1, p2, qe, p, m, param, basis);
+          else geometry(r1, r2, p2, p1, qe, p, m, param, basis);
+          
+          if (param[0] > cutoff) continue;
+          double hsq = param[0] * param[0];
+          double calpha = cos(param[1]);
+          double etamin = calpha * param[2];
+          double dsq1;
+          if (etamin < param[6]) {
+            dsq1 = hsq + pow(sin(param[1]) * param[6], 2);
+            dsq1 += pow(param[2] - calpha * param[6], 2);
+          }
+          else
+            dsq1 = hsq + pow(sin(param[1]) * param[2], 2);
+          
+          etamin = calpha * param[3];
+          double dsq2;
+          if (etamin < param[6]) {
+            dsq2 = hsq + pow(sin(param[1]) * param[6], 2);
+            dsq2 += pow(param[3] - calpha * param[6], 2);
+          }
+          else
+            dsq2 = hsq + pow(sin(param[1]) * param[3], 2);
 
-      } else if (endflag == 1) {
-        geometry(r1, r2, p1, p2, qe, p, m, param, basis);
-        if (param[0] > cutoff) continue;
-        fsemi(param, evdwl, fend, flocal);
+          if (dsq1 > cutoffsq && dsq2 > cutoffsq) continue;
 
-        // semi-infinite CNT case with end at end of chain
+          fsemi(param, evdwl, fend, flocal);
+        }
 
-      } else {
-        geometry(r1, r2, p2, p1, qe, p, m, param, basis);
-        if (param[0] > cutoff) continue;
-        fsemi(param, evdwl, fend, flocal);
-      }
+        if (evdwl == 0.0) continue;
+        evdwl *= 0.5;
 
-      evdwl *= 0.5;
+        // transform to global coordinate system
 
-      // transform to global coordinate system
-
-      matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
-      matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
+        matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
+        matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
       
-      // mean chain velocity and relative velocity
+        // mean chain velocity and relative velocity
 
-      add3(vp1, vp2, vp);
-      scale3(0.5*sumw_inv, vp);
-      sub3(vp, vr, vrel);
+        add3(vp1, vp2, vp);
+        scale3(0.5*sumw_inv, vp);
+        sub3(vp, vr, vrel);
 
-      // friction forces
-
-      vtot = len3(vrel);
-      if (vtot < vswitch) scale3(0.25*a1, vrel, fvisc);
-      else {
-        fvisc_tot = b2 / vtot - a2;
-        if (fvisc_tot < 0.0) zero3(fvisc);
-        else scale3(0.25*fvisc_tot, vrel, fvisc);
-      }
-
-      add3(fvisc, f[i1], f[i1]);
-      add3(fvisc, f[i2], f[i2]);
-
-      // forces acting on approximate chain
-
-      add3(fglobal[0], fglobal[1], ftotal);
-      if (endflag) scaleadd3(fend, m, ftotal, ftotal);
-      scale3(-0.5, ftotal);
-
-      sub3(r1, p, delr1);
-      sub3(r2, p, delr2);
-      cross3(delr1, fglobal[0], t1);
-      cross3(delr2, fglobal[1], t2);
-      add3(t1, t2, torque);
-
-      // additional torque contribution from chain end
-
-      if (endflag) {
-        sub3(qe, p, delqe);
-        cross3(delqe, m, t3);
-        scale3(fend, t3);
-        add3(t3, torque, torque);
-      }
-
-      cross3(torque, m, ftorque);
-      lp = param[5] - param[4];
-      scale3(1.0 / lp, ftorque);
-
-      if (endflag == 2) {
-        add3(ftotal, ftorque, fglobal[3]);
-        sub3(ftotal, ftorque, fglobal[2]);
-      } else {
-        add3(ftotal, ftorque, fglobal[2]);
-        sub3(ftotal, ftorque, fglobal[3]);
-      }
-
-      scale3(0.5, fglobal[0]);
-      scale3(0.5, fglobal[1]);
-      scale3(0.5, fglobal[2]);
-      scale3(0.5, fglobal[3]);
-
-      // weight gradient terms acting on current segment
-
-      outer3(p1, dr1_sumw, temp);
-      minus3(q1_dr1_w, temp, dr1_p1);
-      outer3(p2, dr1_sumw, temp);
-      minus3(q2_dr1_w, temp, dr1_p2);
-      outer3(p1, dr2_sumw, temp);
-      minus3(q1_dr2_w, temp, dr2_p1);
-      outer3(p2, dr2_sumw, temp);
-      minus3(q2_dr2_w, temp, dr2_p2);
-
-      transpose_matvec(dr1_p1, fglobal[2], fgrad_r1_p1);
-      transpose_matvec(dr1_p2, fglobal[3], fgrad_r1_p2);
-      transpose_matvec(dr2_p1, fglobal[2], fgrad_r2_p1);
-      transpose_matvec(dr2_p2, fglobal[3], fgrad_r2_p2);
-
-      // add forces to nodes in current segment
-
-      add3(fglobal[0], f[i1], f[i1]);
-      add3(fglobal[1], f[i2], f[i2]);
-
-      scaleadd3(sumw_inv, fgrad_r1_p1, f[i1], f[i1]);
-      scaleadd3(sumw_inv, fgrad_r1_p2, f[i1], f[i1]);
-      scaleadd3(sumw_inv, fgrad_r2_p1, f[i2], f[i2]);
-      scaleadd3(sumw_inv, fgrad_r2_p2, f[i2], f[i2]);
-
-      // add forces in approximate chain
-
-      for (k = 0; k < clen - 1; k++) {
-        if (w[k] == 0.0) continue;
-        j1 = chain[j][k];
-        j2 = chain[j][k + 1];
-        j1 &= NEIGHMASK;
-        j2 &= NEIGHMASK;
-        scale = w[k] * sumw_inv;
-        scaleadd3(scale, fglobal[2], f[j1], f[j1]);
-        scaleadd3(scale, fglobal[3], f[j2], f[j2]);
-      
         // friction forces
+
+        vtot = len3(vrel);
+        if (vtot < vswitch) scale3(0.25*a1, vrel, fvisc);
+        else {
+          fvisc_tot = b2 / vtot - a2;
+          if (fvisc_tot < 0.0) zero3(fvisc);
+          else scale3(0.25*fvisc_tot, vrel, fvisc);
+        }
+
+        add3(fvisc, f[i1], f[i1]);
+        add3(fvisc, f[i2], f[i2]);
+
+        // forces acting on approximate chain
+
+        add3(fglobal[0], fglobal[1], ftotal);
+        if (endflag) scaleadd3(fend, m, ftotal, ftotal);
+        scale3(-0.5, ftotal);
+
+        sub3(r1, p, delr1);
+        sub3(r2, p, delr2);
+        cross3(delr1, fglobal[0], t1);
+        cross3(delr2, fglobal[1], t2);
+        add3(t1, t2, torque);
+
+        // additional torque contribution from chain end
+
+        if (endflag) {
+          sub3(qe, p, delqe);
+          cross3(delqe, m, t3);
+          scale3(fend, t3);
+          add3(t3, torque, torque);
+        }
+
+        cross3(torque, m, ftorque);
+        lp = param[5] - param[4];
+        scale3(1.0 / lp, ftorque);
+
+        if (endflag == 2) {
+          add3(ftotal, ftorque, fglobal[3]);
+          sub3(ftotal, ftorque, fglobal[2]);
+        } else {
+          add3(ftotal, ftorque, fglobal[2]);
+          sub3(ftotal, ftorque, fglobal[3]);
+        }
+
+        scale3(0.5, fglobal[0]);
+        scale3(0.5, fglobal[1]);
+        scale3(0.5, fglobal[2]);
+        scale3(0.5, fglobal[3]);
+
+        // weight gradient terms acting on current segment
+
+        outer3(p1, dr1_sumw, temp);
+        minus3(q1_dr1_w, temp, dr1_p1);
+        outer3(p2, dr1_sumw, temp);
+        minus3(q2_dr1_w, temp, dr1_p2);
+        outer3(p1, dr2_sumw, temp);
+        minus3(q1_dr2_w, temp, dr2_p1);
+        outer3(p2, dr2_sumw, temp);
+        minus3(q2_dr2_w, temp, dr2_p2);
+
+        transpose_matvec(dr1_p1, fglobal[2], fgrad_r1_p1);
+        transpose_matvec(dr1_p2, fglobal[3], fgrad_r1_p2);
+        transpose_matvec(dr2_p1, fglobal[2], fgrad_r2_p1);
+        transpose_matvec(dr2_p2, fglobal[3], fgrad_r2_p2);
+
+        // add forces to nodes in current segment
+
+        add3(fglobal[0], f[i1], f[i1]);
+        add3(fglobal[1], f[i2], f[i2]);
+
+        scaleadd3(sumw_inv, fgrad_r1_p1, f[i1], f[i1]);
+        scaleadd3(sumw_inv, fgrad_r1_p2, f[i1], f[i1]);
+        scaleadd3(sumw_inv, fgrad_r2_p1, f[i2], f[i2]);
+        scaleadd3(sumw_inv, fgrad_r2_p2, f[i2], f[i2]);
+
+        // add forces in approximate chain
+
+        for (k = 0; k < clen - 1; k++) {
+          if (w[k] == 0.0) continue;
+          j1 = chain[j][k];
+          j2 = chain[j][k + 1];
+          j1 &= NEIGHMASK;
+          j2 &= NEIGHMASK;
+          scale = w[k] * sumw_inv;
+          scaleadd3(scale, fglobal[2], f[j1], f[j1]);
+          scaleadd3(scale, fglobal[3], f[j2], f[j2]);
         
-        scaleadd3(-scale, fvisc, f[j1], f[j1]);
-        scaleadd3(-scale, fvisc, f[j2], f[j2]);
-      }
+          // friction forces
+          
+          scaleadd3(-scale, fvisc, f[j1], f[j1]);
+          scaleadd3(-scale, fvisc, f[j2], f[j2]);
+        }
 
-      // weight gradient terms acting on approximate chain
-      // iterate over nodes instead of segments
+        // weight gradient terms acting on approximate chain
+        // iterate over nodes instead of segments
 
-      for (k = 0; k < clen; k++) {
-        if (wnode[k] == 0.0) continue;
-        j1 = chain[j][k];
-        j1 &= NEIGHMASK;
+        for (k = 0; k < clen; k++) {
+          if (wnode[k] == 0.0) continue;
+          j1 = chain[j][k];
+          j1 &= NEIGHMASK;
 
-        outer3(p1, dq_w[k], temp);
-        minus3(q1_dq_w[k], temp, dq_p1);
-        outer3(p2, dq_w[k], temp);
-        minus3(q2_dq_w[k], temp, dq_p2);
+          outer3(p1, dq_w[k], temp);
+          minus3(q1_dq_w[k], temp, dq_p1);
+          outer3(p2, dq_w[k], temp);
+          minus3(q2_dq_w[k], temp, dq_p2);
 
-        transpose_matvec(dq_p1, fglobal[2], fgrad_q_p1);
-        transpose_matvec(dq_p2, fglobal[3], fgrad_q_p2);
+          transpose_matvec(dq_p1, fglobal[2], fgrad_q_p1);
+          transpose_matvec(dq_p2, fglobal[3], fgrad_q_p2);
 
-        scaleadd3(sumw_inv, fgrad_q_p1, f[j1], f[j1]);
-        scaleadd3(sumw_inv, fgrad_q_p2, f[j1], f[j1]);
-      }
+          scaleadd3(sumw_inv, fgrad_q_p1, f[j1], f[j1]);
+          scaleadd3(sumw_inv, fgrad_q_p2, f[j1], f[j1]);
+        }
 
-      // force on node at CNT end
+        // force on node at CNT end
 
-      if (endflag) scaleadd3(0.5 * fend, m, f[endindex], f[endindex]);
+        if (endflag) scaleadd3(0.5 * fend, m, f[endindex], f[endindex]);
 
-      // compute energy
+        // compute energy
 
-      if (eflag_either) {
-        if (eflag_global) eng_vdwl += evdwl;
-        if (eflag_atom) {
-          eatom[i1] += 0.25 * evdwl;
-          eatom[i2] += 0.25 * evdwl;
-          for (k = 0; k < clen - 1; k++) {
-            if (w[k] == 0.0) continue;
-            j1 = chain[j][k];
-            j2 = chain[j][k + 1];
-            j1 &= NEIGHMASK;
-            j2 &= NEIGHMASK;
-            evdwl_chain = 0.5 * w[k] * sumw_inv * evdwl;
-            eatom[j1] += evdwl_chain;
-            eatom[j2] += evdwl_chain;
+        if (eflag_either) {
+          if (eflag_global) eng_vdwl += evdwl;
+          if (eflag_atom) {
+            eatom[i1] += 0.25 * evdwl;
+            eatom[i2] += 0.25 * evdwl;
+            for (k = 0; k < clen - 1; k++) {
+              if (w[k] == 0.0) continue;
+              j1 = chain[j][k];
+              j2 = chain[j][k + 1];
+              j1 &= NEIGHMASK;
+              j2 &= NEIGHMASK;
+              evdwl_chain = 0.5 * w[k] * sumw_inv * evdwl;
+              eatom[j1] += evdwl_chain;
+              eatom[j2] += evdwl_chain;
+            }
           }
         }
       }