diff --git a/src/MESONT/pair_mesocnt.cpp b/src/MESONT/pair_mesocnt.cpp
index f8df017fc4..a5cb261bbc 100644
--- a/src/MESONT/pair_mesocnt.cpp
+++ b/src/MESONT/pair_mesocnt.cpp
@@ -148,7 +148,7 @@ void PairMesoCNT::compute(int eflag, int vflag)
 
   // update bond neighbor list when necessary
 
-  if (update->ntimestep == neighbor->lastcall) bond_neigh();
+  if (update->ntimestep == neighbor->lastcall) bond_neigh_topo();
 
   // iterate over all bonds
 
@@ -804,10 +804,155 @@ double PairMesoCNT::init_one(int /* i */, int /* j */)
 }
 
 /* ----------------------------------------------------------------------
-   update bond neighbor lists
+   update bond neighbor lists based on atom and mol IDs
 ------------------------------------------------------------------------- */
 
-void PairMesoCNT::bond_neigh()
+void PairMesoCNT::bond_neigh_id()
+{
+  int nlocal = atom->nlocal;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+
+  int *numneigh = list->numneigh;
+  int *numneigh_max;
+  memory->create(numneigh_max, nbondlist, "pair:numneigh_max");
+
+  for (int i = 0; i < nbondlist; i++) {
+    int i1 = bondlist[i][0];
+    int i2 = bondlist[i][1];
+    int numneigh1, numneigh2;
+
+    // prevent ghost atom with undefined neighbors
+
+    if (i1 > nlocal - 1)
+      numneigh1 = 0;
+    else
+      numneigh1 = numneigh[i1];
+    if (i2 > nlocal - 1)
+      numneigh2 = 0;
+    else
+      numneigh2 = numneigh[i2];
+    numneigh_max[i] = numneigh1 + numneigh2 + 2;
+  }
+
+  // create temporary arrays for chain creation
+
+  memory->create(reduced_nlist, nbondlist, "pair:reduced_nlist");
+  memory->create_ragged(reduced_neighlist, nbondlist, numneigh_max, "pair:reduced_neighlist");
+
+  // reduce neighbors to common list and find longest common list size
+
+  for (int i = 0; i < nbondlist; i++) {
+    int i1 = bondlist[i][0];
+    int i2 = bondlist[i][1];
+
+    neigh_common(i1, i2, reduced_nlist[i], reduced_neighlist[i]);
+
+    // sort list according to atom-id
+
+    sort(reduced_neighlist[i], reduced_nlist[i]);
+  }
+
+  // resize chain arrays
+
+  memory->destroy(numchainlist);
+  memory->destroy(nchainlist);
+  memory->destroy(endlist);
+  memory->destroy(chainlist);
+  
+  memory->grow(selfid, nbondlist, "pair:selfid");
+  memory->grow(selfpos, nbondlist, 2, "pair:selfpos");
+
+  // count neighbor chains per bond
+
+  memory->create(numchainlist, nbondlist, "pair:numchainlist");
+
+  int numchain_max = 0;
+  for (int i = 0; i < nbondlist; i++) {
+    numchainlist[i] = count_chains(reduced_neighlist[i], reduced_nlist[i]);
+    if (numchain_max < numchainlist[i]) numchain_max = numchainlist[i];
+  }
+
+  // count neighbor chain lengths per bond
+
+  memory->create_ragged(nchainlist, nbondlist, numchainlist, "pair:nchainlist");
+
+  for (int i = 0; i < nbondlist; i++)
+    chain_lengths(reduced_neighlist[i], reduced_nlist[i], nchainlist[i]);
+
+  // create connected neighbor chains and check for ends
+  // MEMORY: prevent zero-size array creation for chainlist
+
+  memory->create_ragged(endlist, nbondlist, numchainlist, "pair:endlist");
+  if (numchain_max > 0)
+    memory->create_ragged(chainlist, nbondlist, numchainlist, nchainlist, "pair:chainlist");
+  else
+    memory->create(chainlist, 1, 1, 1, "pair:chainlist");
+
+  int nchain_max = 0;
+  for (int i = 0; i < nbondlist; i++) {
+    int *reduced_neigh = reduced_neighlist[i];
+    int *end = endlist[i];
+    int *nchain = nchainlist[i];
+    int **chain = chainlist[i];
+
+    // set up connected chains and check for ends
+
+    chain_split(reduced_neigh, reduced_nlist[i], nchain, chain, end);
+
+    // find longest chain
+
+    for (int j = 0; j < numchainlist[i]; j++)
+      if (nchain_max < nchain[j]) nchain_max = nchain[j];
+
+    // find selfid and selfpos
+    
+    tagint tag1 = atom->tag[bondlist[i][0]];
+    tagint tag2 = atom->tag[bondlist[i][1]];
+    
+    bool found1 = false;
+    bool found2 = false;
+
+    for (int j = 0; j < numchainlist[i]; j++) {
+      for (int k = 0; k < nchain[j]; k++) {
+        tagint temp_tag = atom->tag[chain[j][k]];
+        if (tag1 == temp_tag) {
+          selfid[i] = j;
+          selfpos[i][0] = k;
+          found1 = true;
+        }
+        else if (tag2 == temp_tag) {
+          selfid[i] = j;
+          selfpos[i][1] = k;
+          found2 = true;
+        }
+
+        if (found1 && found2) break;
+      }
+      if (found1 && found2) break;
+    }
+  }
+
+  // resize potential arrays
+
+  memory->grow(w, nchain_max, "pair:w");
+  memory->grow(wnode, nchain_max, "pair:wnode");
+  memory->grow(dq_w, nchain_max, 3, "pair:dq_w");
+  memory->grow(q1_dq_w, nchain_max, 3, 3, "pair:q1_dq_w");
+  memory->grow(q2_dq_w, nchain_max, 3, 3, "pair:q2_dq_w");
+
+  // destroy temporary arrays for chain creation
+
+  memory->destroy(numneigh_max);
+  memory->destroy(reduced_neighlist);
+  memory->destroy(reduced_nlist);
+}
+
+/* ----------------------------------------------------------------------
+   update bond neighbor lists based on bond topology
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::bond_neigh_topo()
 {
   int nlocal = atom->nlocal;
   int nghost = atom->nghost;
@@ -1120,6 +1265,123 @@ void PairMesoCNT::neigh_common(int i1, int i2, int &numred, int *redlist)
   }
 }
 
+/* ----------------------------------------------------------------------
+   count neighbor chains of given bond
+------------------------------------------------------------------------- */
+
+int PairMesoCNT::count_chains(int *redlist, int numred)
+{
+  if (numred == 0) return 0;
+
+  tagint *tag = atom->tag;
+  tagint *mol = atom->molecule;
+  int count = 1;
+
+  // split neighbor list and count chains
+
+  for (int j = 0; j < numred - 1; j++) {
+    int j1 = redlist[j];
+    int j2 = redlist[j + 1];
+    if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) count++;
+  }
+
+  return count;
+}
+
+/* ----------------------------------------------------------------------
+   count lengths of neighbor chains of given bond
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::chain_lengths(int *redlist, int numred, int *nchain)
+{
+  if (numred == 0) return;
+
+  tagint *tag = atom->tag;
+  tagint *mol = atom->molecule;
+  int clen = 0;
+  int cid = 0;
+
+  // split neighbor list into connected chains
+
+  for (int j = 0; j < numred - 1; j++) {
+    int j1 = redlist[j];
+    int j2 = redlist[j + 1];
+    clen++;
+    if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) {
+      nchain[cid++] = clen;
+      clen = 0;
+    }
+  }
+  clen++;
+  nchain[cid++] = clen;
+}
+
+/* ----------------------------------------------------------------------
+   split neighbors into chains and identify ends
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::chain_split(int *redlist, int numred, int *nchain, int **chain, int *end)
+{
+  if (numred == 0) return;
+
+  tagint *tag = atom->tag;
+  tagint *mol = atom->molecule;
+  tagint *type = atom->type;
+  int clen = 0;
+  int cid = 0;
+
+  // split neighbor list into connected chains
+
+  for (int j = 0; j < numred - 1; j++) {
+    int j1 = redlist[j];
+    int j2 = redlist[j + 1];
+    chain[cid][clen++] = j1;
+    if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) {
+      cid++;
+      clen = 0;
+    }
+  }
+  chain[cid][clen++] = redlist[numred - 1];
+  cid++;
+
+  // check for chain ends
+
+  for (int j = 0; j < cid; j++) {
+    int clen = nchain[j];
+    int cstart = chain[j][0];
+    int cend = chain[j][clen - 1];
+
+    bool estart = match_end(type[cstart]);
+    bool eend = match_end(type[cend]);
+
+    if (estart && eend) end[j] = 3;
+    else if (estart) end[j] = 1;
+    else if (eend) end[j] = 2;
+    else end[j] = 0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   insertion sort list according to corresponding atom ID
+------------------------------------------------------------------------- */
+
+void PairMesoCNT::sort(int *list, int size)
+{
+  int i, j, temp1, temp2;
+  tagint *tag = atom->tag;
+  for (i = 1; i < size; i++) {
+    j = i;
+    temp1 = list[j - 1];
+    temp2 = list[j];
+    while (j > 0 && tag[temp1] > tag[temp2]) {
+      list[j] = temp1;
+      list[j - 1] = temp2;
+      j--;
+      temp1 = list[j - 1];
+      temp2 = list[j];
+    }
+  }
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/MESONT/pair_mesocnt.h b/src/MESONT/pair_mesocnt.h
index 8d8f5e6d83..661c02b380 100644
--- a/src/MESONT/pair_mesocnt.h
+++ b/src/MESONT/pair_mesocnt.h
@@ -63,9 +63,16 @@ class PairMesoCNT : public Pair {
 
   bool match_end(int);
 
+  int count_chains(int *, int);
+
   void allocate();
-  void bond_neigh();
+  void bond_neigh_id();
+  void bond_neigh_topo();
   void neigh_common(int, int, int &, int *);
+  void chain_lengths(int *, int, int *);
+  void chain_split(int *, int, int *, int **, int *);
+  void sort(int *, int);
+  
   void read_file(const char *);
   void read_data(PotentialFileReader &, double *, double &, double &, int);
   void read_data(PotentialFileReader &, double **, double &, double &, double &, double &, int);
diff --git a/src/MESONT/pair_mesocnt_viscous.cpp b/src/MESONT/pair_mesocnt_viscous.cpp
index 8475473cbd..6eadfc4129 100644
--- a/src/MESONT/pair_mesocnt_viscous.cpp
+++ b/src/MESONT/pair_mesocnt_viscous.cpp
@@ -82,7 +82,7 @@ void PairMesoCNTViscous::compute(int eflag, int vflag)
 
   // update bond neighbor list when necessary
 
-  if (update->ntimestep == neighbor->lastcall) bond_neigh();
+  if (update->ntimestep == neighbor->lastcall) bond_neigh_topo();
 
   // iterate over all bonds