Revamping public methods for history, adding support for bond react

src/BPM/bond_bpm.cpp

@@ -334,7 +334,8 @@ void BondBPM::process_broken(int i, int j)
         n = num_bond[i];
         bond_type[i][m] = bond_type[i][n-1];
         bond_atom[i][m] = bond_atom[i][n-1];
-        fix_bond_history->delete_bond(i, m);
+        fix_bond_history->shift_history(i, m, n-1);
+        fix_bond_history->delete_history(i, n-1);
         num_bond[i]--;
         break;
       }
@@ -348,7 +349,8 @@ void BondBPM::process_broken(int i, int j)
         n = num_bond[j];
         bond_type[j][m] = bond_type[j][n-1];
         bond_atom[j][m] = bond_atom[j][n-1];
-        fix_bond_history->delete_bond(j, m);
+        fix_bond_history->shift_history(j, m, n-1);
+        fix_bond_history->delete_history(j, n-1);
         num_bond[j]--;
         break;
       }

src/MC/fix_bond_break.cpp

@@ -303,11 +303,11 @@ void FixBondBreak::post_integrate()
           bond_type[i][k] = bond_type[i][k+1];
           if (n_histories > 0)
             for (auto &ihistory: histories)
-              ((FixBondHistory *) ihistory)->shift_bond(i,k,k+1);
+              ((FixBondHistory *) ihistory)->shift_history(i,k,k+1);
         }
         if (n_histories > 0)
           for (auto &ihistory: histories)
-            ((FixBondHistory *) ihistory)->delete_bond(i,num_bond[i]-1);
+            ((FixBondHistory *) ihistory)->delete_history(i,num_bond[i]-1);
         num_bond[i]--;
         break;
       }

src/MC/fix_bond_swap.cpp

@@ -22,6 +22,7 @@
 #include "compute.h"
 #include "domain.h"
 #include "error.h"
+#include "fix_bond_history.h"
 #include "force.h"
 #include "memory.h"
 #include "modify.h"
@@ -450,6 +451,10 @@ void FixBondSwap::post_integrate()
   if (!accept) return;
   naccept++;
 
+  // find instances of bond/history to reset history
+  auto histories = modify->get_fix_by_style("BOND_HISTORY");
+  int n_histories = histories.size();
+
   // change bond partners of affected atoms
   // on atom i: bond i-inext changes to i-jnext
   // on atom j: bond j-jnext changes to j-inext
@@ -457,13 +462,33 @@ void FixBondSwap::post_integrate()
   // on atom jnext: bond jnext-j changes to jnext-i
 
   for (ibond = 0; ibond < num_bond[i]; ibond++)
-    if (bond_atom[i][ibond] == tag[inext]) bond_atom[i][ibond] = tag[jnext];
+    if (bond_atom[i][ibond] == tag[inext]) {
+      if (n_histories > 0)
+        for (auto &ihistory: histories)
+          ((FixBondHistory *) ihistory)->delete_history(i,ibond);
+      bond_atom[i][ibond] = tag[jnext];
+    }
   for (jbond = 0; jbond < num_bond[j]; jbond++)
-    if (bond_atom[j][jbond] == tag[jnext]) bond_atom[j][jbond] = tag[inext];
+    if (bond_atom[j][jbond] == tag[jnext]) {
+      if (n_histories > 0)
+        for (auto &ihistory: histories)
+          ((FixBondHistory *) ihistory)->delete_history(j,jbond);
+      bond_atom[j][jbond] = tag[inext];
+    }
   for (ibond = 0; ibond < num_bond[inext]; ibond++)
-    if (bond_atom[inext][ibond] == tag[i]) bond_atom[inext][ibond] = tag[j];
+    if (bond_atom[inext][ibond] == tag[i]) {
+      if (n_histories > 0)
+        for (auto &ihistory: histories)
+          ((FixBondHistory *) ihistory)->delete_history(inext,ibond);
+      bond_atom[inext][ibond] = tag[j];
+    }
   for (jbond = 0; jbond < num_bond[jnext]; jbond++)
-    if (bond_atom[jnext][jbond] == tag[j]) bond_atom[jnext][jbond] = tag[i];
+    if (bond_atom[jnext][jbond] == tag[j]) {
+      if (n_histories > 0)
+        for (auto &ihistory: histories)
+          ((FixBondHistory *) ihistory)->delete_history(jnext,jbond);
+      bond_atom[jnext][jbond] = tag[i];
+    }
 
   // set global tags of 4 atoms in bonds
 

src/REACTION/fix_bond_react.cpp

@@ -24,6 +24,7 @@ Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "fix_bond_history.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -3091,6 +3092,10 @@ void FixBondReact::update_everything()
     // next let's update bond info
     // cool thing is, newton_bond issues are already taken care of in templates
     // same with class2 improper issues, which is why this fix started in the first place
+    // also need to find any instances of bond history to update histories
+    auto histories = modify->get_fix_by_style("BOND_HISTORY");
+    int n_histories = histories.size();
+
     for (int i = 0; i < update_num_mega; i++) {
       rxnID = update_mega_glove[0][i];
       twomol = atom->molecules[reacted_mol[rxnID]];
@@ -3100,6 +3105,14 @@ void FixBondReact::update_everything()
         if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
           if (landlocked_atoms[j][rxnID] == 1) {
             delta_bonds -= num_bond[atom->map(update_mega_glove[jj+1][i])];
+            // If deleting all bonds, first cache then remove all histories
+            if (n_histories > 0)
+              for (auto &ihistory: histories) {
+                for (int n = 0; n < num_bond[atom->map(update_mega_glove[jj+1][i])]; n++)
+                  ((FixBondHistory *) ihistory)->cache_history(atom->map(update_mega_glove[jj+1][i]), n);
+                for (int n = 0; n < num_bond[atom->map(update_mega_glove[jj+1][i])]; n++)
+                  ((FixBondHistory *) ihistory)->delete_history(atom->map(update_mega_glove[jj+1][i]), 0);
+              }
             num_bond[atom->map(update_mega_glove[jj+1][i])] = 0;
           }
           if (landlocked_atoms[j][rxnID] == 0) {
@@ -3107,10 +3120,21 @@ void FixBondReact::update_everything()
               for (int n = 0; n < twomol->natoms; n++) {
                 int nn = equivalences[n][1][rxnID]-1;
                 if (n!=j && bond_atom[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] && landlocked_atoms[n][rxnID] == 1) {
+                  // Cache history information, shift history, then delete final element
+                  if (n_histories > 0)
+                    for (auto &ihistory: histories)
+                      ((FixBondHistory *) ihistory)->cache_history(atom->map(update_mega_glove[jj+1][i]), p);
                   for (int m = p; m < num_bond[atom->map(update_mega_glove[jj+1][i])]-1; m++) {
                     bond_type[atom->map(update_mega_glove[jj+1][i])][m] = bond_type[atom->map(update_mega_glove[jj+1][i])][m+1];
                     bond_atom[atom->map(update_mega_glove[jj+1][i])][m] = bond_atom[atom->map(update_mega_glove[jj+1][i])][m+1];
+                    if (n_histories > 0)
+                      for (auto &ihistory: histories)
+                        ((FixBondHistory *) ihistory)->shift_history(atom->map(update_mega_glove[jj+1][i]),m,m+1);
                   }
+                  if (n_histories > 0)
+                    for (auto &ihistory: histories)
+                      ((FixBondHistory *) ihistory)->delete_history(atom->map(update_mega_glove[jj+1][i]),
+                                                                 num_bond[atom->map(update_mega_glove[jj+1][i])]-1);
                   num_bond[atom->map(update_mega_glove[jj+1][i])]--;
                   delta_bonds--;
                 }
@@ -3129,6 +3153,10 @@ void FixBondReact::update_everything()
             for (int p = 0; p < twomol->num_bond[j]; p++) {
               bond_type[atom->map(update_mega_glove[jj+1][i])][p] = twomol->bond_type[j][p];
               bond_atom[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
+              // Check cached history data to see if bond regenerated
+              if (n_histories > 0)
+                for (auto &ihistory: histories)
+                  ((FixBondHistory *) ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), p);
             }
           }
           if (landlocked_atoms[j][rxnID] == 0) {
@@ -3137,6 +3165,10 @@ void FixBondReact::update_everything()
                 insert_num = num_bond[atom->map(update_mega_glove[jj+1][i])];
                 bond_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->bond_type[j][p];
                 bond_atom[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
+                // Check cached history data to see if bond regenerated
+                if (n_histories > 0)
+                  for (auto &ihistory: histories)
+                    ((FixBondHistory *) ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), insert_num);
                 num_bond[atom->map(update_mega_glove[jj+1][i])]++;
                 if (num_bond[atom->map(update_mega_glove[jj+1][i])] > atom->bond_per_atom)
                   error->one(FLERR,"Bond/react topology/atom exceed system topology/atom");
@@ -3148,6 +3180,10 @@ void FixBondReact::update_everything()
       }
     }
 
+    if (n_histories > 0)
+      for (auto &ihistory: histories)
+        ((FixBondHistory *) ihistory)->clear_cache();
+
     // Angles! First let's delete all angle info:
     if (force->angle && twomol->angleflag) {
       int *num_angle = atom->num_angle;

src/delete_bonds.cpp

@@ -338,8 +338,10 @@ void DeleteBonds::command(int narg, char **arg)
               atom->bond_type[i][m] = atom->bond_type[i][n-1];
               atom->bond_atom[i][m] = atom->bond_atom[i][n-1];
               if (n_histories > 0)
-                for (auto &ihistory: histories)
+                for (auto &ihistory: histories) {
-                  ((FixBondHistory *) ihistory)->delete_bond(i,m);
+                  ((FixBondHistory *) ihistory)->shift_history(i,m,n-1);
+                  ((FixBondHistory *) ihistory)->delete_history(i,n-1);
+                }
               atom->num_bond[i]--;
             } else m++;
           } else m++;

src/dump.cpp

@@ -920,7 +920,7 @@ void Dump::modify_params(int narg, char **arg)
 
     } else if (strcmp(arg[iarg],"header") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
-      header_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      write_header_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"every") == 0) {

src/fix_bond_history.cpp

@@ -17,15 +17,15 @@
 #include "atom.h"
 #include "comm.h"
 #include "error.h"
-#include "force.h"
 #include "group.h"
 #include "memory.h"
 #include "modify.h"
 #include "neighbor.h"
-#include "string.h"
 
-#include <cstring>
+#include <map>
-#include <mpi.h>
+#include <string>
+#include <utility>
+#include <vector>
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -284,28 +284,80 @@ void FixBondHistory::set_arrays(int i)
 }
 
 /* ----------------------------------------------------------------------
-   Remove all data for row by compressing data - moving last element
+   Delete bond by zeroing data
 ------------------------------------------------------------------------- */
 
-void FixBondHistory::delete_bond(int i, int m)
+void FixBondHistory::delete_history(int i, int m)
 {
   double **stored = atom->darray[index];
-  int n = atom->num_bond[i];
-  if (m != n-1) shift_bond(i, m, n-1);
-
   for (int idata = 0; idata < ndata; idata ++)
-    stored[i][(n-1)*ndata+idata] = 0.0;
+    stored[i][m*ndata+idata] = 0.0;
 }
 
 /* ----------------------------------------------------------------------
    Shift bond data to a new location
-   Used by some fixes that delete bonds which don't move last element
 ------------------------------------------------------------------------- */
 
-void FixBondHistory::shift_bond(int i, int m, int k)
+void FixBondHistory::shift_history(int i, int m, int k)
 {
+  if (m == k) return;
+
   double **stored = atom->darray[index];
   int n = atom->num_bond[i];
   for (int idata = 0; idata < ndata; idata ++)
     stored[i][m*ndata+idata] = stored[i][k*ndata+idata];
 }
+
+/* ----------------------------------------------------------------------
+   Temporarily caches history for a deleted bond which
+   could be recreated before the cache is emptied
+   NOTE: the cache methods still need to be tested, need an example first
+------------------------------------------------------------------------- */
+
+void FixBondHistory::cache_history(int i, int m)
+{
+  // Order tags to create a unique key pair
+  tagint max_tag = MAX(atom->tag[i], atom->bond_atom[i][m]);
+  tagint min_tag = MIN(atom->tag[i], atom->bond_atom[i][m]);
+  std::pair<tagint, tagint> key = std::make_pair(min_tag, max_tag);
+
+  // Copy data to vector
+  double **stored = atom->darray[index];
+  std::vector <double> data;
+  for (int idata = 0; idata < ndata; idata ++)
+    data.push_back(stored[i][m*ndata+idata]);
+
+  // Add data to cache
+  cached_histories.insert(std::make_pair(key,data));
+}
+
+/* ----------------------------------------------------------------------
+   Checks to see if a newly created bond has cached history
+------------------------------------------------------------------------- */
+
+void FixBondHistory::check_cache(int i, int m)
+{
+  // Order tags to create a unique key pair
+  tagint max_tag = MAX(atom->tag[i], atom->bond_atom[i][m]);
+  tagint min_tag = MIN(atom->tag[i], atom->bond_atom[i][m]);
+  std::pair<tagint, tagint> key = std::make_pair(min_tag, max_tag);
+
+  // Check if it exists, if so, copy data
+  double **stored = atom->darray[index];
+  std::vector <double> data;
+  auto pos = cached_histories.find(key);
+  if (pos != cached_histories.end()) {
+    data = pos->second;
+    for (int idata = 0; idata < ndata; idata ++)
+      stored[i][m*ndata+idata] = data[idata];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   Delete saved memory
+------------------------------------------------------------------------- */
+
+void FixBondHistory::clear_cache()
+{
+  cached_histories.clear();
+}

src/fix_bond_history.h

@@ -22,6 +22,10 @@ FixStyle(BOND_HISTORY,FixBondHistory)
 
 #include "fix.h"
 
+#include <map>
+#include <utility>
+#include <vector>
+
 namespace LAMMPS_NS {
 
 class FixBondHistory : public Fix {
@@ -41,8 +45,19 @@ class FixBondHistory : public Fix {
 
   void update_atom_value(int, int, int, double);
   double get_atom_value(int, int, int);
-  void delete_bond(int,int);
+
-  void shift_bond(int,int,int);
+  // methods to reorder/delete elements of atom->bond_atom
+  void delete_history(int,int);
+  void shift_history(int,int,int);
+  void cache_history(int,int);
+  void check_cache(int,int);
+  void clear_cache();
+
+  // if data is temporarily stored while the bond_atom array
+  // is being reordered, use map of vectors with pairs for keys
+  // to enable quick look up
+  std::map<std::pair<tagint, tagint>, std::vector <double>> cached_histories;
+
   double **bondstore;
   int stored_flag;