Merge pull request #3337 from vpalkar/pair-srp-react

Adding pair style srp/react

src/.gitignore

@@ -1522,6 +1522,8 @@
 /fix_smd_wall_surface.h
 /fix_srp.cpp
 /fix_srp.h
+/fix_srp_react.cpp
+/fix_srp_react.h
 /fix_tfmc.cpp
 /fix_tfmc.h
 /fix_ttm.cpp
@@ -1562,6 +1564,8 @@
 /pair_smd_ulsph.h
 /pair_srp.cpp
 /pair_srp.h
+/pair_srp_react.cpp
+/pair_srp_react.h
 /pair_thole.cpp
 /pair_thole.h
 /pair_buck_mdf.cpp

src/MC/fix_bond_break.h

@@ -25,6 +25,7 @@ FixStyle(bond/break,FixBondBreak);
 namespace LAMMPS_NS {
 
 class FixBondBreak : public Fix {
+  friend class FixSRPREACT;
  public:
   FixBondBreak(class LAMMPS *, int, char **);
   ~FixBondBreak() override;

src/MC/fix_bond_create.h

@@ -25,6 +25,7 @@ FixStyle(bond/create,FixBondCreate);
 namespace LAMMPS_NS {
 
 class FixBondCreate : public Fix {
+  friend class FixSRPREACT;
  public:
   FixBondCreate(class LAMMPS *, int, char **);
   ~FixBondCreate() override;

src/MISC/fix_srp.cpp

@@ -68,6 +68,8 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   btype = -1;
   bptype = -1;
 
+  pair_name = "srp";
+
   // zero
   for (int i = 0; i < atom->nmax; i++)
     for (int m = 0; m < 2; m++)
@@ -101,15 +103,11 @@ int FixSRP::setmask()
 
 void FixSRP::init()
 {
-  if (force->pair_match("hybrid",1) == nullptr && force->pair_match("hybrid/overlay",1) == nullptr)
-    error->all(FLERR,"Cannot use pair srp without pair_style hybrid");
+  if (!force->pair_match("^hybrid",0))
+    error->all(FLERR,"Cannot use pair {} without pair_style hybrid", pair_name);
 
-  int has_rigid = 0;
-  for (int i = 0; i < modify->nfix; i++)
-    if (utils::strmatch(modify->fix[i]->style,"^rigid")) ++has_rigid;
-
-  if (has_rigid > 0)
-    error->all(FLERR,"Pair srp is not compatible with rigid fixes.");
+  if (modify->get_fix_by_style("^rigid").size() > 0)
+    error->all(FLERR,"Pair {} is not compatible with rigid fixes.", pair_name);
 
   if ((bptype < 1) || (bptype > atom->ntypes))
     error->all(FLERR,"Illegal bond particle type");
@@ -118,11 +116,10 @@ void FixSRP::init()
   // because this fix's pre_exchange() creates per-atom data structure
   // that data must be current for atom migration to carry it along
 
-  for (int i = 0; i < modify->nfix; i++) {
-    if (modify->fix[i] == this) break;
-    if (modify->fix[i]->pre_exchange_migrate)
-      error->all(FLERR,"Fix SRP comes after a fix which "
-                 "migrates atoms in pre_exchange");
+  for (auto ifix : modify->get_fix_list()) {
+    if (ifix == this) break;
+    if (ifix->pre_exchange_migrate)
+      error->all(FLERR,"Fix {} comes after a fix which migrates atoms in pre_exchange", style);
   }
 
   // setup neigh exclusions for diff atom types
@@ -274,8 +271,8 @@ void FixSRP::setup_pre_force(int /*zz*/)
 
   // stop if cutghost is insufficient
   if (cutneighmax_srp > cutghostmin)
-    error->all(FLERR,"Communication cutoff too small for fix srp. "
-            "Need {:.8}, current {:.8}", cutneighmax_srp, cutghostmin);
+    error->all(FLERR,"Communication cutoff too small for fix {}. "
+               "Need {:.8}, current {:.8}", style, cutneighmax_srp, cutghostmin);
 
   // assign tags for new atoms, update map
   atom->tag_extend();
@@ -345,11 +342,11 @@ void FixSRP::pre_exchange()
     if (atom->type[ii] != bptype) continue;
 
     i = atom->map(static_cast<tagint>(array[ii][0]));
-    if (i < 0) error->all(FLERR,"Fix SRP failed to map atom");
+    if (i < 0) error->all(FLERR,"Fix {} failed to map atom", style);
     i = domain->closest_image(ii,i);
 
     j = atom->map(static_cast<tagint>(array[ii][1]));
-    if (j < 0) error->all(FLERR,"Fix SRP failed to map atom");
+    if (j < 0) error->all(FLERR,"Fix {} failed to map atom", style);
     j = domain->closest_image(ii,j);
 
     // position of bond particle ii

src/MISC/fix_srp.h

@@ -54,9 +54,10 @@ class FixSRP : public Fix {
 
   double **array;
 
- private:
+ protected:
   int btype;
   int bptype;
+  std::string pair_name;
 };
 
 }    // namespace LAMMPS_NS

src/MISC/fix_srp_react.cpp

@@ -0,0 +1,145 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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: Vaibhav Palkar (Kuksenok Lab, Clemson University)
+   based on the pair srp code by Timothy Sirk (ARL)
+
+This pair style interfaces the pair style srp with the stochastic reaction
+fixes bond/break and bond/create by updating pseudo beads corresponding to
+bonds as bond breaking and formation takes place. This is useful in
+simulation of reactions in polymers with soft potentials such as DPD.
+
+See the doc page for pair_style srp/react command for usage instructions.
+
+There is an example script for this package in examples/PACKAGES/srp_react/.
+------------------------------------------------------------------------- */
+
+#include "fix_srp_react.h"
+
+#include "modify.h"
+#include "neighbor.h"
+
+#include "fix_bond_break.h"
+#include "fix_bond_create.h"
+
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+// clang-format on
+/* ---------------------------------------------------------------------- */
+
+FixSRPREACT::FixSRPREACT(LAMMPS *lmp, int narg, char **arg) :
+    FixSRP(lmp, narg, arg), f_bb(nullptr), idbreak(nullptr), f_bc(nullptr), idcreate(nullptr)
+{
+  pair_name = "srp/react";
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixSRPREACT::~FixSRPREACT()
+{
+  // free memory
+  delete[] idbreak;
+  delete[] idcreate;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixSRPREACT::setmask()
+{
+  int mask = 0;
+  mask |= PRE_FORCE;
+  mask |= PRE_EXCHANGE;
+  mask |= POST_RUN;
+  mask |= POST_NEIGHBOR;
+
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSRPREACT::init()
+{
+  FixSRP::init();
+
+  // find fix bond break
+  if (idbreak) f_bb = (FixBondBreak *) modify->get_fix_by_id(idbreak);
+
+  // find fix bond create
+  if (idcreate) f_bc = (FixBondCreate *) modify->get_fix_by_id(idcreate);
+}
+
+/* ----------------------------------------------------------------------
+   rebuild bond particle array
+------------------------------------------------------------------------- */
+void FixSRPREACT::post_neighbor()
+{
+  // store ncalls as it is reset in fix srp setup pre force
+  int ncalls = neighbor->ncalls;
+
+  if (f_bb) {
+    if (f_bb->breakcount) {
+      setup_pre_force(0);
+
+      //reset break count before exiting
+      // not reseting breakcount would lead to redundant rebuilds
+      f_bb->breakcount = 0;
+
+      // count additional call during setup_pre_force
+      neighbor->ncalls = ncalls + 1;
+    }
+  }
+
+  if (f_bc) {
+    if (f_bc->createcount) {
+      setup_pre_force(0);
+
+      //reset create count before exiting
+      // not reseting createcount would lead to redundant rebuilds
+      f_bc->createcount = 0;
+
+      // count additional call during setup_pre_force
+      neighbor->ncalls = ncalls + 1;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   interface with pair class
+   in addition to bond type and bond particle type,
+     pair srp react sets id of either fix bond break or bond create
+------------------------------------------------------------------------- */
+
+int FixSRPREACT::modify_param(int /*narg*/, char **arg)
+{
+  if (strcmp(arg[0], "btype") == 0) {
+    btype = utils::inumeric(FLERR, arg[1], false, lmp);
+    return 2;
+  }
+  if (strcmp(arg[0], "bptype") == 0) {
+    bptype = utils::inumeric(FLERR, arg[1], false, lmp);
+    return 2;
+  }
+  if (strcmp(arg[0], "bond/break") == 0) {
+    idbreak = utils::strdup(arg[1]);
+    return 2;
+  }
+  if (strcmp(arg[0], "bond/create") == 0) {
+    idcreate = utils::strdup(arg[1]);
+    return 2;
+  }
+  return 0;
+}

src/MISC/fix_srp_react.h

@@ -0,0 +1,44 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(SRPREACT,FixSRPREACT);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_SRP_REACT_H
+#define LMP_FIX_SRP_REACT_H
+
+#include "fix_srp.h"
+
+namespace LAMMPS_NS {
+
+class FixSRPREACT : public FixSRP {
+ public:
+  FixSRPREACT(class LAMMPS *, int, char **);
+  ~FixSRPREACT() override;
+  int setmask() override;
+  void init() override;
+  void post_neighbor() override;
+  int modify_param(int, char **) override;
+
+ private:
+  class FixBondBreak *f_bb;
+  char *idbreak;
+  class FixBondCreate *f_bc;
+  char *idcreate;
+};
+}    // namespace LAMMPS_NS
+#endif
+#endif

src/MISC/pair_srp.cpp

@@ -33,7 +33,6 @@ Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix.h"
 #include "fix_srp.h"
 #include "force.h"
 #include "memory.h"
@@ -84,7 +83,7 @@ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp), fix_id(nullptr)
   //   will be invoked before other fixes that migrate atoms
   //   this is checked for in FixSRP
 
-  f_srp = dynamic_cast<FixSRP *>( modify->add_fix(fmt::format("{:02d}_FIX_SRP all SRP",srp_instance)));
+  f_srp = dynamic_cast<FixSRP *>(modify->add_fix(fmt::format("{:02d}_FIX_SRP all SRP", srp_instance)));
   ++srp_instance;
 }
 
@@ -126,7 +125,7 @@ PairSRP::~PairSRP()
   }
 
   // check nfix in case all fixes have already been deleted
-  if (modify->nfix) modify->delete_fix(f_srp->id);
+  if (modify->nfix && modify->get_fix_by_id(f_srp->id)!=nullptr) modify->delete_fix(f_srp->id);
 }
 
 /* ----------------------------------------------------------------------
@@ -439,7 +438,7 @@ void PairSRP::coeff(int narg, char **arg)
 void PairSRP::init_style()
 {
   if (!force->newton_pair)
-    error->all(FLERR,"PairSRP: Pair srp requires newton pair on");
+    error->all(FLERR,"Pair srp requires newton pair on");
 
   // verify that fix SRP is still defined and has not been changed.
 

src/MISC/pair_srp.h

@@ -40,7 +40,7 @@ class PairSRP : public Pair {
   void write_restart_settings(FILE *) override;
   void read_restart_settings(FILE *) override;
 
- private:
+ protected:
   inline void onetwoexclude(int *&, int &, int *&, int *&, int **&);
   inline void remapBonds(int &);
   void allocate();

src/MISC/pair_srp_react.cpp

@@ -0,0 +1,230 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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: Vaibhav Palkar (Kuksenok Lab, Clemson University)
+   based on the pair srp code by Timothy Sirk (ARL)
+
+This pair style interfaces the pair style srp with the stochastic reaction
+fixes bond/break and bond/create by updating pseudo beads corresponding to
+bonds as bond breaking and formation takes place. This is useful in
+simulation of reactions in polymers with soft potentials such as DPD.
+
+See the doc page for pair_style srp/react command for usage instructions.
+
+There is an example script for this package in examples/PACKAGES/srp_react/.
+------------------------------------------------------------------------- */
+
+#include "pair_srp_react.h"
+
+#include "atom.h"
+#include "citeme.h"
+#include "comm.h"
+#include "error.h"
+#include "fix_srp_react.h"
+#include "force.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+#include "output.h"
+#include "thermo.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+static const char cite_srpreact[] =
+  "@Article{palkar2022\n"
+  " author = {Palkar, Vaibhav and Kuksenok, Olga},\n"
+  " title = {Controlling Degradation and Erosion of Polymer Networks: Insights from Mesoscale Modeling},\n"
+  " journal = {J. Phys. Chem. B},\n"
+  " year = {2022},\n"
+  " volume = {126},\n"
+  " pages = {336-346}\n"
+  "}\n\n";
+
+static int srp_instance = 0;
+
+/* ----------------------------------------------------------------------
+ constructor
+ ---------------------------------------------------------------------- */
+
+PairSRPREACT::PairSRPREACT(LAMMPS *lmp) :
+  PairSRP(lmp), idbreak(nullptr), idcreate(nullptr), bond_break(false), bond_create(false)
+{
+
+  if (lmp->citeme) lmp->citeme->add(cite_srpreact);
+
+  // pair srp/react has its own fix, hence delete fix srp instance
+  // created in the constructor of pair srp
+  for (auto ifix : modify->get_fix_by_style("SRP"))
+    modify->delete_fix(ifix->id);
+
+  // similar to fix SRP, create fix SRP REACT instance here with unique fix id
+  f_srp = (FixSRPREACT *) modify->add_fix(fmt::format("{:02d}_FIX_SRP_REACT all SRPREACT",srp_instance));
+  ++srp_instance;
+}
+
+PairSRPREACT::~PairSRPREACT()
+{
+  delete[] idbreak;
+  delete[] idcreate;
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+ ------------------------------------------------------------------------- */
+
+void PairSRPREACT::settings(int narg, char **arg)
+{
+  if (narg < 3 || narg > 8)
+    error->all(FLERR,"Illegal pair_style command");
+
+  if (atom->tag_enable == 0)
+    error->all(FLERR,"Pair_style srp/react requires atom IDs");
+
+  cut_global = utils::numeric(FLERR,arg[0],false,lmp);
+  // wildcard
+  if (strcmp(arg[1],"*") == 0) {
+    btype = 0;
+  } else {
+    btype = utils::inumeric(FLERR,arg[1],false,lmp);
+    if ((btype > atom->nbondtypes) || (btype <= 0))
+      error->all(FLERR,"Illegal pair_style command");
+  }
+
+  // settings
+  midpoint = false;
+  min = false;
+
+  if (strcmp(arg[2],"min") == 0) min = true;
+  else if (strcmp(arg[2],"mid") == 0) midpoint = true;
+  else
+    error->all(FLERR,"Illegal pair_style command");
+
+  // default for bond/break and bond/create settings
+  bond_create=false;
+  bond_break=false;
+  idbreak = nullptr;
+  idcreate= nullptr;
+
+  // find whether id is of bond/break or bond/create
+  const char *reactid = arg[3];
+  if (utils::strmatch(modify->get_fix_by_id(reactid)->style,"^bond/break")) {
+    bond_break = true;
+    idbreak = utils::strdup(reactid);
+  } else if (utils::strmatch(modify->get_fix_by_id(reactid)->style,"^bond/create")) {
+    bond_create = true;
+    idcreate = utils::strdup(reactid);
+  } else error->all(FLERR,"Illegal pair_style command");
+
+  int iarg = 4;
+  // default exclude 1-2
+  // scaling for 1-2, etc not supported
+  exclude = 1;
+
+  // use last atom type by default for bond particles
+  bptype = atom->ntypes;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"exclude") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command");
+      exclude = utils::logical(FLERR, arg[iarg+1], false, lmp);
+      if (min && !exclude) error->all(FLERR,"Illegal exclude option in pair srp command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"bptype") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command");
+      bptype = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+      if ((bptype < 1) || (bptype > atom->ntypes))
+        error->all(FLERR,"Illegal bond particle type for srp");
+      iarg += 2;
+   } else error->all(FLERR,"Illegal pair srp command");
+  }
+
+  // reset cutoffs if explicitly set
+  if (allocated) {
+    int i,j;
+    for (i = 1; i <= bptype; i++)
+      for (j = i; j <= bptype; j++)
+        if (setflag[i][j]) cut[i][j] = cut_global;
+  }
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+ ------------------------------------------------------------------------- */
+
+void PairSRPREACT::init_style()
+{
+  if (!force->newton_pair)
+    error->all(FLERR,"Pair srp/react requires newton pair on");
+
+  // verify that fix SRP is still defined and has not been changed.
+
+  if (strcmp(f_srp->style,"SRPREACT") != 0)
+    error->all(FLERR,"Fix SRPREACT has been changed unexpectedly");
+
+  if (comm->me == 0)
+    utils::logmesg(lmp,"Using type {} for bond particles\n",bptype);
+
+  // set bond and bond particle types in fix srp
+  // bonds of this type will be represented by bond particles
+  // if bond type is 0, then all bonds have bond particles
+  // btype = bond type
+
+  char c0[20];
+  char* arg0[2];
+  sprintf(c0, "%d", btype);
+  arg0[0] = (char *) "btype";
+  arg0[1] = c0;
+  f_srp->modify_params(2, arg0);
+
+  // bptype = bond particle type
+  sprintf(c0, "%d", bptype);
+  arg0[0] = (char *) "bptype";
+  arg0[1] = c0;
+  f_srp->modify_params(2, arg0);
+
+  // if using fix bond/break, set id of fix bond/break in fix srp
+  // idbreak = id of fix bond break
+  if (bond_break) {
+    sprintf(c0, "%s", idbreak);
+    arg0[0] = (char *) "bond/break";
+    arg0[1] = c0;
+    f_srp->modify_params(2, arg0);
+  }
+
+  // if using fix bond/create, set id of fix bond/create in fix srp
+  // idcreate = id of fix bond break
+  if (bond_create) {
+    sprintf(c0, "%s", idcreate);
+    arg0[0] = (char *) "bond/create";
+    arg0[1] = c0;
+    f_srp->modify_params(2, arg0);
+  }
+
+  // bond particles do not contribute to energy or virial
+  // bond particles do not belong to group all
+  // but thermo normalization is by nall
+  // therefore should turn off normalization
+  char *arg1[2];
+  arg1[0] = (char *) "norm";
+  arg1[1] = (char *) "no";
+  output->thermo->modify_params(2, arg1);
+  if (comm->me == 0) error->message(FLERR,"Thermo normalization turned off by pair srp/react");
+
+  neighbor->request(this,instance_me);
+}

src/MISC/pair_srp_react.h

@@ -0,0 +1,41 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(srp/react,PairSRPREACT);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_SRP_REACT_H
+#define LMP_PAIR_SRP_REACT_H
+
+#include "pair_srp.h"
+
+namespace LAMMPS_NS {
+
+class PairSRPREACT : public PairSRP {
+ public:
+  PairSRPREACT(class LAMMPS *);
+  ~PairSRPREACT() override;
+  void settings(int, char **) override;
+  void init_style() override;
+
+ private:
+  char *idbreak;
+  char *idcreate;
+  bool bond_create, bond_break;
+};
+}    // namespace LAMMPS_NS
+#endif
+#endif