From 0236cabce99bc31b68b4fd2787c0feb2c810cb8c Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 21 Sep 2020 10:29:42 -0400
Subject: [PATCH] bond/react:add-modify-create-keyword

---
 src/USER-REACTION/fix_bond_react.cpp | 22 +++++++++++++++++++---
 src/USER-REACTION/fix_bond_react.h   |  1 +
 2 files changed, 20 insertions(+), 3 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 995c3a53da..93cf55000a 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -205,6 +205,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   memory->create(limit_duration,nreacts,"bond/react:limit_duration");
   memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
   memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid");
+  memory->create(modify_create_fragid,nreacts,"bond/react:modify_create_fragid");
   memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");
   memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag");
   memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id");
@@ -225,6 +226,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     max_rxn[i] = INT_MAX;
     stabilize_steps_flag[i] = 0;
     custom_charges_fragid[i] = -1;
+    modify_create_fragid[i] = -1;
     // set default limit duration to 60 timesteps
     limit_duration[i] = 60;
     reaction_count[i] = 0;
@@ -380,6 +382,16 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
                                                          "'custom_charges' keyword does not exist");
         }
         iarg += 2;
+      } else if (strcmp(arg[iarg],"modify_create") == 0) {
+        if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
+                                      "'modify_create' has too few arguments");
+        if (strcmp(arg[iarg+1],"no") == 0) modify_create_fragid[rxn] = -1; //default
+        else {
+          modify_create_fragid[rxn] = atom->molecules[unreacted_mol[rxn]]->findfragment(arg[iarg+1]);
+          if (modify_create_fragid[rxn] < 0) error->one(FLERR,"Bond/react: Molecule fragment for "
+                                                         "'modify_create' keyword does not exist");
+        }
+        iarg += 2;
       } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
     }
   }
@@ -561,6 +573,7 @@ FixBondReact::~FixBondReact()
   memory->destroy(var_id);
   memory->destroy(stabilize_steps_flag);
   memory->destroy(custom_charges_fragid);
+  memory->destroy(modify_create_fragid);
 
   memory->destroy(iatomtype);
   memory->destroy(jatomtype);
@@ -2605,8 +2618,7 @@ void FixBondReact::update_everything()
 
   // used when deleting atoms
   int ndel,ndelone;
-  int *mark = new int[nlocal];
-  for (int i = 0; i < nlocal; i++) mark[i] = 0;
+  int *mark;
   tagint *tag = atom->tag;
   AtomVec *avec = atom->avec;
 
@@ -2653,6 +2665,8 @@ void FixBondReact::update_everything()
     delete [] iskip;
 
     // mark to-delete atoms
+    mark = new int[nlocal];
+    for (int i = 0; i < nlocal; i++) mark[i] = 0;
     for (int i = 0; i < update_num_mega; i++) {
       rxnID = update_mega_glove[0][i];
       onemol = atom->molecules[unreacted_mol[rxnID]];
@@ -3150,6 +3164,8 @@ void FixBondReact::update_everything()
   atom->nimpropers += Tdelta_imprp;
 }
 
+
+
 /* ----------------------------------------------------------------------
 read superimpose file
 ------------------------------------------------------------------------- */
@@ -3258,7 +3274,7 @@ void FixBondReact::Equivalences(char *line, int myrxn)
   for (int i = 0; i < nequivalent; i++) {
     readline(line);
     sscanf(line,"%d %d",&tmp1,&tmp2);
-    if (tmp1 > onemol->natoms || tmp2 > onemol->natoms)
+    if (tmp1 > onemol->natoms || tmp2 > twomol->natoms)
       error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
     //equivalences is-> clmn 1: post-reacted, clmn 2: pre-reacted
     equivalences[tmp2-1][0][myrxn] = tmp2;
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 2696d84768..506a98c4ec 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -65,6 +65,7 @@ class FixBondReact : public Fix {
   int custom_exclude_flag;
   int *stabilize_steps_flag;
   int *custom_charges_fragid;
+  int *modify_create_fragid;
   int nconstraints;
   int narrhenius;
   double **constraints;