From 1b1cb5568d1f475ce0316a6d0487d2cf83e875cf Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 14:48:47 -0400 Subject: [PATCH 01/35] non-'bond/react ' changes --- src/compute.h | 5 +++++ src/compute_bond_local.cpp | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/src/compute.h b/src/compute.h index 90b88c7d1d..1d30510de8 100644 --- a/src/compute.h +++ b/src/compute.h @@ -16,6 +16,9 @@ #include "pointers.h" // IWYU pragma: export +#include +#include + namespace LAMMPS_NS { class Compute : protected Pointers { @@ -94,6 +97,8 @@ class Compute : protected Pointers { double dof; // degrees-of-freedom for temperature + std::map, int> atoms2bond; // maps atom pair to index of local bond array + int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index 11a8ed7148..cb9d278d5d 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -236,6 +236,7 @@ int ComputeBondLocal::compute_bonds(int flag) double mvv2e; double engpot,engtrans,engvib,engrot,fbond; double *ptr; + std::set aset; double **x = atom->x; double **v = atom->v; @@ -379,6 +380,11 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } + // populate map to access per-bond values from atom tags + + aset = {tag[atom1], tag[atom2]}; + atoms2bond.insert(std::make_pair(aset,m)); + if (nvalues == 1) ptr = &vlocal[m]; else ptr = alocal[m]; From 6b9b5daa0d0c27ff43f32470dff0b60d29e10739 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 15:50:27 -0400 Subject: [PATCH 02/35] bond/react: per-bond custom constraint --- src/REACTION/fix_bond_react.cpp | 62 +++++++++++++++++++++++++++------ src/REACTION/fix_bond_react.h | 5 ++- 2 files changed, 56 insertions(+), 11 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 6efb9630a2..4eecb68a88 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -137,10 +137,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : status = PROCEED; // reaction functions used by 'custom' constraint - nrxnfunction = 2; + nrxnfunction = 3; rxnfunclist.resize(nrxnfunction); + peratomflag.resize(nrxnfunction); rxnfunclist[0] = "rxnsum"; + peratomflag[0] = 1; rxnfunclist[1] = "rxnave"; + peratomflag[1] = 1; + rxnfunclist[2] = "rxnbond"; + peratomflag[2] = 0; nvvec = 0; ncustomvars = 0; vvec = nullptr; @@ -2105,7 +2110,7 @@ double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col) } /* ---------------------------------------------------------------------- -get per-atom variable names used by custom constraint +get per-atom variable names used by custom constraint ------------------------------------------------------------------------- */ void FixBondReact::customvarnames() @@ -2126,6 +2131,7 @@ void FixBondReact::customvarnames() // find next reaction special function occurrence pos1 = INT_MAX; for (int i = 0; i < nrxnfunction; i++) { + if (peratomflag[i] == 0) continue; pos = varstr.find(rxnfunclist[i],prev3+1); if (pos == std::string::npos) continue; if (pos < pos1) pos1 = pos; @@ -2251,12 +2257,55 @@ double FixBondReact::custom_constraint(const std::string& varstr) } /* ---------------------------------------------------------------------- -currently two 'rxn' functions: rxnsum and rxnave +currently three 'rxn' functions: rxnsum, rxnave, and rxnbond ------------------------------------------------------------------------- */ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& varid, const std::string& fragid) { + int ifrag = -1; + if (fragid != "all") { + ifrag = onemol->findfragment(fragid.c_str()); + if (ifrag < 0) error->one(FLERR,"Bond/react: Molecule fragment " + "in reaction special function does not exist"); + } + + // start with 'rxnbond' per-bond function + // for 'rxnbond', varid corresponds to 'compute bond/local' name, + // and fragid is a pre-reaction fragment containing the two atoms in the bond + if (rxnfunc == "rxnbond") { + int icompute,ibond,nsum; + double perbondval; + std::set aset; + std::string computeid = varid; + + if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " + "name should begin with 'c_'"); + computeid = computeid.substr(2); + icompute = modify->find_compute(computeid); + if (icompute < 0) error->one(FLERR,"Bond/react: Reaction special function compute name does not exist"); + cperbond = modify->compute[icompute]; + std::string compute_style = cperbond->style; + if (compute_style != "bond/local") error->one(FLERR,"Bond/react: Compute used by reaction " + "special function 'rxnbond' must be of style 'bond/local'"); + if (cperbond->size_local_cols > 0) error->one(FLERR,"Bond/react: 'Compute bond/local' used by reaction " + "special function 'rxnbond' must compute one value"); + + nsum = 0; + for (int i = 0; i < onemol->natoms; i++) { + if (onemol->fragmentmask[ifrag][i]) { + aset.insert(glove[i][1]); + nsum++; + } + } + if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " + "must contain exactly two atoms"); + + ibond = cperbond->atoms2bond.find(aset)->second; + perbondval = cperbond->vector_local[ibond]; + return perbondval; + } + int ivar = -1; for (int i = 0; i < ncustomvars; i++) { if (varid == customvarstrs[i]) { @@ -2270,13 +2319,6 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& error->one(FLERR,"Bond/react: Reaction special function variable " "name does not exist"); - int ifrag = -1; - if (fragid != "all") { - ifrag = onemol->findfragment(fragid.c_str()); - if (ifrag < 0) error->one(FLERR,"Bond/react: Molecule fragment " - "in reaction special function does not exist"); - } - int iatom; int nsum = 0; double sumvvec = 0; diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index f38fa1f136..c134da0bb4 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -25,6 +25,7 @@ FixStyle(bond/react,FixBondReact); #define LMP_FIX_BOND_REACT_H #include "fix.h" +#include "compute.h" namespace LAMMPS_NS { @@ -74,7 +75,8 @@ class FixBondReact : public Fix { int *nconstraints; char **constraintstr; int nrxnfunction; - std::vector rxnfunclist; + std::vector rxnfunclist; // lists current special rxn function + std::vector peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond) int narrhenius; int **var_flag, **var_id; // for keyword values with variable inputs int status; @@ -222,6 +224,7 @@ class FixBondReact : public Fix { std::vector customvarstrs; int nvvec; double **vvec; // per-atom vector to store variable constraint atom-style variable values + Compute *cperbond; // pointer to 'compute bond/local' used by custom constraint ('rxnbond' function) std::vector> constraints; // DEBUG From 7bbae300c8693ba40752ddd9c42dbe01c1cd290c Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:02:13 -0400 Subject: [PATCH 03/35] actually evaluate bond/local compute value (even when not printed on that timestep) --- src/REACTION/fix_bond_react.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 4eecb68a88..9f1338e641 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2301,6 +2301,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); + cperbond->compute_local(); ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From f0f39ced6f7aa5bf4919241d46bb22ae57fbb856 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:51:52 -0400 Subject: [PATCH 04/35] prevent multiple compute evaluations on a timestep --- src/REACTION/fix_bond_react.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 9f1338e641..94dd997198 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2301,7 +2301,9 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); - cperbond->compute_local(); + if (cperbond->invoked_local != lmp->update->ntimestep) + cperbond->compute_local(); + ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From 3f87b88d579f7866fd48de7157453039a8f1a1ba Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 29 May 2022 10:43:36 -0400 Subject: [PATCH 05/35] check that bond is actually in map --- src/REACTION/fix_bond_react.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 94dd997198..bd056e2667 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2278,6 +2278,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& double perbondval; std::set aset; std::string computeid = varid; + std::map,int>::iterator it; if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " "name should begin with 'c_'"); @@ -2304,7 +2305,10 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (cperbond->invoked_local != lmp->update->ntimestep) cperbond->compute_local(); - ibond = cperbond->atoms2bond.find(aset)->second; + it = cperbond->atoms2bond.find(aset); + if (it == cperbond->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]; return perbondval; } From dc378b8ffbfe9dfe4508ccffef5d7c700c704cab Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 14:48:47 -0400 Subject: [PATCH 06/35] non-'bond/react ' changes --- src/compute.h | 5 +++++ src/compute_bond_local.cpp | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/src/compute.h b/src/compute.h index 90b88c7d1d..1d30510de8 100644 --- a/src/compute.h +++ b/src/compute.h @@ -16,6 +16,9 @@ #include "pointers.h" // IWYU pragma: export +#include +#include + namespace LAMMPS_NS { class Compute : protected Pointers { @@ -94,6 +97,8 @@ class Compute : protected Pointers { double dof; // degrees-of-freedom for temperature + std::map, int> atoms2bond; // maps atom pair to index of local bond array + int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index 11a8ed7148..cb9d278d5d 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -236,6 +236,7 @@ int ComputeBondLocal::compute_bonds(int flag) double mvv2e; double engpot,engtrans,engvib,engrot,fbond; double *ptr; + std::set aset; double **x = atom->x; double **v = atom->v; @@ -379,6 +380,11 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } + // populate map to access per-bond values from atom tags + + aset = {tag[atom1], tag[atom2]}; + atoms2bond.insert(std::make_pair(aset,m)); + if (nvalues == 1) ptr = &vlocal[m]; else ptr = alocal[m]; From 930d0d756b2072372118355bd1cbc667bbab220e Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 15:50:27 -0400 Subject: [PATCH 07/35] bond/react: per-bond custom constraint --- src/REACTION/fix_bond_react.cpp | 62 +++++++++++++++++++++++++++------ src/REACTION/fix_bond_react.h | 5 ++- 2 files changed, 56 insertions(+), 11 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 0405971bdd..29a4ace892 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -137,10 +137,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : status = PROCEED; // reaction functions used by 'custom' constraint - nrxnfunction = 2; + nrxnfunction = 3; rxnfunclist.resize(nrxnfunction); + peratomflag.resize(nrxnfunction); rxnfunclist[0] = "rxnsum"; + peratomflag[0] = 1; rxnfunclist[1] = "rxnave"; + peratomflag[1] = 1; + rxnfunclist[2] = "rxnbond"; + peratomflag[2] = 0; nvvec = 0; ncustomvars = 0; vvec = nullptr; @@ -2112,7 +2117,7 @@ double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col) } /* ---------------------------------------------------------------------- -get per-atom variable names used by custom constraint +get per-atom variable names used by custom constraint ------------------------------------------------------------------------- */ void FixBondReact::customvarnames() @@ -2134,6 +2139,7 @@ void FixBondReact::customvarnames() // find next reaction special function occurrence pos1 = std::string::npos; for (int i = 0; i < nrxnfunction; i++) { + if (peratomflag[i] == 0) continue; pos = varstr.find(rxnfunclist[i],prev3+1); if (pos == std::string::npos) continue; if (pos < pos1) pos1 = pos; @@ -2260,12 +2266,55 @@ double FixBondReact::custom_constraint(const std::string& varstr) } /* ---------------------------------------------------------------------- -currently two 'rxn' functions: rxnsum and rxnave +currently three 'rxn' functions: rxnsum, rxnave, and rxnbond ------------------------------------------------------------------------- */ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& varid, const std::string& fragid) { + int ifrag = -1; + if (fragid != "all") { + ifrag = onemol->findfragment(fragid.c_str()); + if (ifrag < 0) error->one(FLERR,"Bond/react: Molecule fragment " + "in reaction special function does not exist"); + } + + // start with 'rxnbond' per-bond function + // for 'rxnbond', varid corresponds to 'compute bond/local' name, + // and fragid is a pre-reaction fragment containing the two atoms in the bond + if (rxnfunc == "rxnbond") { + int icompute,ibond,nsum; + double perbondval; + std::set aset; + std::string computeid = varid; + + if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " + "name should begin with 'c_'"); + computeid = computeid.substr(2); + icompute = modify->find_compute(computeid); + if (icompute < 0) error->one(FLERR,"Bond/react: Reaction special function compute name does not exist"); + cperbond = modify->compute[icompute]; + std::string compute_style = cperbond->style; + if (compute_style != "bond/local") error->one(FLERR,"Bond/react: Compute used by reaction " + "special function 'rxnbond' must be of style 'bond/local'"); + if (cperbond->size_local_cols > 0) error->one(FLERR,"Bond/react: 'Compute bond/local' used by reaction " + "special function 'rxnbond' must compute one value"); + + nsum = 0; + for (int i = 0; i < onemol->natoms; i++) { + if (onemol->fragmentmask[ifrag][i]) { + aset.insert(glove[i][1]); + nsum++; + } + } + if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " + "must contain exactly two atoms"); + + ibond = cperbond->atoms2bond.find(aset)->second; + perbondval = cperbond->vector_local[ibond]; + return perbondval; + } + int ivar = -1; for (int i = 0; i < ncustomvars; i++) { if (varid == customvarstrs[i]) { @@ -2279,13 +2328,6 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& error->one(FLERR,"Fix bond/react: Reaction special function variable " "name does not exist"); - int ifrag = -1; - if (fragid != "all") { - ifrag = onemol->findfragment(fragid.c_str()); - if (ifrag < 0) error->one(FLERR,"Fix bond/react: Molecule fragment " - "in reaction special function does not exist"); - } - int iatom; int nsum = 0; double sumvvec = 0; diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index f38fa1f136..c134da0bb4 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -25,6 +25,7 @@ FixStyle(bond/react,FixBondReact); #define LMP_FIX_BOND_REACT_H #include "fix.h" +#include "compute.h" namespace LAMMPS_NS { @@ -74,7 +75,8 @@ class FixBondReact : public Fix { int *nconstraints; char **constraintstr; int nrxnfunction; - std::vector rxnfunclist; + std::vector rxnfunclist; // lists current special rxn function + std::vector peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond) int narrhenius; int **var_flag, **var_id; // for keyword values with variable inputs int status; @@ -222,6 +224,7 @@ class FixBondReact : public Fix { std::vector customvarstrs; int nvvec; double **vvec; // per-atom vector to store variable constraint atom-style variable values + Compute *cperbond; // pointer to 'compute bond/local' used by custom constraint ('rxnbond' function) std::vector> constraints; // DEBUG From deb892d7cd25ba444fdcad4548bad4af5a69bea6 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:02:13 -0400 Subject: [PATCH 08/35] actually evaluate bond/local compute value (even when not printed on that timestep) --- src/REACTION/fix_bond_react.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 29a4ace892..50df6c3f87 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2310,6 +2310,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); + cperbond->compute_local(); ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From 32406aab06be2ec32ec4bd3ab3fe0fe315915287 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:51:52 -0400 Subject: [PATCH 09/35] prevent multiple compute evaluations on a timestep --- src/REACTION/fix_bond_react.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 50df6c3f87..e67d6432db 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2310,7 +2310,9 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); - cperbond->compute_local(); + if (cperbond->invoked_local != lmp->update->ntimestep) + cperbond->compute_local(); + ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From 07ab4dd5da64747ef2a51cf5a24c52c6865a377f Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 29 May 2022 10:43:36 -0400 Subject: [PATCH 10/35] check that bond is actually in map --- src/REACTION/fix_bond_react.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index e67d6432db..21287031d5 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2287,6 +2287,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& double perbondval; std::set aset; std::string computeid = varid; + std::map,int>::iterator it; if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " "name should begin with 'c_'"); @@ -2313,7 +2314,10 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (cperbond->invoked_local != lmp->update->ntimestep) cperbond->compute_local(); - ibond = cperbond->atoms2bond.find(aset)->second; + it = cperbond->atoms2bond.find(aset); + if (it == cperbond->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]; return perbondval; } From 1b59a83ff363e2272dca0e7ee9208127ca05d21e Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 29 Jul 2022 17:24:07 -0400 Subject: [PATCH 11/35] bond:react per-bond custom constraint docs --- doc/src/fix_bond_react.rst | 69 +++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 19 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index a48675dcf2..f9d7160e52 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -506,22 +506,32 @@ The constraint of type 'custom' has the following syntax: custom *varstring* -where 'custom' is the required keyword, and *varstring* is a -variable expression. The expression must be a valid equal-style -variable formula that can be read by the :doc:`variable ` command, +where 'custom' is the required keyword, and *varstring* is a variable +expression. The expression must be a valid equal-style variable +formula that can be read by the :doc:`variable ` command, after any special reaction functions are evaluated. If the resulting expression is zero, the reaction is prevented from occurring; -otherwise, it is permitted to occur. There are two special reaction -functions available, 'rxnsum' and 'rxnave'. These functions operate -over the atoms in a given reaction site, and have one mandatory -argument and one optional argument. The mandatory argument is the -identifier for an atom-style variable. The second, optional argument -is the name of a molecule fragment in the pre-reaction template, and -can be used to operate over a subset of atoms in the reaction site. -The 'rxnsum' function sums the atom-style variable over the reaction -site, while the 'rxnave' returns the average value. For example, a -constraint on the total potential energy of atoms involved in the -reaction can be imposed as follows: +otherwise, it is permitted to occur. There are three special reaction +functions available, 'rxnbond', 'rxnsum', and 'rxnave'. The 'rxnbond' +function allows per-bond values to be included in the variable strings +of the custom constraint. The 'rxnbond' function has two mandatory +arguments. The first argument is the ID of a previously defined +'compute bond/local' command. This 'compute bond/local' must compute +only one value, e.g. bond force. This value is returned by the +'rxnbond' function. The second argument is the name of a molecule +fragment in the pre-reaction template. The fragment must contain +exactly two atoms, corresponding to the atoms involved in the bond +whose value should be calculated. An example of a constraint that uses +the force experienced by a bond in the reaction site is provided +below. The 'rxnsum' and 'rxnave' functions operate over the atoms in a +given reaction site, and have one mandatory argument and one optional +argument. The mandatory argument is the identifier for an atom-style +variable. The second, optional argument is the name of a molecule +fragment in the pre-reaction template, and can be used to operate over +a subset of atoms in the reaction site. The 'rxnsum' function sums the +atom-style variable over the reaction site, while the 'rxnave' returns +the average value. For example, a constraint on the total potential +energy of atoms involved in the reaction can be imposed as follows: .. code-block:: LAMMPS @@ -533,11 +543,32 @@ reaction can be imposed as follows: custom "rxnsum(v_my_pe) > 100" # in Constraints section of map file The above example prevents the reaction from occurring unless the -total potential energy of the reaction site is above 100. The variable -expression can be interpreted as the probability of the reaction -occurring by using an inequality and the 'random(x,y,z)' function -available as an equal-style variable input, similar to the 'arrhenius' -constraint above. +total potential energy of the reaction site is above 100. As a second +example, this time using the 'rxnbond' function, consider a modified +Arrhenius constraint that depends on the bond force of a specific bond: + +.. code-block:: LAMMPS + + # in LAMMPS input script + + compute bondforce all bond/local force + + compute ke_atom all ke/atom + variable ke atom c_ke_atom + + variable E_a equal 100.0 # activation energy + variable natoms equal 10 # number of atoms in reaction site + + +.. code-block:: LAMMPS + + # in Constraints section of map file + + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3*rxnsum(v_ke)*v_natoms)) < random(0,1,12345)" + +By using an inequality and the 'random(x,y,z)' function, the left-hand +side can be interpreted as the probability of the reaction occurring, +similar to the 'arrhenius' constraint above. By default, all constraints must be satisfied for the reaction to occur. In other words, constraints are evaluated as a series of From 3e7ef39d64ecb54c044835232d2ff3caf56aaa01 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 29 Jul 2022 23:34:37 -0400 Subject: [PATCH 12/35] simplify, formatting --- doc/src/fix_bond_react.rst | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index f9d7160e52..3ed4e6ef26 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -522,16 +522,16 @@ only one value, e.g. bond force. This value is returned by the fragment in the pre-reaction template. The fragment must contain exactly two atoms, corresponding to the atoms involved in the bond whose value should be calculated. An example of a constraint that uses -the force experienced by a bond in the reaction site is provided -below. The 'rxnsum' and 'rxnave' functions operate over the atoms in a -given reaction site, and have one mandatory argument and one optional -argument. The mandatory argument is the identifier for an atom-style -variable. The second, optional argument is the name of a molecule -fragment in the pre-reaction template, and can be used to operate over -a subset of atoms in the reaction site. The 'rxnsum' function sums the -atom-style variable over the reaction site, while the 'rxnave' returns -the average value. For example, a constraint on the total potential -energy of atoms involved in the reaction can be imposed as follows: +the force experienced by a bond is provided below. The 'rxnsum' and +'rxnave' functions operate over the atoms in a given reaction site, +and have one mandatory argument and one optional argument. The +mandatory argument is the identifier for an atom-style variable. The +second, optional argument is the name of a molecule fragment in the +pre-reaction template, and can be used to operate over a subset of +atoms in the reaction site. The 'rxnsum' function sums the atom-style +variable over the reaction site, while the 'rxnave' returns the +average value. For example, a constraint on the total potential energy +of atoms involved in the reaction can be imposed as follows: .. code-block:: LAMMPS @@ -564,7 +564,7 @@ Arrhenius constraint that depends on the bond force of a specific bond: # in Constraints section of map file - custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3*rxnsum(v_ke)*v_natoms)) < random(0,1,12345)" + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3\*rxnsum(v_ke)\*v_natoms)) < random(0,1,12345)" By using an inequality and the 'random(x,y,z)' function, the left-hand side can be interpreted as the probability of the reaction occurring, From b19f40a855c3e5938bc76e117107f7833f8e9e39 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 26 Aug 2022 13:43:47 -0400 Subject: [PATCH 13/35] correct, simplify rxnbond example --- doc/src/fix_bond_react.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 3ed4e6ef26..c84644f39e 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -557,14 +557,14 @@ Arrhenius constraint that depends on the bond force of a specific bond: variable ke atom c_ke_atom variable E_a equal 100.0 # activation energy - variable natoms equal 10 # number of atoms in reaction site + variable l0 equal 1.0 # characteristic length .. code-block:: LAMMPS # in Constraints section of map file - custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3\*rxnsum(v_ke)\*v_natoms)) < random(0,1,12345)" + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag)*v_l0)/(2/3*rxnave(v_ke))) < random(0,1,12345)" By using an inequality and the 'random(x,y,z)' function, the left-hand side can be interpreted as the probability of the reaction occurring, From 9e79a1ef76c9f6a183f5146f73e5921df58530dd Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 14:48:47 -0400 Subject: [PATCH 14/35] non-'bond/react ' changes --- src/compute.h | 5 +++++ src/compute_bond_local.cpp | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/src/compute.h b/src/compute.h index 3ffe60d18c..068f86ca9f 100644 --- a/src/compute.h +++ b/src/compute.h @@ -16,6 +16,9 @@ #include "pointers.h" // IWYU pragma: export +#include +#include + namespace LAMMPS_NS { class Compute : protected Pointers { @@ -94,6 +97,8 @@ class Compute : protected Pointers { double dof; // degrees-of-freedom for temperature + std::map, int> atoms2bond; // maps atom pair to index of local bond array + int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index f2603e8cdd..ef07456070 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -247,6 +247,7 @@ int ComputeBondLocal::compute_bonds(int flag) double mvv2e; double engpot,engtrans,engvib,engrot,fbond; double *ptr; + std::set aset; double **x = atom->x; double **v = atom->v; @@ -390,6 +391,11 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } + // populate map to access per-bond values from atom tags + + aset = {tag[atom1], tag[atom2]}; + atoms2bond.insert(std::make_pair(aset,m)); + if (nvalues == 1) ptr = &vlocal[m]; else ptr = alocal[m]; From d090fce26268eb5e8608f0c47cd707298e78897f Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 15:50:27 -0400 Subject: [PATCH 15/35] bond/react: per-bond custom constraint --- src/REACTION/fix_bond_react.cpp | 62 +++++++++++++++++++++++++++------ src/REACTION/fix_bond_react.h | 5 ++- 2 files changed, 56 insertions(+), 11 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 5dc5967c29..6cd91aad2e 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -138,10 +138,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : status = PROCEED; // reaction functions used by 'custom' constraint - nrxnfunction = 2; + nrxnfunction = 3; rxnfunclist.resize(nrxnfunction); + peratomflag.resize(nrxnfunction); rxnfunclist[0] = "rxnsum"; + peratomflag[0] = 1; rxnfunclist[1] = "rxnave"; + peratomflag[1] = 1; + rxnfunclist[2] = "rxnbond"; + peratomflag[2] = 0; nvvec = 0; ncustomvars = 0; vvec = nullptr; @@ -2118,7 +2123,7 @@ double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col) } /* ---------------------------------------------------------------------- -get per-atom variable names used by custom constraint +get per-atom variable names used by custom constraint ------------------------------------------------------------------------- */ void FixBondReact::customvarnames() @@ -2140,6 +2145,7 @@ void FixBondReact::customvarnames() // find next reaction special function occurrence pos1 = std::string::npos; for (int i = 0; i < nrxnfunction; i++) { + if (peratomflag[i] == 0) continue; pos = varstr.find(rxnfunclist[i],prev3+1); if (pos == std::string::npos) continue; if (pos < pos1) pos1 = pos; @@ -2261,12 +2267,55 @@ double FixBondReact::custom_constraint(const std::string& varstr) } /* ---------------------------------------------------------------------- -currently two 'rxn' functions: rxnsum and rxnave +currently three 'rxn' functions: rxnsum, rxnave, and rxnbond ------------------------------------------------------------------------- */ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& varid, const std::string& fragid) { + int ifrag = -1; + if (fragid != "all") { + ifrag = onemol->findfragment(fragid.c_str()); + if (ifrag < 0) error->one(FLERR,"Bond/react: Molecule fragment " + "in reaction special function does not exist"); + } + + // start with 'rxnbond' per-bond function + // for 'rxnbond', varid corresponds to 'compute bond/local' name, + // and fragid is a pre-reaction fragment containing the two atoms in the bond + if (rxnfunc == "rxnbond") { + int icompute,ibond,nsum; + double perbondval; + std::set aset; + std::string computeid = varid; + + if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " + "name should begin with 'c_'"); + computeid = computeid.substr(2); + icompute = modify->find_compute(computeid); + if (icompute < 0) error->one(FLERR,"Bond/react: Reaction special function compute name does not exist"); + cperbond = modify->compute[icompute]; + std::string compute_style = cperbond->style; + if (compute_style != "bond/local") error->one(FLERR,"Bond/react: Compute used by reaction " + "special function 'rxnbond' must be of style 'bond/local'"); + if (cperbond->size_local_cols > 0) error->one(FLERR,"Bond/react: 'Compute bond/local' used by reaction " + "special function 'rxnbond' must compute one value"); + + nsum = 0; + for (int i = 0; i < onemol->natoms; i++) { + if (onemol->fragmentmask[ifrag][i]) { + aset.insert(glove[i][1]); + nsum++; + } + } + if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " + "must contain exactly two atoms"); + + ibond = cperbond->atoms2bond.find(aset)->second; + perbondval = cperbond->vector_local[ibond]; + return perbondval; + } + int ivar = -1; for (int i = 0; i < ncustomvars; i++) { if (varid == customvarstrs[i]) { @@ -2280,13 +2329,6 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& error->one(FLERR,"Fix bond/react: Reaction special function variable " "name does not exist"); - int ifrag = -1; - if (fragid != "all") { - ifrag = onemol->findfragment(fragid.c_str()); - if (ifrag < 0) error->one(FLERR,"Fix bond/react: Molecule fragment " - "in reaction special function does not exist"); - } - int iatom; int nsum = 0; double sumvvec = 0; diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index e9ca7ec362..5923ad485b 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -25,6 +25,7 @@ FixStyle(bond/react,FixBondReact); #define LMP_FIX_BOND_REACT_H #include "fix.h" +#include "compute.h" namespace LAMMPS_NS { @@ -74,7 +75,8 @@ class FixBondReact : public Fix { int *nconstraints; char **constraintstr; int nrxnfunction; - std::vector rxnfunclist; + std::vector rxnfunclist; // lists current special rxn function + std::vector peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond) int narrhenius; int **var_flag, **var_id; // for keyword values with variable inputs int status; @@ -222,6 +224,7 @@ class FixBondReact : public Fix { std::vector customvarstrs; int nvvec; double **vvec; // per-atom vector to store variable constraint atom-style variable values + Compute *cperbond; // pointer to 'compute bond/local' used by custom constraint ('rxnbond' function) std::vector> constraints; // DEBUG From 5e9b4d86788589fe5e1701951c7f628a7deee10b Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:02:13 -0400 Subject: [PATCH 16/35] actually evaluate bond/local compute value (even when not printed on that timestep) --- src/REACTION/fix_bond_react.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 6cd91aad2e..c2a03b0515 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2311,6 +2311,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); + cperbond->compute_local(); ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From fe6cf36101eb36d1c1e42b69a4d3f2669cf9fe05 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 22:51:52 -0400 Subject: [PATCH 17/35] prevent multiple compute evaluations on a timestep --- src/REACTION/fix_bond_react.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index c2a03b0515..d439a236b9 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2311,7 +2311,9 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (nsum != 2) error->one(FLERR,"Bond/react: Molecule fragment of reaction special function 'rxnbond' " "must contain exactly two atoms"); - cperbond->compute_local(); + if (cperbond->invoked_local != lmp->update->ntimestep) + cperbond->compute_local(); + ibond = cperbond->atoms2bond.find(aset)->second; perbondval = cperbond->vector_local[ibond]; return perbondval; From 23353208f28e80e24d6106ec31724adefe212ad3 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 29 May 2022 10:43:36 -0400 Subject: [PATCH 18/35] check that bond is actually in map --- src/REACTION/fix_bond_react.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index d439a236b9..8970ad9b0c 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -2288,6 +2288,7 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& double perbondval; std::set aset; std::string computeid = varid; + std::map,int>::iterator it; if (computeid.substr(0,2) != "c_") error->one(FLERR,"Bond/react: Reaction special function compute " "name should begin with 'c_'"); @@ -2314,7 +2315,10 @@ double FixBondReact::rxnfunction(const std::string& rxnfunc, const std::string& if (cperbond->invoked_local != lmp->update->ntimestep) cperbond->compute_local(); - ibond = cperbond->atoms2bond.find(aset)->second; + it = cperbond->atoms2bond.find(aset); + if (it == cperbond->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]; return perbondval; } From f0a421eb2572bff2dc07540c68a6fc37ea304437 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 29 Jul 2022 17:24:07 -0400 Subject: [PATCH 19/35] bond:react per-bond custom constraint docs --- doc/src/fix_bond_react.rst | 69 +++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 19 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index c557d2617d..2066b6e1c4 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -520,22 +520,32 @@ The constraint of type "custom" has the following syntax: custom *varstring* -where "custom" is the required keyword, and *varstring* is a -variable expression. The expression must be a valid equal-style -variable formula that can be read by the :doc:`variable ` command, +where 'custom' is the required keyword, and *varstring* is a variable +expression. The expression must be a valid equal-style variable +formula that can be read by the :doc:`variable ` command, after any special reaction functions are evaluated. If the resulting expression is zero, the reaction is prevented from occurring; -otherwise, it is permitted to occur. There are two special reaction -functions available, "rxnsum" and "rxnave". These functions operate -over the atoms in a given reaction site, and have one mandatory -argument and one optional argument. The mandatory argument is the -identifier for an atom-style variable. The second, optional argument -is the name of a molecule fragment in the pre-reaction template, and -can be used to operate over a subset of atoms in the reaction site. -The "rxnsum" function sums the atom-style variable over the reaction -site, while the "rxnave" returns the average value. For example, a -constraint on the total potential energy of atoms involved in the -reaction can be imposed as follows: +otherwise, it is permitted to occur. There are three special reaction +functions available, 'rxnbond', 'rxnsum', and 'rxnave'. The 'rxnbond' +function allows per-bond values to be included in the variable strings +of the custom constraint. The 'rxnbond' function has two mandatory +arguments. The first argument is the ID of a previously defined +'compute bond/local' command. This 'compute bond/local' must compute +only one value, e.g. bond force. This value is returned by the +'rxnbond' function. The second argument is the name of a molecule +fragment in the pre-reaction template. The fragment must contain +exactly two atoms, corresponding to the atoms involved in the bond +whose value should be calculated. An example of a constraint that uses +the force experienced by a bond in the reaction site is provided +below. The 'rxnsum' and 'rxnave' functions operate over the atoms in a +given reaction site, and have one mandatory argument and one optional +argument. The mandatory argument is the identifier for an atom-style +variable. The second, optional argument is the name of a molecule +fragment in the pre-reaction template, and can be used to operate over +a subset of atoms in the reaction site. The 'rxnsum' function sums the +atom-style variable over the reaction site, while the 'rxnave' returns +the average value. For example, a constraint on the total potential +energy of atoms involved in the reaction can be imposed as follows: .. code-block:: LAMMPS @@ -547,11 +557,32 @@ reaction can be imposed as follows: custom "rxnsum(v_my_pe) > 100" # in Constraints section of map file The above example prevents the reaction from occurring unless the -total potential energy of the reaction site is above 100. The variable -expression can be interpreted as the probability of the reaction -occurring by using an inequality and the :doc:`random(x,y,z) ` -function available for equal-style variables, similar to the 'arrhenius' -constraint above. +total potential energy of the reaction site is above 100. As a second +example, this time using the 'rxnbond' function, consider a modified +Arrhenius constraint that depends on the bond force of a specific bond: + +.. code-block:: LAMMPS + + # in LAMMPS input script + + compute bondforce all bond/local force + + compute ke_atom all ke/atom + variable ke atom c_ke_atom + + variable E_a equal 100.0 # activation energy + variable natoms equal 10 # number of atoms in reaction site + + +.. code-block:: LAMMPS + + # in Constraints section of map file + + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3*rxnsum(v_ke)*v_natoms)) < random(0,1,12345)" + +By using an inequality and the 'random(x,y,z)' function, the left-hand +side can be interpreted as the probability of the reaction occurring, +similar to the 'arrhenius' constraint above. By default, all constraints must be satisfied for the reaction to occur. In other words, constraints are evaluated as a series of From ed6fe96909cc565f751989f9a599bb07a1ade215 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 29 Jul 2022 23:34:37 -0400 Subject: [PATCH 20/35] simplify, formatting --- doc/src/fix_bond_react.rst | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 2066b6e1c4..3ad9e64252 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -536,16 +536,16 @@ only one value, e.g. bond force. This value is returned by the fragment in the pre-reaction template. The fragment must contain exactly two atoms, corresponding to the atoms involved in the bond whose value should be calculated. An example of a constraint that uses -the force experienced by a bond in the reaction site is provided -below. The 'rxnsum' and 'rxnave' functions operate over the atoms in a -given reaction site, and have one mandatory argument and one optional -argument. The mandatory argument is the identifier for an atom-style -variable. The second, optional argument is the name of a molecule -fragment in the pre-reaction template, and can be used to operate over -a subset of atoms in the reaction site. The 'rxnsum' function sums the -atom-style variable over the reaction site, while the 'rxnave' returns -the average value. For example, a constraint on the total potential -energy of atoms involved in the reaction can be imposed as follows: +the force experienced by a bond is provided below. The 'rxnsum' and +'rxnave' functions operate over the atoms in a given reaction site, +and have one mandatory argument and one optional argument. The +mandatory argument is the identifier for an atom-style variable. The +second, optional argument is the name of a molecule fragment in the +pre-reaction template, and can be used to operate over a subset of +atoms in the reaction site. The 'rxnsum' function sums the atom-style +variable over the reaction site, while the 'rxnave' returns the +average value. For example, a constraint on the total potential energy +of atoms involved in the reaction can be imposed as follows: .. code-block:: LAMMPS @@ -578,7 +578,7 @@ Arrhenius constraint that depends on the bond force of a specific bond: # in Constraints section of map file - custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3*rxnsum(v_ke)*v_natoms)) < random(0,1,12345)" + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3\*rxnsum(v_ke)\*v_natoms)) < random(0,1,12345)" By using an inequality and the 'random(x,y,z)' function, the left-hand side can be interpreted as the probability of the reaction occurring, From 189f803ef619bbf1275fc0031229b534de762ccf Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 26 Aug 2022 13:43:47 -0400 Subject: [PATCH 21/35] correct, simplify rxnbond example --- doc/src/fix_bond_react.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 3ad9e64252..b66b8a8bc2 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -571,14 +571,14 @@ Arrhenius constraint that depends on the bond force of a specific bond: variable ke atom c_ke_atom variable E_a equal 100.0 # activation energy - variable natoms equal 10 # number of atoms in reaction site + variable l0 equal 1.0 # characteristic length .. code-block:: LAMMPS # in Constraints section of map file - custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag))/(2/3\*rxnsum(v_ke)\*v_natoms)) < random(0,1,12345)" + custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag)*v_l0)/(2/3*rxnave(v_ke))) < random(0,1,12345)" By using an inequality and the 'random(x,y,z)' function, the left-hand side can be interpreted as the probability of the reaction occurring, From 736b420a490483c2d1c5acde16f283a91e4cf4a8 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 1 Nov 2022 18:03:30 -0400 Subject: [PATCH 22/35] reaction rate limit option --- src/REACTION/fix_bond_react.cpp | 153 +++++++++++++++++++++----------- src/REACTION/fix_bond_react.h | 9 +- 2 files changed, 107 insertions(+), 55 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 8970ad9b0c..859abe135a 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -80,7 +80,7 @@ static const char cite_fix_bond_react[] = #define DELTA 16 #define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm #define MAXCONARGS 14 // max # of arguments for any type of constraint + rxnID -#define NUMVARVALS 4 // max # of keyword values that have variables as input +#define NUMVARVALS 5 // max # of keyword values that have variables as input // various statuses of superimpose algorithm: // ACCEPT: site successfully matched to pre-reacted template @@ -98,7 +98,7 @@ enum{DISTANCE,ANGLE,DIHEDRAL,ARRHENIUS,RMSD,CUSTOM}; enum{ATOM,FRAG}; // keyword values that accept variables as input -enum{NEVERY,RMIN,RMAX,PROB}; +enum{NEVERY,RMIN,RMAX,PROB,NRATE}; // flag for one-proc vs shared reaction sites enum{LOCAL,GLOBAL}; @@ -229,6 +229,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips"); memory->create(seed,nreacts,"bond/react:seed"); memory->create(limit_duration,nreacts,"bond/react:limit_duration"); + memory->create(rate_limit,3,nreacts,"bond/react:rate_limit"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid"); memory->create(create_atoms_flag,nreacts,"bond/react:create_atoms_flag"); @@ -254,6 +255,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fraction[i] = 1; seed[i] = 12345; max_rxn[i] = INT_MAX; + for (int j = 0; j < 3; j++) + rate_limit[j][i] = 0; stabilize_steps_flag[i] = 0; custom_charges_fragid[i] = -1; create_atoms_flag[i] = 0; @@ -291,55 +294,31 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (groupid == -1) error->all(FLERR,"Could not find fix group ID"); groupbits[rxn] = group->bitmask[groupid]; - if (strncmp(arg[iarg],"v_",2) == 0) { - const char *str = &arg[iarg][2]; - var_id[NEVERY][rxn] = input->variable->find(str); - if (var_id[NEVERY][rxn] < 0) - error->all(FLERR,"Fix bond/react: Variable name {} does not exist", str); - if (!input->variable->equalstyle(var_id[NEVERY][rxn])) - error->all(FLERR,"Fix bond/react: Variable {} is not equal-style", str); - var_flag[NEVERY][rxn] = 1; - } else { + if (strncmp(arg[iarg],"v_",2) == 0) read_variable_keyword(&arg[iarg][2],NEVERY,rxn); + else { nevery[rxn] = utils::inumeric(FLERR,arg[iarg],false,lmp); if (nevery[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command: " "'Nevery' must be a positive integer"); } iarg++; + double cutoff; if (strncmp(arg[iarg],"v_",2) == 0) { - const char *str = &arg[iarg][2]; - var_id[RMIN][rxn] = input->variable->find(str); - if (var_id[RMIN][rxn] < 0) - error->all(FLERR,"Fix bond/react: Variable name {} does not exist", str); - if (!input->variable->equalstyle(var_id[RMIN][rxn])) - error->all(FLERR,"Fix bond/react: Variable {} is not equal-style", str); - double cutoff = input->variable->compute_equal(var_id[RMIN][rxn]); - cutsq[rxn][0] = cutoff*cutoff; - var_flag[RMIN][rxn] = 1; - } else { - double cutoff = utils::numeric(FLERR,arg[iarg],false,lmp); + read_variable_keyword(&arg[iarg][2],RMIN,rxn); + cutoff = input->variable->compute_equal(var_id[RMIN][rxn]); + } else cutoff = utils::numeric(FLERR,arg[iarg],false,lmp); if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/react command: " "'Rmin' cannot be negative"); cutsq[rxn][0] = cutoff*cutoff; - } iarg++; if (strncmp(arg[iarg],"v_",2) == 0) { - const char *str = &arg[iarg][2]; - var_id[RMAX][rxn] = input->variable->find(str); - if (var_id[RMAX][rxn] < 0) - error->all(FLERR,"Fix bond/react: Variable name {} does not exist", str); - if (!input->variable->equalstyle(var_id[RMAX][rxn])) - error->all(FLERR,"Fix bond/react: Variable is {} not equal-style", str); - double cutoff = input->variable->compute_equal(var_id[RMAX][rxn]); - cutsq[rxn][1] = cutoff*cutoff; - var_flag[RMAX][rxn] = 1; - } else { - double cutoff = utils::numeric(FLERR,arg[iarg],false,lmp); + read_variable_keyword(&arg[iarg][2],RMAX,rxn); + cutoff = input->variable->compute_equal(var_id[RMAX][rxn]); + } else cutoff = utils::numeric(FLERR,arg[iarg],false,lmp); if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/react command:" "'Rmax' cannot be negative"); cutsq[rxn][1] = cutoff*cutoff; - } iarg++; unreacted_mol[rxn] = atom->find_molecule(arg[iarg++]); @@ -349,7 +328,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (reacted_mol[rxn] == -1) error->all(FLERR,"Reacted molecule template ID for " "fix bond/react does not exist"); - // read superimpose file + //read map file files[rxn] = utils::strdup(arg[iarg]); iarg++; @@ -359,14 +338,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : "'prob' keyword has too few arguments"); // check if probability is a variable if (strncmp(arg[iarg+1],"v_",2) == 0) { - const char *str = &arg[iarg+1][2]; - var_id[PROB][rxn] = input->variable->find(str); - if (var_id[PROB][rxn] < 0) - error->all(FLERR,"Fix bond/react: Variable name {} does not exist", str); - if (!input->variable->equalstyle(var_id[PROB][rxn])) - error->all(FLERR,"Fix bond/react: Variable {} is not equal-style", str); + read_variable_keyword(&arg[iarg+1][2],PROB,rxn); fraction[rxn] = input->variable->compute_equal(var_id[PROB][rxn]); - var_flag[PROB][rxn] = 1; } else { // otherwise probability should be a number fraction[rxn] = utils::numeric(FLERR,arg[iarg+1],false,lmp); @@ -385,6 +358,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (max_rxn[rxn] < 0) error->all(FLERR,"Illegal fix bond/react command: " "'max_rxn' cannot be negative"); iarg += 2; + } else if (strcmp(arg[iarg],"rate_limit") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'rate_limit' has too few arguments"); + rate_limit[0][rxn] = 1; // serves as flag for rate_limit keyword + if (strncmp(arg[iarg+1],"v_",2) == 0) read_variable_keyword(&arg[iarg+1][2],NRATE,rxn); + else rate_limit[1][rxn] = utils::numeric(FLERR,arg[iarg+1],false,lmp); + rate_limit[2][rxn] = utils::numeric(FLERR,arg[iarg+2],false,lmp); + iarg += 3; } else if (strcmp(arg[iarg],"stabilize_steps") == 0) { if (stabilization_flag == 0) error->all(FLERR,"Stabilize_steps keyword " "used without stabilization keyword"); @@ -438,21 +419,24 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } max_natoms = 0; // the number of atoms in largest molecule template + int max_rate_limit_steps = 0; for (int myrxn = 0; myrxn < nreacts; myrxn++) { twomol = atom->molecules[reacted_mol[myrxn]]; max_natoms = MAX(max_natoms,twomol->natoms); + max_rate_limit_steps = MAX(max_rate_limit_steps,rate_limit[2][myrxn]); } memory->create(equivalences,max_natoms,2,nreacts,"bond/react:equivalences"); memory->create(reverse_equiv,max_natoms,2,nreacts,"bond/react:reverse_equiv"); memory->create(edge,max_natoms,nreacts,"bond/react:edge"); memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms"); + memory->create(store_rxn_count,max_rate_limit_steps,nreacts,"bond/react:store_rxn_count"); memory->create(custom_charges,max_natoms,nreacts,"bond/react:custom_charges"); memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms"); memory->create(create_atoms,max_natoms,nreacts,"bond/react:create_atoms"); memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms"); - for (int j = 0; j < nreacts; j++) + for (int j = 0; j < nreacts; j++) { for (int i = 0; i < max_natoms; i++) { edge[i][j] = 0; custom_charges[i][j] = 1; // update all partial charges by default @@ -467,6 +451,10 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : equivalences[i][m][j] = i+1; } } + for (int i = 0; i < max_rate_limit_steps; i++) { + store_rxn_count[i][j] = -1; + } + } // read all map files afterward for (int i = 0; i < nreacts; i++) { @@ -476,7 +464,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : onemol->check_attributes(); twomol->check_attributes(); get_molxspecials(); - read(i); + read_map_file(i); fclose(fp); if (ncreate == 0 && onemol->natoms != twomol->natoms) error->all(FLERR,"Fix bond/react: Reaction templates must contain the same number of atoms"); @@ -603,6 +591,7 @@ FixBondReact::~FixBondReact() memory->destroy(equivalences); memory->destroy(reverse_equiv); memory->destroy(landlocked_atoms); + memory->destroy(store_rxn_count); memory->destroy(custom_charges); memory->destroy(delete_atoms); memory->destroy(create_atoms); @@ -622,6 +611,7 @@ FixBondReact::~FixBondReact() memory->destroy(limit_duration); memory->destroy(var_flag); memory->destroy(var_id); + memory->destroy(rate_limit); memory->destroy(stabilize_steps_flag); memory->destroy(custom_charges_fragid); memory->destroy(molecule_keyword); @@ -824,6 +814,15 @@ void FixBondReact::init_list(int /*id*/, NeighList *ptr) void FixBondReact::post_integrate() { + // update store_rxn_count on every step + for (int myrxn = 0; myrxn < nreacts; myrxn++) { + if (rate_limit[0][myrxn] == 1) { + for (int i = rate_limit[2][myrxn]-1; i > 0; i--) + store_rxn_count[i][myrxn] = store_rxn_count[i-1][myrxn]; + store_rxn_count[0][myrxn] = reaction_count_total[myrxn]; + } + } + // check if any reactions could occur on this timestep int nevery_check = 1; for (int i = 0; i < nreacts; i++) { @@ -916,8 +915,22 @@ void FixBondReact::post_integrate() int j; for (rxnID = 0; rxnID < nreacts; rxnID++) { + int rate_limit_flag = 1; + if (rate_limit[0][rxnID] == 1) { + int myrxn_count = store_rxn_count[rate_limit[2][rxnID]-1][rxnID]; + if (myrxn_count == -1) rate_limit_flag = 0; + else { + int nrxns_delta = reaction_count_total[rxnID] - myrxn_count; + int my_nrate; + if (var_flag[NRATE][rxnID] == 1) { + my_nrate = input->variable->compute_equal(var_id[NRATE][rxnID]); + } else my_nrate = rate_limit[1][rxnID]; + if (nrxns_delta > my_nrate) rate_limit_flag = 0; + } + } if ((update->ntimestep % nevery[rxnID]) || - (max_rxn[rxnID] <= reaction_count_total[rxnID])) continue; + (max_rxn[rxnID] <= reaction_count_total[rxnID]) || + (rate_limit_flag == 0)) continue; for (int ii = 0; ii < nall; ii++) { partner[ii] = 0; finalpartner[ii] = 0; @@ -1390,9 +1403,25 @@ void FixBondReact::superimpose_algorithm() std::random_device rnd; std::minstd_rand park_rng(rnd()); - // check if we overstepped our reaction limit + // check if we overstepped our reaction limit, via either max_rxn or rate_limit for (int i = 0; i < nreacts; i++) { - if (reaction_count_total[i] > max_rxn[i]) { + int overstep = 0; + int max_rxn_overstep = reaction_count_total[i] - max_rxn[i]; + overstep = MAX(overstep,max_rxn_overstep); + if (rate_limit[0][i] == 1) { + int myrxn_count = store_rxn_count[rate_limit[2][i]-1][i]; + if (myrxn_count != -1) { + int nrxn_delta = reaction_count_total[i] - myrxn_count; + int my_nrate; + if (var_flag[NRATE][i] == 1) { + my_nrate = input->variable->compute_equal(var_id[NRATE][i]); + } else my_nrate = rate_limit[1][i]; + int rate_limit_overstep = nrxn_delta - my_nrate; + overstep = MAX(overstep,rate_limit_overstep); + } + } + + if (overstep > 0) { // let's randomly choose rxns to skip, unbiasedly from local and ghostly int *local_rxncounts; int *all_localskips; @@ -1400,8 +1429,11 @@ void FixBondReact::superimpose_algorithm() memory->create(all_localskips,nprocs,"bond/react:all_localskips"); MPI_Gather(&local_rxn_count[i],1,MPI_INT,local_rxncounts,1,MPI_INT,0,world); if (me == 0) { - int overstep = reaction_count_total[i] - max_rxn[i]; int delta_rxn = reaction_count[i] + ghostly_rxn_count[i]; + // when using variable input for rate_limit, rate_limit_overstep could be > delta_rxn (below) + // we need to limit overstep to the number of reactions on this timestep + // essentially skipping all reactions, would be more efficient to use a skip_all flag + if (overstep > delta_rxn) overstep = delta_rxn; int *rxn_by_proc; memory->create(rxn_by_proc,delta_rxn,"bond/react:rxn_by_proc"); for (int j = 0; j < delta_rxn; j++) @@ -1419,14 +1451,15 @@ void FixBondReact::superimpose_algorithm() else all_localskips[rxn_by_proc[j]]++; } memory->destroy(rxn_by_proc); + reaction_count_total[i] -= overstep; } - reaction_count_total[i] = max_rxn[i]; MPI_Scatter(&all_localskips[0],1,MPI_INT,&nlocalskips[i],1,MPI_INT,0,world); MPI_Bcast(&nghostlyskips[i],1,MPI_INT,0,world); memory->destroy(local_rxncounts); memory->destroy(all_localskips); } } + MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world); // this updates topology next step next_reneighbor = update->ntimestep; @@ -3034,6 +3067,8 @@ void FixBondReact::update_everything() } delete [] iskip; + if (update_num_mega == 0) continue; // can i do this? + // if inserted atoms and global map exists, reset map now instead // of waiting for comm since other pre-exchange fixes may use it // invoke map_init() b/c atom count has grown @@ -3844,10 +3879,24 @@ int FixBondReact::insert_atoms(tagint **my_mega_glove, int iupdate) } /* ---------------------------------------------------------------------- -read superimpose file +add equal-style variable to keyword argument list ------------------------------------------------------------------------- */ -void FixBondReact::read(int myrxn) +void FixBondReact::read_variable_keyword(const char *myarg, int keyword, int myrxn) +{ + var_id[keyword][myrxn] = input->variable->find(myarg); + if (var_id[keyword][myrxn] < 0) + error->all(FLERR,"Fix bond/react: Variable name {} does not exist",myarg); + if (!input->variable->equalstyle(var_id[keyword][myrxn])) + error->all(FLERR,"Fix bond/react: Variable {} is not equal-style",myarg); + var_flag[keyword][myrxn] = 1; +} + +/* ---------------------------------------------------------------------- +read map file +------------------------------------------------------------------------- */ + +void FixBondReact::read_map_file(int myrxn) { char line[MAXLINE],keyword[MAXLINE]; char *eof,*ptr; diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index 5923ad485b..5a6e62f81a 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -65,6 +65,8 @@ class FixBondReact : public Fix { int stabilization_flag; int reset_mol_ids_flag; int custom_exclude_flag; + int **rate_limit; + int **store_rxn_count; int *stabilize_steps_flag; int *custom_charges_fragid; int *create_atoms_flag; @@ -162,7 +164,8 @@ class FixBondReact : public Fix { // but whose first neighbors haven't int glove_counter; // used to determine when to terminate Superimpose Algorithm - void read(int); + void read_variable_keyword(const char *, int, int); + void read_map_file(int); void EdgeIDs(char *, int); void Equivalences(char *, int); void DeleteAtoms(char *, int); @@ -200,7 +203,7 @@ class FixBondReact : public Fix { void ghost_glovecast(); void update_everything(); int insert_atoms(tagint **, int); - void unlimit_bond(); + void unlimit_bond(); // removes atoms from stabilization, and other post-reaction every-step operations void limit_bond(int); void dedup_mega_gloves(int); //dedup global mega_glove void write_restart(FILE *) override; @@ -223,7 +226,7 @@ class FixBondReact : public Fix { int ncustomvars; std::vector customvarstrs; int nvvec; - double **vvec; // per-atom vector to store variable constraint atom-style variable values + 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::vector> constraints; From eb01f816ec4a2f68764d021ef37e021e483d7089 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Thu, 3 Nov 2022 11:09:44 -0400 Subject: [PATCH 23/35] rate_limit docs --- doc/src/fix_bond_react.rst | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index b66b8a8bc2..e6e137e1ea 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -42,13 +42,16 @@ Syntax * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates * zero or more individual keyword/value pairs may be appended to each react argument -* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create* +* individual_keyword = *prob* or *rate_limit* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create* .. parsed-literal:: *prob* values = fraction seed fraction = initiate reaction with this probability if otherwise eligible seed = random number seed (positive integer) + *rate_limit* = Nlimit Nsteps + Nlimit = maximum number of reactions allowed to occur within interval + Nsteps = the interval (number of timesteps) over which to count reactions *max_rxn* value = N N = maximum number of reactions allowed to occur *stabilize_steps* value = timesteps @@ -629,6 +632,12 @@ eligible reaction only occurs if the random number is less than the fraction. Up to :math:`N` reactions are permitted to occur, as optionally specified by the *max_rxn* keyword. +The *rate_limit* keyword can enforce an upper limit on the overall +rate of the reaction. The number of reaction occurences is limited to +Nlimit within an interval of Nsteps timesteps. At the beginning of a +run, including when reading from a restart file, no reactions are +permitted within the first Nsteps timesteps. + The *stabilize_steps* keyword allows for the specification of how many time steps a reaction site is stabilized before being returned to the overall system thermostat. In order to produce the most physical From 1a271b48703aee43b4a2ad68d4d6ac17b10652e0 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Thu, 3 Nov 2022 11:21:16 -0400 Subject: [PATCH 24/35] add mention of Nlimit variable support --- doc/src/fix_bond_react.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index e6e137e1ea..5c31dddb14 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -636,7 +636,8 @@ The *rate_limit* keyword can enforce an upper limit on the overall rate of the reaction. The number of reaction occurences is limited to Nlimit within an interval of Nsteps timesteps. At the beginning of a run, including when reading from a restart file, no reactions are -permitted within the first Nsteps timesteps. +permitted within the first Nsteps timesteps. Nlimit can be specified +with an equal-style :doc:`variable `. The *stabilize_steps* keyword allows for the specification of how many time steps a reaction site is stabilized before being returned to the From 0b7c391dfd069cd2b8ac3f5e5a4674568be36fb8 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 4 Nov 2022 00:49:26 -0400 Subject: [PATCH 25/35] smooth restart when using rate_limit breaks backward compatibility of restart files (probably can be avoided) --- src/REACTION/fix_bond_react.cpp | 31 ++++++++++++++++++++++++++----- src/REACTION/fix_bond_react.h | 2 ++ 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 859abe135a..1c91dd6ebd 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -419,7 +419,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } max_natoms = 0; // the number of atoms in largest molecule template - int max_rate_limit_steps = 0; + max_rate_limit_steps = 0; for (int myrxn = 0; myrxn < nreacts; myrxn++) { twomol = atom->molecules[reacted_mol[myrxn]]; max_natoms = MAX(max_natoms,twomol->natoms); @@ -817,8 +817,9 @@ void FixBondReact::post_integrate() // update store_rxn_count on every step for (int myrxn = 0; myrxn < nreacts; myrxn++) { if (rate_limit[0][myrxn] == 1) { - for (int i = rate_limit[2][myrxn]-1; i > 0; i--) + for (int i = rate_limit[2][myrxn]-1; i > 0; i--) { store_rxn_count[i][myrxn] = store_rxn_count[i-1][myrxn]; + } store_rxn_count[0][myrxn] = reaction_count_total[myrxn]; } } @@ -3067,7 +3068,7 @@ void FixBondReact::update_everything() } delete [] iskip; - if (update_num_mega == 0) continue; // can i do this? + if (update_num_mega == 0) continue; // if inserted atoms and global map exists, reset map now instead // of waiting for comm since other pre-exchange fixes may use it @@ -4398,17 +4399,25 @@ void FixBondReact::unpack_reverse_comm(int n, int *list, double *buf) void FixBondReact::write_restart(FILE *fp) { set[0].nreacts = nreacts; + set[0].max_rate_limit_steps = max_rate_limit_steps; for (int i = 0; i < nreacts; i++) { set[i].reaction_count_total = reaction_count_total[i]; strncpy(set[i].rxn_name,rxn_name[i],MAXLINE-1); set[i].rxn_name[MAXLINE-1] = '\0'; } + int *rbuf; + int rbufcount = max_rate_limit_steps*nreacts; + memory->create(rbuf,rbufcount,"bond/react:rbuf"); + memcpy(rbuf,&store_rxn_count[0][0],sizeof(int)*rbufcount); + if (me == 0) { - int size = nreacts*sizeof(Set); + int size = nreacts*sizeof(Set)+rbufcount*sizeof(int); fwrite(&size,sizeof(int),1,fp); fwrite(set,sizeof(Set),nreacts,fp); + fwrite(rbuf,sizeof(int),rbufcount,fp); } + memory->destroy(rbuf); } /* ---------------------------------------------------------------------- @@ -4418,13 +4427,25 @@ void FixBondReact::write_restart(FILE *fp) void FixBondReact::restart(char *buf) { Set *set_restart = (Set *) buf; - for (int i = 0; i < set_restart[0].nreacts; i++) { + int r_nreacts = set_restart[0].nreacts; + int r_max_rate_limit_steps = set_restart[0].max_rate_limit_steps; + int **ibuf; + int ibufcount = r_max_rate_limit_steps*r_nreacts; + memory->create(ibuf,r_max_rate_limit_steps,r_nreacts,"bond/react:ibuf"); + memcpy(&ibuf[0][0],&buf[nreacts*sizeof(Set)],sizeof(int)*ibufcount); + int n2cpy = r_max_rate_limit_steps; + if (max_rate_limit_steps < n2cpy) n2cpy = max_rate_limit_steps; + for (int i = 0; i < r_nreacts; i++) { for (int j = 0; j < nreacts; j++) { if (strcmp(set_restart[i].rxn_name,rxn_name[j]) == 0) { reaction_count_total[j] = set_restart[i].reaction_count_total; + // read rate_limit restart information + for (int k = 0; k < n2cpy; k++) + store_rxn_count[k][j] = ibuf[k][i]; } } } + memory->destroy(ibuf); } /* ---------------------------------------------------------------------- diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index 5a6e62f81a..452338ece4 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -90,6 +90,7 @@ class FixBondReact : public Fix { int *reaction_count_total; int nmax; // max num local atoms int max_natoms; // max natoms in a molecule template + int max_rate_limit_steps; // max rate limit interval tagint *partner, *finalpartner; double **distsq; int *nattempt; @@ -213,6 +214,7 @@ class FixBondReact : public Fix { int nreacts; char rxn_name[MAXLINE]; int reaction_count_total; + int max_rate_limit_steps; }; Set *set; From d5929a5cf80eff406b5ed66f43998671640322fb Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 4 Nov 2022 23:52:13 -0400 Subject: [PATCH 26/35] reflect rate_limit restart support in docs --- doc/src/fix_bond_react.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 5c31dddb14..36de92e50a 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -634,10 +634,10 @@ specified by the *max_rxn* keyword. The *rate_limit* keyword can enforce an upper limit on the overall rate of the reaction. The number of reaction occurences is limited to -Nlimit within an interval of Nsteps timesteps. At the beginning of a -run, including when reading from a restart file, no reactions are -permitted within the first Nsteps timesteps. Nlimit can be specified -with an equal-style :doc:`variable `. +Nlimit within an interval of Nsteps timesteps. No reactions are +permitted to occur within the first Nsteps timesteps of the first run +after reading a data file. Nlimit can be specified with an equal-style +:doc:`variable `. The *stabilize_steps* keyword allows for the specification of how many time steps a reaction site is stabilized before being returned to the From 565853ee070465ddae3e9375bc0459adae8f5e93 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 5 Nov 2022 01:05:45 -0400 Subject: [PATCH 27/35] rescale_charges keyword --- src/REACTION/fix_bond_react.cpp | 39 ++++++++++++++++++++++++++++++++- src/REACTION/fix_bond_react.h | 1 + 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 1c91dd6ebd..07de44d7f7 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -232,6 +232,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(rate_limit,3,nreacts,"bond/react:rate_limit"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid"); + memory->create(rescale_charges_flag,nreacts,"bond/react:rescale_charges_flag"); memory->create(create_atoms_flag,nreacts,"bond/react:create_atoms_flag"); memory->create(modify_create_fragid,nreacts,"bond/react:modify_create_fragid"); memory->create(overlapsq,nreacts,"bond/react:overlapsq"); @@ -259,6 +260,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : rate_limit[j][i] = 0; stabilize_steps_flag[i] = 0; custom_charges_fragid[i] = -1; + rescale_charges_flag[i] = 0; create_atoms_flag[i] = 0; modify_create_fragid[i] = -1; overlapsq[i] = 0; @@ -384,6 +386,13 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : "'custom_charges' keyword does not exist"); } iarg += 2; + } else if (strcmp(arg[iarg],"rescale_charges") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'rescale_charges' has too few arguments"); + if (strcmp(arg[iarg+1],"no") == 0) rescale_charges_flag[rxn] = 0; //default + else if (strcmp(arg[iarg+1],"yes") == 0) rescale_charges_flag[rxn] = 1; + else error->one(FLERR,"Bond/react: Illegal option for 'rescale_charges' keyword"); + iarg += 2; } else if (strcmp(arg[iarg],"molecule") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " "'molecule' has too few arguments"); @@ -614,6 +623,7 @@ FixBondReact::~FixBondReact() memory->destroy(rate_limit); memory->destroy(stabilize_steps_flag); memory->destroy(custom_charges_fragid); + memory->destroy(rescale_charges_flag); memory->destroy(molecule_keyword); memory->destroy(nconstraints); memory->destroy(constraintstr); @@ -3097,6 +3107,33 @@ void FixBondReact::update_everything() } } + // update charges and types of landlocked atoms + double charge_rescale_addend = 0; + if (rescale_charges_flag[rxnID] == 1) { + double sim_total_charge = 0; + double mol_total_charge = 0; + int n_custom_charge = 0; + 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 (atom->map(update_mega_glove[jj+1][i]) >= 0 && + atom->map(update_mega_glove[jj+1][i]) < nlocal) { + if (landlocked_atoms[j][rxnID] == 1) + type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; + if (twomol->qflag && atom->q_flag && custom_charges[jj][rxnID] == 1) { + double *q = atom->q; + sim_total_charge += q[atom->map(update_mega_glove[jj+1][i])]; + mol_total_charge += twomol->q[j]; + n_custom_charge++; + } + } + } + } + charge_rescale_addend = (sim_total_charge-mol_total_charge)/n_custom_charge; + } + // update charges and types of landlocked atoms for (int i = 0; i < update_num_mega; i++) { rxnID = update_mega_glove[0][i]; @@ -3109,7 +3146,7 @@ void FixBondReact::update_everything() type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; if (twomol->qflag && atom->q_flag && custom_charges[jj][rxnID] == 1) { double *q = atom->q; - q[atom->map(update_mega_glove[jj+1][i])] = twomol->q[j]; + q[atom->map(update_mega_glove[jj+1][i])] = twomol->q[j]+charge_rescale_addend; } } } diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index 452338ece4..628a71649b 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -69,6 +69,7 @@ class FixBondReact : public Fix { int **store_rxn_count; int *stabilize_steps_flag; int *custom_charges_fragid; + int *rescale_charges_flag; int *create_atoms_flag; int *modify_create_fragid; double *overlapsq; From b2652a4542d215ca33f2704858290b22d264a87d Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 5 Nov 2022 01:22:54 -0400 Subject: [PATCH 28/35] fix writing restarts without rate_limit --- src/REACTION/fix_bond_react.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 07de44d7f7..6d19a8a7a9 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -3107,7 +3107,7 @@ void FixBondReact::update_everything() } } - // update charges and types of landlocked atoms + // get charge rescale delta double charge_rescale_addend = 0; if (rescale_charges_flag[rxnID] == 1) { double sim_total_charge = 0; @@ -4437,24 +4437,28 @@ void FixBondReact::write_restart(FILE *fp) { set[0].nreacts = nreacts; set[0].max_rate_limit_steps = max_rate_limit_steps; + + if (set[0].max_rate_limit_steps > 0) for (int i = 0; i < nreacts; i++) { set[i].reaction_count_total = reaction_count_total[i]; strncpy(set[i].rxn_name,rxn_name[i],MAXLINE-1); set[i].rxn_name[MAXLINE-1] = '\0'; } - int *rbuf; int rbufcount = max_rate_limit_steps*nreacts; - memory->create(rbuf,rbufcount,"bond/react:rbuf"); - memcpy(rbuf,&store_rxn_count[0][0],sizeof(int)*rbufcount); + int *rbuf; + if (rbufcount) { + memory->create(rbuf,rbufcount,"bond/react:rbuf"); + memcpy(rbuf,&store_rxn_count[0][0],sizeof(int)*rbufcount); + } if (me == 0) { int size = nreacts*sizeof(Set)+rbufcount*sizeof(int); fwrite(&size,sizeof(int),1,fp); fwrite(set,sizeof(Set),nreacts,fp); - fwrite(rbuf,sizeof(int),rbufcount,fp); + if (rbufcount) fwrite(rbuf,sizeof(int),rbufcount,fp); } - memory->destroy(rbuf); + if (rbufcount) memory->destroy(rbuf); } /* ---------------------------------------------------------------------- From 3b47afa69f1833822e0350bb321f847174126761 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 5 Nov 2022 02:03:02 -0400 Subject: [PATCH 29/35] Update fix_bond_react.rst --- doc/src/fix_bond_react.rst | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 36de92e50a..575dbce528 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -42,7 +42,7 @@ Syntax * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates * zero or more individual keyword/value pairs may be appended to each react argument -* individual_keyword = *prob* or *rate_limit* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create* +* individual_keyword = *prob* or *rate_limit* or *max_rxn* or *stabilize_steps* or *custom_charges* or *rescale_charges* or *molecule* or *modify_create* .. parsed-literal:: @@ -59,6 +59,9 @@ Syntax *custom_charges* value = *no* or fragment-ID *no* = update all atomic charges (default) fragment-ID = ID of molecule fragment whose charges are updated + *rescale_charges* value = *no* or *yes* + *no* = do not rescale atomic charges (default) + *yes* = rescale charges such that total charge does not change during reaction *molecule* value = *off* or *inter* or *intra* *off* = allow both inter- and intramolecular reactions (default) *inter* = search for reactions between molecules with different IDs @@ -657,6 +660,17 @@ charges are updated to those specified by the post-reaction template fragment defined in the pre-reaction molecule template. In this case, only the atomic charges of atoms in the molecule fragment are updated. +The *rescale_charges* keyword can be used to ensure the total charge +of the system does not change as reactions occur. When the argument is +set to *yes*\ , a fixed value is added to the charges of post-reaction +atoms such that their total charge equals that of the pre-reaction +site. If only a subset of atomic charges are updated via the +*custom_charges* keyword, this rescaling is applied to the subset. +This keyword could be useful for systems that contain different +molecules with the same reactive site, if the partial charges on the +reaction site vary from molecule to molecule, or when removing +reaction by-products. + The *molecule* keyword can be used to force the reaction to be intermolecular, intramolecular or either. When the value is set to *off*\ , molecule IDs are not considered when searching for reactions From f9d07d8932277edceea4eb97d879d9cc9d41c2c3 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 9 Nov 2022 17:40:57 -0500 Subject: [PATCH 30/35] revert non-bond-react changes --- src/compute.h | 5 ----- src/compute_bond_local.cpp | 6 ------ 2 files changed, 11 deletions(-) diff --git a/src/compute.h b/src/compute.h index 068f86ca9f..3ffe60d18c 100644 --- a/src/compute.h +++ b/src/compute.h @@ -16,9 +16,6 @@ #include "pointers.h" // IWYU pragma: export -#include -#include - namespace LAMMPS_NS { class Compute : protected Pointers { @@ -97,8 +94,6 @@ class Compute : protected Pointers { double dof; // degrees-of-freedom for temperature - std::map, int> atoms2bond; // maps atom pair to index of local bond array - int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index ef07456070..f2603e8cdd 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -247,7 +247,6 @@ int ComputeBondLocal::compute_bonds(int flag) double mvv2e; double engpot,engtrans,engvib,engrot,fbond; double *ptr; - std::set aset; double **x = atom->x; double **v = atom->v; @@ -391,11 +390,6 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } - // populate map to access per-bond values from atom tags - - aset = {tag[atom1], tag[atom2]}; - atoms2bond.insert(std::make_pair(aset,m)); - if (nvalues == 1) ptr = &vlocal[m]; else ptr = alocal[m]; From aaa8e9d21953eaede0b81cc8f5bcf82623964209 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 9 Nov 2022 17:51:21 -0500 Subject: [PATCH 31/35] 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 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 +#include + namespace LAMMPS_NS { class FixBondReact : public Fix { @@ -79,9 +82,10 @@ class FixBondReact : public Fix { char **constraintstr; int nrxnfunction; std::vector rxnfunclist; // lists current special rxn function - std::vector peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond) + std::vector 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, int> atoms2bond; // maps atom pair to index of local bond array std::vector> constraints; // DEBUG From 0890bc026e9f79650d3c05fccb426b0cb69b45d0 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 11 Nov 2022 15:28:36 -0500 Subject: [PATCH 32/35] ensure backward compatibility of restarts --- src/REACTION/fix_bond_react.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index a2a0c00301..6f3af8953d 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -4503,18 +4503,23 @@ void FixBondReact::write_restart(FILE *fp) /* ---------------------------------------------------------------------- use selected state info from restart file to restart the Fix + bond/react restart format was updated after LAMMPS version 3 Nov 2022 ------------------------------------------------------------------------- */ void FixBondReact::restart(char *buf) { - Set *set_restart = (Set *) buf; - int r_nreacts = set_restart[0].nreacts; - int r_max_rate_limit_steps = set_restart[0].max_rate_limit_steps; + int r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy; int **ibuf; - int ibufcount = r_max_rate_limit_steps*r_nreacts; - memory->create(ibuf,r_max_rate_limit_steps,r_nreacts,"bond/react:ibuf"); - memcpy(&ibuf[0][0],&buf[nreacts*sizeof(Set)],sizeof(int)*ibufcount); - int n2cpy = r_max_rate_limit_steps; + Set *set_restart = (Set *) buf; + + r_nreacts = set_restart[0].nreacts; + if (lmp->restart_ver > utils::date2num("3 Nov 2022")) { + r_max_rate_limit_steps = set_restart[0].max_rate_limit_steps; + ibufcount = r_max_rate_limit_steps*r_nreacts; + memory->create(ibuf,r_max_rate_limit_steps,r_nreacts,"bond/react:ibuf"); + memcpy(&ibuf[0][0],&buf[nreacts*sizeof(Set)],sizeof(int)*ibufcount); + n2cpy = r_max_rate_limit_steps; + } else n2cpy = 0; if (max_rate_limit_steps < n2cpy) n2cpy = max_rate_limit_steps; for (int i = 0; i < r_nreacts; i++) { for (int j = 0; j < nreacts; j++) { @@ -4526,7 +4531,7 @@ void FixBondReact::restart(char *buf) } } } - memory->destroy(ibuf); + if (lmp->restart_ver > utils::date2num("3 Nov 2022")) memory->destroy(ibuf); } /* ---------------------------------------------------------------------- From 86172eb75f66a9045bd6cb66c200b636e07c5e99 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 11 Nov 2022 17:16:37 -0500 Subject: [PATCH 33/35] bond/react restart revision numbers --- src/REACTION/fix_bond_react.cpp | 19 ++++++++++++++----- src/REACTION/fix_bond_react.h | 1 + 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 6f3af8953d..e6aa5a9cfb 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -4475,6 +4475,7 @@ void FixBondReact::unpack_reverse_comm(int n, int *list, double *buf) void FixBondReact::write_restart(FILE *fp) { + int revision = 1; set[0].nreacts = nreacts; set[0].max_rate_limit_steps = max_rate_limit_steps; @@ -4493,8 +4494,9 @@ void FixBondReact::write_restart(FILE *fp) } if (me == 0) { - int size = nreacts*sizeof(Set)+rbufcount*sizeof(int); + int size = nreacts*sizeof(Set)+(rbufcount+1)*sizeof(int); fwrite(&size,sizeof(int),1,fp); + fwrite(&revision,sizeof(int),1,fp); fwrite(set,sizeof(Set),nreacts,fp); if (rbufcount) fwrite(rbuf,sizeof(int),rbufcount,fp); } @@ -4503,23 +4505,30 @@ void FixBondReact::write_restart(FILE *fp) /* ---------------------------------------------------------------------- use selected state info from restart file to restart the Fix - bond/react restart format was updated after LAMMPS version 3 Nov 2022 + bond/react restart revisions numbers added after LAMMPS version 3 Nov 2022 ------------------------------------------------------------------------- */ void FixBondReact::restart(char *buf) { - int r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy; + int revision,r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy; int **ibuf; Set *set_restart = (Set *) buf; + if (lmp->restart_ver > utils::date2num("3 Nov 2022")) { + revision = buf[0]; + buf++; + } else revision = 0; + r_nreacts = set_restart[0].nreacts; - if (lmp->restart_ver > utils::date2num("3 Nov 2022")) { + + if (revision > 0) { r_max_rate_limit_steps = set_restart[0].max_rate_limit_steps; ibufcount = r_max_rate_limit_steps*r_nreacts; memory->create(ibuf,r_max_rate_limit_steps,r_nreacts,"bond/react:ibuf"); memcpy(&ibuf[0][0],&buf[nreacts*sizeof(Set)],sizeof(int)*ibufcount); n2cpy = r_max_rate_limit_steps; } else n2cpy = 0; + if (max_rate_limit_steps < n2cpy) n2cpy = max_rate_limit_steps; for (int i = 0; i < r_nreacts; i++) { for (int j = 0; j < nreacts; j++) { @@ -4531,7 +4540,7 @@ void FixBondReact::restart(char *buf) } } } - if (lmp->restart_ver > utils::date2num("3 Nov 2022")) memory->destroy(ibuf); + if (revision > 0) memory->destroy(ibuf); } /* ---------------------------------------------------------------------- diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index 41908e40d4..b55608b07e 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -216,6 +216,7 @@ class FixBondReact : public Fix { void write_restart(FILE *) override; void restart(char *buf) override; + // store restart data struct Set { int nreacts; char rxn_name[MAXLINE]; From 72752931fac264831266130ccfa01e17db17fdea Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 11 Nov 2022 17:25:49 -0500 Subject: [PATCH 34/35] version tags --- doc/src/fix_bond_react.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 575dbce528..15098761cc 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -520,6 +520,8 @@ example, the molecule fragment could consist of only the backbone atoms of a polymer chain. This constraint can be used to enforce a specific relative position and orientation between reacting molecules. +.. versionchanged:: TBD + The constraint of type "custom" has the following syntax: .. parsed-literal:: @@ -635,6 +637,8 @@ eligible reaction only occurs if the random number is less than the fraction. Up to :math:`N` reactions are permitted to occur, as optionally specified by the *max_rxn* keyword. +.. versionadded:: TBD + The *rate_limit* keyword can enforce an upper limit on the overall rate of the reaction. The number of reaction occurences is limited to Nlimit within an interval of Nsteps timesteps. No reactions are @@ -660,6 +664,8 @@ charges are updated to those specified by the post-reaction template fragment defined in the pre-reaction molecule template. In this case, only the atomic charges of atoms in the molecule fragment are updated. +.. versionadded:: TBD + The *rescale_charges* keyword can be used to ensure the total charge of the system does not change as reactions occur. When the argument is set to *yes*\ , a fixed value is added to the charges of post-reaction From eee6862f2529ff63edfce30dc0642811c3496e77 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 12 Nov 2022 02:13:39 -0500 Subject: [PATCH 35/35] copy error + other fixes --- src/REACTION/fix_bond_react.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index e6aa5a9cfb..8857590c4f 100644 --- a/src/REACTION/fix_bond_react.cpp +++ b/src/REACTION/fix_bond_react.cpp @@ -4479,9 +4479,9 @@ void FixBondReact::write_restart(FILE *fp) set[0].nreacts = nreacts; set[0].max_rate_limit_steps = max_rate_limit_steps; - if (set[0].max_rate_limit_steps > 0) for (int i = 0; i < nreacts; i++) { set[i].reaction_count_total = reaction_count_total[i]; + strncpy(set[i].rxn_name,rxn_name[i],MAXLINE-1); set[i].rxn_name[MAXLINE-1] = '\0'; } @@ -4510,22 +4510,21 @@ void FixBondReact::write_restart(FILE *fp) void FixBondReact::restart(char *buf) { - int revision,r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy; + int n,revision,r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy; int **ibuf; - Set *set_restart = (Set *) buf; - if (lmp->restart_ver > utils::date2num("3 Nov 2022")) { - revision = buf[0]; - buf++; - } else revision = 0; + n = 0; + if (lmp->restart_ver > utils::date2num("3 Nov 2022")) revision = buf[n++]; + else revision = 0; + Set *set_restart = (Set *) &buf[n*sizeof(int)]; r_nreacts = set_restart[0].nreacts; if (revision > 0) { r_max_rate_limit_steps = set_restart[0].max_rate_limit_steps; ibufcount = r_max_rate_limit_steps*r_nreacts; memory->create(ibuf,r_max_rate_limit_steps,r_nreacts,"bond/react:ibuf"); - memcpy(&ibuf[0][0],&buf[nreacts*sizeof(Set)],sizeof(int)*ibufcount); + memcpy(&ibuf[0][0],&buf[sizeof(int)+r_nreacts*sizeof(Set)],sizeof(int)*ibufcount); n2cpy = r_max_rate_limit_steps; } else n2cpy = 0;