From aaa8e9d21953eaede0b81cc8f5bcf82623964209 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Wed, 9 Nov 2022 17:51:21 -0500
Subject: [PATCH] populate atoms2bond in bond/react

---
 src/REACTION/fix_bond_react.cpp | 44 +++++++++++++++++++++++++++++++--
 src/REACTION/fix_bond_react.h   | 10 ++++++--
 2 files changed, 50 insertions(+), 4 deletions(-)

diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp
index 6d19a8a7a9..a2a0c00301 100644
--- a/src/REACTION/fix_bond_react.cpp
+++ b/src/REACTION/fix_bond_react.cpp
@@ -888,6 +888,8 @@ void FixBondReact::post_integrate()
   for (int i = 0; i < nreacts; i++) {
     nattempt[i] = 0;
   }
+  // reset per-bond compute map flag
+  atoms2bondflag = 0;
 
   int nlocal = atom->nlocal;
   int nall = atom->nlocal + atom->nghost;
@@ -2346,6 +2348,11 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string&
     if (cperbond->size_local_cols > 0) error->one(FLERR,"Bond/react: 'Compute bond/local' used by reaction "
                                          "special function 'rxnbond' must compute one value");
 
+    if (atoms2bondflag == 0) {
+      atoms2bondflag = 1;
+      get_atoms2bond(cperbond->groupbit);
+    }
+
     nsum = 0;
     for (int i = 0; i < onemol->natoms; i++) {
       if (onemol->fragmentmask[ifrag][i]) {
@@ -2359,8 +2366,8 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string&
     if (cperbond->invoked_local != lmp->update->ntimestep)
       cperbond->compute_local();
 
-    it = cperbond->atoms2bond.find(aset);
-    if (it == cperbond->atoms2bond.end()) error->one(FLERR,"Bond/react: Unable to locate bond referenced by "
+    it = atoms2bond.find(aset);
+    if (it == atoms2bond.end()) error->one(FLERR,"Bond/react: Unable to locate bond referenced by "
                                             "reaction special function 'rxnbond'");
     ibond = it->second;
     perbondval = cperbond->vector_local[ibond];
@@ -2406,6 +2413,39 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string&
   return 0.0;
 }
 
+/* ----------------------------------------------------------------------
+populate map to get bond index from atom IDs
+------------------------------------------------------------------------- */
+
+void FixBondReact::get_atoms2bond(int cgroupbit)
+{
+  int i,m,atom1,atom2,btype,nb;
+  std::set<tagint> aset;
+
+  int nlocal = atom->nlocal;
+  tagint *tag = atom->tag;
+  int *num_bond = atom->num_bond;
+  tagint **bond_atom = atom->bond_atom;
+  int **bond_type = atom->bond_type;
+  int *mask = atom->mask;
+
+  m = 0;
+  atoms2bond.clear();
+  for (atom1 = 0; atom1 < nlocal; atom1++) {
+    if (!(mask[atom1] & cgroupbit)) continue;
+    nb = num_bond[atom1];
+    for (i = 0; i < nb; i++) {
+      btype = bond_type[atom1][i];
+      atom2 = atom->map(bond_atom[atom1][i]);
+      if (atom2 < 0 || !(mask[atom2] & cgroupbit)) continue;
+      if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
+      if (btype == 0) continue;
+      aset = {tag[atom1], tag[atom2]};
+      atoms2bond.insert(std::make_pair(aset,m++));
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
 return handedness (1 or -1) of a chiral center, given ordered set of coordinates
 ------------------------------------------------------------------------- */
diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h
index 628a71649b..41908e40d4 100644
--- a/src/REACTION/fix_bond_react.h
+++ b/src/REACTION/fix_bond_react.h
@@ -27,6 +27,9 @@ FixStyle(bond/react,FixBondReact);
 #include "fix.h"
 #include "compute.h"
 
+#include <map>
+#include <set>
+
 namespace LAMMPS_NS {
 
 class FixBondReact : public Fix {
@@ -79,9 +82,10 @@ class FixBondReact : public Fix {
   char **constraintstr;
   int nrxnfunction;
   std::vector<std::string> rxnfunclist;     // lists current special rxn function
-  std::vector<int> peratomflag;     // 1 if special rxn function uses per-atom variable (vs. per-bond)
+  std::vector<int> peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond)
+  int atoms2bondflag;           // 1 if atoms2bond map has been populated on this timestep
   int narrhenius;
-  int **var_flag, **var_id;    // for keyword values with variable inputs
+  int **var_flag, **var_id;     // for keyword values with variable inputs
   int status;
   int *groupbits;
 
@@ -191,6 +195,7 @@ class FixBondReact : public Fix {
   double custom_constraint(const std::string &);    // evaulate expression for custom constraint
   double rxnfunction(const std::string &, const std::string &,
                      const std::string &);    // eval rxn_sum and rxn_ave
+  void get_atoms2bond(int);
   int get_chirality(double[12]);              // get handedness given an ordered set of coordinates
 
   void open(char *);
@@ -231,6 +236,7 @@ class FixBondReact : public Fix {
   int nvvec;
   double **vvec;    // per-atom vector to store custom constraint atom-style variable values
   Compute *cperbond;    // pointer to 'compute bond/local' used by custom constraint ('rxnbond' function)
+  std::map<std::set<tagint>, int> atoms2bond;    // maps atom pair to index of local bond array
   std::vector<std::vector<Constraint>> constraints;
 
   // DEBUG