adjust a few more fix styles for fix bond/history being optional

This commit is contained in:
Axel Kohlmeyer
2024-01-13 04:50:28 -05:00
parent a6b0c349d8
commit 497c48bd80
3 changed files with 42 additions and 3 deletions

View File

@ -25,7 +25,9 @@
#include "respa.h" #include "respa.h"
#include "update.h" #include "update.h"
#if defined(LMP_BPM)
#include "fix_bond_history.h" #include "fix_bond_history.h"
#endif
#include <cstring> #include <cstring>
@ -261,9 +263,11 @@ void FixBondBreak::post_integrate()
commflag = 1; commflag = 1;
comm->forward_comm(this,2); comm->forward_comm(this,2);
#if defined(LMP_BPM)
// find instances of bond history to delete data // find instances of bond history to delete data
auto histories = modify->get_fix_by_style("BOND_HISTORY"); auto histories = modify->get_fix_by_style("BOND_HISTORY");
int n_histories = histories.size(); int n_histories = histories.size();
#endif
// break bonds // break bonds
// if both atoms list each other as winning bond partner // if both atoms list each other as winning bond partner
@ -299,13 +303,17 @@ void FixBondBreak::post_integrate()
for (k = m; k < num_bond[i]-1; k++) { for (k = m; k < num_bond[i]-1; k++) {
bond_atom[i][k] = bond_atom[i][k+1]; bond_atom[i][k] = bond_atom[i][k+1];
bond_type[i][k] = bond_type[i][k+1]; bond_type[i][k] = bond_type[i][k+1];
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->shift_history(i,k,k+1); dynamic_cast<FixBondHistory *>(ihistory)->shift_history(i,k,k+1);
#endif
} }
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i,num_bond[i]-1); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i,num_bond[i]-1);
#endif
num_bond[i]--; num_bond[i]--;
break; break;
} }

View File

@ -22,7 +22,6 @@
#include "compute.h" #include "compute.h"
#include "domain.h" #include "domain.h"
#include "error.h" #include "error.h"
#include "fix_bond_history.h"
#include "force.h" #include "force.h"
#include "memory.h" #include "memory.h"
#include "modify.h" #include "modify.h"
@ -32,6 +31,10 @@
#include "random_mars.h" #include "random_mars.h"
#include "update.h" #include "update.h"
#if defined(LMP_BPM)
#include "fix_bond_history.h"
#endif
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
@ -467,9 +470,11 @@ void FixBondSwap::post_integrate()
if (!accept) return; if (!accept) return;
naccept++; naccept++;
#if defined(LMP_BPM)
// find instances of bond/history to reset history // find instances of bond/history to reset history
auto histories = modify->get_fix_by_style("BOND_HISTORY"); auto histories = modify->get_fix_by_style("BOND_HISTORY");
int n_histories = histories.size(); int n_histories = histories.size();
#endif
// change bond partners of affected atoms // change bond partners of affected atoms
// on atom i: bond i-inext changes to i-jnext // on atom i: bond i-inext changes to i-jnext
@ -479,30 +484,38 @@ void FixBondSwap::post_integrate()
for (ibond = 0; ibond < num_bond[i]; ibond++) for (ibond = 0; ibond < num_bond[i]; ibond++)
if (bond_atom[i][ibond] == tag[inext]) { if (bond_atom[i][ibond] == tag[inext]) {
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i,ibond); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i,ibond);
#endif
bond_atom[i][ibond] = tag[jnext]; bond_atom[i][ibond] = tag[jnext];
} }
for (jbond = 0; jbond < num_bond[j]; jbond++) for (jbond = 0; jbond < num_bond[j]; jbond++)
if (bond_atom[j][jbond] == tag[jnext]) { if (bond_atom[j][jbond] == tag[jnext]) {
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(j,jbond); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(j,jbond);
#endif
bond_atom[j][jbond] = tag[inext]; bond_atom[j][jbond] = tag[inext];
} }
for (ibond = 0; ibond < num_bond[inext]; ibond++) for (ibond = 0; ibond < num_bond[inext]; ibond++)
if (bond_atom[inext][ibond] == tag[i]) { if (bond_atom[inext][ibond] == tag[i]) {
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(inext,ibond); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(inext,ibond);
#endif
bond_atom[inext][ibond] = tag[j]; bond_atom[inext][ibond] = tag[j];
} }
for (jbond = 0; jbond < num_bond[jnext]; jbond++) for (jbond = 0; jbond < num_bond[jnext]; jbond++)
if (bond_atom[jnext][jbond] == tag[j]) { if (bond_atom[jnext][jbond] == tag[j]) {
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(jnext,jbond); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(jnext,jbond);
#endif
bond_atom[jnext][jbond] = tag[i]; bond_atom[jnext][jbond] = tag[i];
} }

View File

@ -25,7 +25,6 @@ Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
#include "compute.h" #include "compute.h"
#include "domain.h" #include "domain.h"
#include "error.h" #include "error.h"
#include "fix_bond_history.h"
#include "force.h" #include "force.h"
#include "group.h" #include "group.h"
#include "input.h" #include "input.h"
@ -43,6 +42,10 @@ Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
#include "update.h" #include "update.h"
#include "variable.h" #include "variable.h"
#if defined(LMP_BPM)
#include "fix_bond_history.h"
#endif
#include "superpose3d.h" #include "superpose3d.h"
#include <cctype> #include <cctype>
@ -3256,10 +3259,11 @@ void FixBondReact::update_everything()
// next let's update bond info // next let's update bond info
// cool thing is, newton_bond issues are already taken care of in templates // 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 // same with class2 improper issues, which is why this fix started in the first place
#if defined(LMP_BPM)
// also need to find any instances of bond history to update histories // also need to find any instances of bond history to update histories
auto histories = modify->get_fix_by_style("BOND_HISTORY"); auto histories = modify->get_fix_by_style("BOND_HISTORY");
int n_histories = histories.size(); int n_histories = histories.size();
#endif
for (int i = 0; i < update_num_mega; i++) { for (int i = 0; i < update_num_mega; i++) {
rxnID = update_mega_glove[0][i]; rxnID = update_mega_glove[0][i];
twomol = atom->molecules[reacted_mol[rxnID]]; twomol = atom->molecules[reacted_mol[rxnID]];
@ -3269,6 +3273,7 @@ void FixBondReact::update_everything()
if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) { 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) { if (landlocked_atoms[j][rxnID] == 1) {
delta_bonds -= num_bond[atom->map(update_mega_glove[jj+1][i])]; delta_bonds -= num_bond[atom->map(update_mega_glove[jj+1][i])];
#if defined(LMP_BPM)
// If deleting all bonds, first cache then remove all histories // If deleting all bonds, first cache then remove all histories
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) { for (auto &ihistory: histories) {
@ -3277,6 +3282,7 @@ void FixBondReact::update_everything()
for (int n = 0; n < num_bond[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++)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(atom->map(update_mega_glove[jj+1][i]), 0); dynamic_cast<FixBondHistory *>(ihistory)->delete_history(atom->map(update_mega_glove[jj+1][i]), 0);
} }
#endif
num_bond[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) { if (landlocked_atoms[j][rxnID] == 0) {
@ -3284,21 +3290,27 @@ void FixBondReact::update_everything()
for (int n = 0; n < twomol->natoms; n++) { for (int n = 0; n < twomol->natoms; n++) {
int nn = equivalences[n][1][rxnID]-1; 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) { 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) {
#if defined(LMP_BPM)
// Cache history information, shift history, then delete final element // Cache history information, shift history, then delete final element
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->cache_history(atom->map(update_mega_glove[jj+1][i]), p); dynamic_cast<FixBondHistory *>(ihistory)->cache_history(atom->map(update_mega_glove[jj+1][i]), p);
#endif
for (int m = p; m < num_bond[atom->map(update_mega_glove[jj+1][i])]-1; m++) { 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_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]; bond_atom[atom->map(update_mega_glove[jj+1][i])][m] = bond_atom[atom->map(update_mega_glove[jj+1][i])][m+1];
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->shift_history(atom->map(update_mega_glove[jj+1][i]),m,m+1); dynamic_cast<FixBondHistory *>(ihistory)->shift_history(atom->map(update_mega_glove[jj+1][i]),m,m+1);
#endif
} }
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->delete_history(atom->map(update_mega_glove[jj+1][i]), dynamic_cast<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])]-1);
#endif
num_bond[atom->map(update_mega_glove[jj+1][i])]--; num_bond[atom->map(update_mega_glove[jj+1][i])]--;
delta_bonds--; delta_bonds--;
} }
@ -3317,10 +3329,12 @@ void FixBondReact::update_everything()
for (int p = 0; p < twomol->num_bond[j]; p++) { 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_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]; bond_atom[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
#if defined(LMP_BPM)
// Check cached history data to see if bond regenerated // Check cached history data to see if bond regenerated
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), p); dynamic_cast<FixBondHistory *>(ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), p);
#endif
} }
} }
if (landlocked_atoms[j][rxnID] == 0) { if (landlocked_atoms[j][rxnID] == 0) {
@ -3329,10 +3343,12 @@ void FixBondReact::update_everything()
insert_num = num_bond[atom->map(update_mega_glove[jj+1][i])]; 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_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]; 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];
#if defined(LMP_BPM)
// Check cached history data to see if bond regenerated // Check cached history data to see if bond regenerated
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), insert_num); dynamic_cast<FixBondHistory *>(ihistory)->check_cache(atom->map(update_mega_glove[jj+1][i]), insert_num);
#endif
num_bond[atom->map(update_mega_glove[jj+1][i])]++; 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) if (num_bond[atom->map(update_mega_glove[jj+1][i])] > atom->bond_per_atom)
error->one(FLERR,"Fix bond/react topology/atom exceed system topology/atom"); error->one(FLERR,"Fix bond/react topology/atom exceed system topology/atom");
@ -3344,9 +3360,11 @@ void FixBondReact::update_everything()
} }
} }
#if defined(LMP_BPM)
if (n_histories > 0) if (n_histories > 0)
for (auto &ihistory: histories) for (auto &ihistory: histories)
dynamic_cast<FixBondHistory *>(ihistory)->clear_cache(); dynamic_cast<FixBondHistory *>(ihistory)->clear_cache();
#endif
// Angles! First let's delete all angle info: // Angles! First let's delete all angle info:
if (force->angle && twomol->angleflag) { if (force->angle && twomol->angleflag) {