diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 04521cd0d8..6cdb5d7ee2 100644 --- a/src/BPM/bond_bpm.cpp +++ b/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; } diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 7973add89d..53d37df5a1 100644 --- a/src/MC/fix_bond_break.cpp +++ b/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; } diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 5f7ffcbbe5..a3a09c180d 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/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 diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 31a09567ca..08eff6dc48 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/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; diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 682a4476b5..ca933f3c97 100644 --- a/src/delete_bonds.cpp +++ b/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) - ((FixBondHistory *) ihistory)->delete_bond(i,m); + for (auto &ihistory: histories) { + ((FixBondHistory *) ihistory)->shift_history(i,m,n-1); + ((FixBondHistory *) ihistory)->delete_history(i,n-1); + } atom->num_bond[i]--; } else m++; } else m++; diff --git a/src/dump.cpp b/src/dump.cpp index 4c02aa3070..7eca4a938c 100644 --- a/src/dump.cpp +++ b/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) { diff --git a/src/fix_bond_history.cpp b/src/fix_bond_history.cpp index dfadd4201d..1f6ac6fe48 100644 --- a/src/fix_bond_history.cpp +++ b/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 -#include +#include +#include +#include +#include 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 key = std::make_pair(min_tag, max_tag); + + // Copy data to vector + double **stored = atom->darray[index]; + std::vector 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 key = std::make_pair(min_tag, max_tag); + + // Check if it exists, if so, copy data + double **stored = atom->darray[index]; + std::vector 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(); +} \ No newline at end of file diff --git a/src/fix_bond_history.h b/src/fix_bond_history.h index 970b1429d6..b1651a4dbe 100644 --- a/src/fix_bond_history.h +++ b/src/fix_bond_history.h @@ -22,6 +22,10 @@ FixStyle(BOND_HISTORY,FixBondHistory) #include "fix.h" +#include +#include +#include + 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::vector > cached_histories; + double **bondstore; int stored_flag;