added buckled flag as custom atom property to angle_mesocnt, used by pair_mesocnt to determine if substitute chain heuristic can be used

src/MESONT/angle_mesocnt.cpp

@@ -20,6 +20,7 @@
 #include "force.h"
 #include "math_const.h"
 #include "memory.h"
+#include "modify.h"
 #include "neighbor.h"
 
 #include <cmath>
@@ -72,6 +73,10 @@ void AngleMesoCNT::compute(int eflag, int vflag)
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
+  int flag, cols;
+  int index = atom->find_custom("buckled", flag, cols);
+  int *buckled = atom->ivector[index];
+
   for (n = 0; n < nanglelist; n++) {
     i1 = anglelist[n][0];
     i2 = anglelist[n][1];
@@ -117,12 +122,16 @@ void AngleMesoCNT::compute(int eflag, int vflag)
       tk = kh[type] * dtheta;
       if (eflag) eangle = tk * dtheta;
       a = -2.0 * tk * s;
+
+      buckled[i2] = 0;
     }
     // bending buckling
     else {
       if (eflag) eangle = kb[type] * fabs(dtheta) + thetab[type] * (kh[type] * thetab[type] - kb[type]);
       if (dtheta < 0) a = kb[type] * s;
       else a = -kb[type] * s;
+
+      buckled[i2] = 1;
     }
     a11 = a * c / rsq1;
     a12 = -a / (r1 * r2);
@@ -208,6 +217,13 @@ void AngleMesoCNT::coeff(int narg, char **arg)
   if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
 }
 
+void AngleMesoCNT::init_style()
+{
+  char *id_fix = utils::strdup("angle_mesocnt_buckled");
+  if (modify->find_fix(id_fix) < 0)
+    modify->add_fix(std::string(id_fix)+" all property/atom i_buckled ghost yes");
+}
+
 /* ---------------------------------------------------------------------- */
 
 double AngleMesoCNT::equilibrium_angle(int i)

src/MESONT/angle_mesocnt.h

@@ -28,6 +28,7 @@ class AngleMesoCNT : public Angle {
   ~AngleMesoCNT() override;
   void compute(int, int) override;
   void coeff(int, char **) override;
+  void init_style() override;
   double equilibrium_angle(int) override;
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;

src/MESONT/pair_mesocnt.cpp

@@ -33,6 +33,7 @@
 
 #include <cmath>
 #include <cstring>
+#include <string>
 #include <unordered_map>
 
 using namespace LAMMPS_NS;
@@ -83,6 +84,7 @@ PairMesoCNT::~PairMesoCNT()
     memory->destroy(numchainlist);
     memory->destroy(nchainlist);
     memory->destroy(endlist);
+    memory->destroy(modelist);
     memory->destroy(chainlist);
 
     memory->destroy(w);
@@ -134,6 +136,7 @@ void PairMesoCNT::allocate()
   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");
@@ -563,6 +566,9 @@ void PairMesoCNT::bond_neigh()
   int **nspecial = atom->nspecial;
   tagint **special = atom->special;
 
+  int flag, cols;
+  int buckled_index = atom->find_custom("buckled", flag, cols);
+
   comm->forward_comm(this);
 
   // create version of atom->special with local ids and correct images
@@ -623,6 +629,7 @@ void PairMesoCNT::bond_neigh()
   memory->destroy(numchainlist);
   memory->destroy(nchainlist);
   memory->destroy(endlist);
+  memory->destroy(modelist);
   memory->destroy(chainlist);
 
   // split neighbor list into neighbor chains based on bond topology
@@ -749,6 +756,7 @@ void PairMesoCNT::bond_neigh()
   // MEMORY: prevent zero-size array creation for chainlist
 
   memory->create_ragged(endlist, nbondlist, numchainlist, "pair:endlist");
+  memory->create_ragged(modelist, nbondlist, numchainlist, "pair:modelist");
   if (empty_neigh)
     memory->create(chainlist, 1, 1, 1, "pair:chainlist");
   else
@@ -777,9 +785,21 @@ void PairMesoCNT::bond_neigh()
       int cstart = chainlist[i][j][0];
       int cend = chainlist[i][j][clen-1];
 
-      if (match_end(type[cstart])) endlist[i][j] = 1;
-      else if (match_end(type[cend])) endlist[i][j] = 2;
+      modelist[i][j] = 0;
+
+      if (nspecial[cstart][0] == 1 && nspecial[cend][0] == 1) modelist[i][j] = 1;
+      else if (nspecial[cstart][0] == 1) endlist[i][j] = 1;
+      else if (nspecial[cend][0] == 1) endlist[i][j] = 2;
       else endlist[i][j] = 0;
+      
+      if (buckled_index != -1) {
+        for (int k = 0; k < nchainlist[i][j]; k++) {
+          if (atom->ivector[buckled_index][chainlist[i][j][k]]) {
+            modelist[i][j] = 1;
+            break;
+          }
+        }
+      }
     }
   }
 

src/MESONT/pair_mesocnt.h

@@ -39,7 +39,7 @@ class PairMesoCNT : public Pair {
   int nend_types;
   int uinf_points, gamma_points, phi_points, usemi_points;
   int *end_types, *reduced_nlist, *numchainlist;
-  int **reduced_neighlist, **nchainlist, **endlist;
+  int **reduced_neighlist, **nchainlist, **endlist, **modelist;
   int ***chainlist;
 
   double ang, ang_inv, eunit, funit;