From 3faecc4d2817d25be4c4015d88a0b315ea632e2c Mon Sep 17 00:00:00 2001 From: jrgissing Date: Tue, 30 Oct 2018 22:11:52 -0600 Subject: [PATCH] add option to update all atoms' atomic charges option to update all atomic charges, even when edge atoms are defined --- doc/src/fix_bond_react.txt | 25 ++++++++++++++++++++++--- src/USER-MISC/fix_bond_react.cpp | 12 ++++++++++-- src/USER-MISC/fix_bond_react.h | 1 + 3 files changed, 33 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 8e7cb1bdae..f806da3324 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -41,7 +41,10 @@ react = mandatory argument indicating new reaction specification :l fraction = initiate reaction with this probability if otherwise eligible seed = random number seed (positive integer) {stabilize_steps} value = timesteps - timesteps = number of timesteps to apply internally created nve/limit.html :pre + timesteps = number of timesteps to apply internally created nve/limit.html + {update_edges} value = {none} or {charges} :l + none = do not update topology near the edges of reaction templates + charges = update atomic charges of all atoms in reaction templates :pre :ule [Examples:] @@ -155,7 +158,17 @@ Some atoms in the pre-reacted template that are not reacting may have missing topology with respect to the simulation. For example, the pre-reacted template may contain an atom that would connect to the rest of a long polymer chain. These are referred to as edge atoms, and -are also specified in the map file. +are also specified in the map file. When the pre-reaction template +contains edge atoms, not all atoms, bonds, charges, etc. specified in +the reaction templates will be updated. Specifically, topology that +involves only atoms that are 'too near' to template edges will not be +updated. The definition of 'too near the edge' depends on which +interactions are defined in the simulation. If the simulation has +defined dihedrals, atoms within two bonds of edge atoms are considered +'too near the edge.' If the simulation defines angles, but not +dihedrals, atoms within one bond of edge atoms are considered 'too +near the edge.' If just bonds are defined, only edge atoms are +considered 'too near the edge.' Note that some care must be taken when a building a molecule template for a given simulation. All atom types in the pre-reacted template @@ -255,6 +268,12 @@ The {stabilize_steps} keyword allows for the specification of how many timesteps a reaction site is stabilized before being returned to the overall system thermostat. +The {update_edges} keyword can increase the number of atoms whose +atomic charges are updated, when the pre-reaction template contains +edge atoms. When the value is set to 'charges,' all atoms' atomic +charges are updated to those specified by the post-reaction template, +including atoms near the edge of reaction templates. + In order to produce the most physical behavior, this 'reaction site equilibration time' should be tuned to be as small as possible while retaining stability for a given system or reaction step. After a @@ -323,7 +342,7 @@ bond/break"_fix_bond_break.html, "fix bond/swap"_fix_bond_swap.html, [Default:] -The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60 +The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, update_edges = none :line diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index b30f1b36c6..705bb8ff70 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -163,6 +163,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(seed,nreacts,"bond/react:seed"); memory->create(limit_duration,nreacts,"bond/react:limit_duration"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); + memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag"); memory->create(iatomtype,nreacts,"bond/react:iatomtype"); memory->create(jatomtype,nreacts,"bond/react:jatomtype"); memory->create(ibonding,nreacts,"bond/react:ibonding"); @@ -178,6 +179,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fraction[i] = 1; seed[i] = 12345; stabilize_steps_flag[i] = 0; + update_edges_flag[i] = 0; // set default limit duration to 60 timesteps limit_duration[i] = 60; reaction_count[i] = 0; @@ -249,6 +251,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : limit_duration[rxn] = force->numeric(FLERR,arg[iarg+1]); stabilize_steps_flag[rxn] = 1; iarg += 2; + } else if (strcmp(arg[iarg],"update_edges") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'update_edges' has too few arguments"); + if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1; + iarg += 2; } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } } @@ -379,6 +386,7 @@ FixBondReact::~FixBondReact() memory->destroy(seed); memory->destroy(limit_duration); memory->destroy(stabilize_steps_flag); + memory->destroy(update_edges_flag); memory->destroy(iatomtype); memory->destroy(jatomtype); @@ -2022,13 +2030,13 @@ void FixBondReact::update_everything() } // update charges and types of landlocked atoms - // here, add check for charge instead of requiring it for (int i = 0; i < update_num_mega; i++) { rxnID = update_mega_glove[0][i]; twomol = atom->molecules[reacted_mol[rxnID]]; for (int j = 0; j < twomol->natoms; j++) { int jj = equivalences[j][1][rxnID]-1; - if (landlocked_atoms[j][rxnID] == 1 && atom->map(update_mega_glove[jj+1][i]) >= 0 && + if ((landlocked_atoms[j][rxnID] == 1 || update_edges_flag[rxnID] == 1) && + atom->map(update_mega_glove[jj+1][i]) >= 0 && atom->map(update_mega_glove[jj+1][i]) < nlocal) { type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; if (twomol->qflag && atom->q_flag) { diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index 08756f8131..b960ec8e73 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -58,6 +58,7 @@ class FixBondReact : public Fix { tagint lastcheck; int stabilization_flag; int *stabilize_steps_flag; + int *update_edges_flag; int status; int *groupbits;