Initial commit for pair style srp-react files

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.h

@@ -54,7 +54,7 @@ class FixSRP : public Fix {
 
   double **array;
 
- private:
+ protected:
   int btype;
   int bptype;
 };

src/MISC/fix_srp_react.cpp

@@ -0,0 +1,186 @@
+// 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 "fix_srp_react.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "modify.h"
+#include "neighbor.h"
+
+#include "fix_bond_break.h"
+#include "fix_bond_create.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixSRPREACT::FixSRPREACT(LAMMPS *lmp, int narg, char **arg) : FixSRP(lmp, narg, arg)
+{
+  // default idbreak and idcreate = NULL
+  idbreak = nullptr;
+  idcreate = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixSRPREACT::setmask()
+{
+  int mask = 0;
+  mask |= PRE_FORCE;
+  mask |= PRE_EXCHANGE;
+  mask |= POST_RUN;
+  mask |= POST_NEIGHBOR;
+
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSRPREACT::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");
+
+  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 ((bptype < 1) || (bptype > atom->ntypes))
+    error->all(FLERR,"Illegal bond particle type");
+
+  // this fix must come before any fix which migrates atoms in its pre_exchange()
+  // 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");
+  }
+
+  // setup neigh exclusions for diff atom types
+  // bond particles do not interact with other types
+  // type bptype only interacts with itself
+
+  for (int z = 1; z < atom->ntypes; z++) {
+    if (z == bptype)
+      continue;
+    neighbor->modify_params(fmt::format("exclude type {} {}",z,bptype));
+  }
+
+  // find fix bond break
+  if( idbreak != nullptr )
+    f_bb = (FixBondBreak *) modify->get_fix_by_id(idbreak);
+  
+  // find fix bond create
+  if( idcreate != nullptr )
+    f_bc = (FixBondCreate *) modify->get_fix_by_id(idcreate);
+
+  // free memory
+  delete [] idbreak;
+  delete [] 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( idbreak != nullptr)
+    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( idcreate != nullptr)
+    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) {
+    int n = strlen(arg[1]) + 1;
+    idbreak = new char[n];
+    strcpy(idbreak,arg[1]);
+    return 2;
+  }
+  if (strcmp(arg[0],"bond/create") == 0) {
+    int n = strlen(arg[1]) + 1;
+    idcreate = new char[n];
+    strcpy(idcreate,arg[1]);
+    return 2;
+  }
+  return 0;
+}

src/MISC/fix_srp_react.h

@@ -0,0 +1,50 @@
+/* -*- 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 **);
+  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
+
+/* ERROR/WARNING messages:
+
+*/

src/MISC/pair_srp.cpp

@@ -126,7 +126,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);
 }
 
 /* ----------------------------------------------------------------------

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,239 @@
+// 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 "domain.h"
+#include "error.h"
+#include "fix.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)
+{
+
+  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( int ifix = 0; ifix<modify->nfix; ifix++)
+    if( strcmp(modify->get_fix_by_index(ifix)->style, "SRP") == 0)
+      modify->delete_fix(ifix);
+
+  // 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;
+}
+
+/* ----------------------------------------------------------------------
+ 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 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(strcmp(modify->get_fix_by_id(reactid)->style,"bond/break") == 0)
+  {
+    bond_break = true;
+    int n = strlen(reactid) + 1;
+    idbreak = new char[n];
+    strcpy(idbreak,reactid);
+  }
+  else if(strcmp(modify->get_fix_by_id(reactid)->style,"bond/create") == 0)
+  {
+    bond_create = true;
+    int n = strlen(reactid) + 1;
+    idcreate = new char[n];
+    strcpy(idcreate,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,"PairSRPREACT: Pair srp 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 SRP 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);
+    delete [] idbreak;
+  }
+  
+  // 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);
+    delete [] idcreate;
+  }
+
+  // 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");
+
+  neighbor->request(this,instance_me);
+}

src/MISC/pair_srp_react.h

@@ -0,0 +1,46 @@
+/* -*- 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 *);
+  void settings(int, char **) override;
+  void init_style() override;
+
+ private:
+  char *idbreak;
+  char *idcreate;
+  bool bond_create, bond_break;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/