From 1b1cb5568d1f475ce0316a6d0487d2cf83e875cf Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 28 May 2022 14:48:47 -0400 Subject: [PATCH 01/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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/58] 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 0de50f29f73805c48f62aba59176c9fc0e5861b7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 4 Nov 2022 21:56:06 -0400 Subject: [PATCH 26/58] add option to use the MS-MPI SDK to cross-compile Windows binaries --- cmake/Modules/MPI4WIN.cmake | 107 ++++++++++++++++++++++++------------ 1 file changed, 71 insertions(+), 36 deletions(-) diff --git a/cmake/Modules/MPI4WIN.cmake b/cmake/Modules/MPI4WIN.cmake index aa0c9e1833..02db6d4744 100644 --- a/cmake/Modules/MPI4WIN.cmake +++ b/cmake/Modules/MPI4WIN.cmake @@ -1,39 +1,74 @@ -# Download and configure custom MPICH files for Windows -message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") -set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") -set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") -set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") -set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") -mark_as_advanced(MPICH2_WIN64_DEVEL_URL) -mark_as_advanced(MPICH2_WIN32_DEVEL_URL) -mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) -mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) +# Download and configure MinGW compatible MPICH development files for Windows +option(USE_MSMPI "Use Microsoft's MS-MPI SDK instead of MPICH2-1.4.1" OFF) -include(ExternalProject) -if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN64_DEVEL_URL} - URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) +if(USE_MSMPI) + message(STATUS "Downloading and configuring MS-MPI 10.1 for Windows cross-compilation") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/msmpi-win64-devel.tar.gz" CACHE STRING "URL for MS-MPI (win64) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "86314daf1bffb809f1fcbefb8a547490" CACHE STRING "MD5 checksum of MS-MPI (win64) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) + + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmsmpi.a) + else() + message(FATAL_ERROR "Only x86 64-bit builds are supported with MS-MPI") + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmsmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmsmpi.a") else() - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN32_DEVEL_URL} - URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) + message(STATUS "Downloading and configuring MPICH2-1.4.1 for Windows cross-compilation") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN32_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) + mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) + + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + else() + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN32_DEVEL_URL} + URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") endif() - -ExternalProject_get_property(mpi4win_build SOURCE_DIR) -file(MAKE_DIRECTORY "${SOURCE_DIR}/include") -add_library(MPI::MPI_CXX UNKNOWN IMPORTED) -set_target_properties(MPI::MPI_CXX PROPERTIES - IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" - INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" - INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") -add_dependencies(MPI::MPI_CXX mpi4win_build) - -# set variables for status reporting at the end of CMake run -set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") -set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") -set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") From 710197cd88402f66225db23e663ed1eed5fe1272 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 4 Nov 2022 23:00:34 -0400 Subject: [PATCH 27/58] add MS-MPI support to CMake support for plugins --- cmake/Modules/LAMMPSInterfacePlugin.cmake | 103 ++++++++++++++-------- 1 file changed, 67 insertions(+), 36 deletions(-) diff --git a/cmake/Modules/LAMMPSInterfacePlugin.cmake b/cmake/Modules/LAMMPSInterfacePlugin.cmake index 95bb707e86..e20ae5f5cb 100644 --- a/cmake/Modules/LAMMPSInterfacePlugin.cmake +++ b/cmake/Modules/LAMMPSInterfacePlugin.cmake @@ -112,45 +112,76 @@ if(BUILD_MPI) set(MPI_CXX_SKIP_MPICXX TRUE) # We use a non-standard procedure to cross-compile with MPI on Windows if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND CMAKE_CROSSCOMPILING) - # Download and configure custom MPICH files for Windows - message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") - set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") - set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") - set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") - set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") - mark_as_advanced(MPICH2_WIN64_DEVEL_URL) - mark_as_advanced(MPICH2_WIN32_DEVEL_URL) - mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) - mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) + # Download and configure MinGW compatible MPICH development files for Windows + option(USE_MSMPI "Use Microsoft's MS-MPI SDK instead of MPICH2-1.4.1" OFF) + if(USE_MSMPI) + message(STATUS "Downloading and configuring MS-MPI 10.1 for Windows cross-compilation") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/msmpi-win64-devel.tar.gz" CACHE STRING "URL for MS-MPI (win64) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "86314daf1bffb809f1fcbefb8a547490" CACHE STRING "MD5 checksum of MS-MPI (win64) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) - include(ExternalProject) - if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN64_DEVEL_URL} - URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmsmpi.a) + else() + message(FATAL_ERROR "Only x86 64-bit builds are supported with MS-MPI") + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmsmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmsmpi.a") else() - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN32_DEVEL_URL} - URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) + # Download and configure custom MPICH files for Windows + message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) + + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + else() + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN32_DEVEL_URL} + URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") endif() - - ExternalProject_get_property(mpi4win_build SOURCE_DIR) - file(MAKE_DIRECTORY "${SOURCE_DIR}/include") - add_library(MPI::MPI_CXX UNKNOWN IMPORTED) - set_target_properties(MPI::MPI_CXX PROPERTIES - IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" - INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" - INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") - add_dependencies(MPI::MPI_CXX mpi4win_build) - - # set variables for status reporting at the end of CMake run - set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") - set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") - set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") else() find_package(MPI REQUIRED) option(LAMMPS_LONGLONG_TO_LONG "Workaround if your system or MPI version does not recognize 'long long' data types" OFF) From bd5ea1b896d856433dd72940ca9bd0f676bd5bd8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 4 Nov 2022 23:32:40 -0400 Subject: [PATCH 28/58] add MS-MPI support to ML-PACE plugin and demo plugins --- examples/PACKAGES/pace/plugin/CMakeLists.txt | 12 ++- examples/plugins/LAMMPSInterfaceCXX.cmake | 107 ++++++++++++------- 2 files changed, 80 insertions(+), 39 deletions(-) diff --git a/examples/PACKAGES/pace/plugin/CMakeLists.txt b/examples/PACKAGES/pace/plugin/CMakeLists.txt index 6ad9c791ba..25fa877ffc 100644 --- a/examples/PACKAGES/pace/plugin/CMakeLists.txt +++ b/examples/PACKAGES/pace/plugin/CMakeLists.txt @@ -38,9 +38,15 @@ elseif(CMAKE_SYSTEM_NAME STREQUAL "Windows") execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_SOURCE_DIR}/lammps.ico ${CMAKE_SOURCE_DIR}/lammps-text-logo-wide.bmp ${CMAKE_SOURCE_DIR}/paceplugin.nsis ${CMAKE_BINARY_DIR}) if(BUILD_MPI) - add_custom_target(package ${MAKENSIS_PATH} -V1 -DVERSION=${LAMMPS_VERSION}-MPI paceplugin.nsis - DEPENDS paceplugin - BYPRODUCTS LAMMPS-ML-PACE-plugin-${LAMMPS_VERSION}-MPI.exe) + if(USE_MSMPI) + add_custom_target(package ${MAKENSIS_PATH} -V1 -DVERSION=${LAMMPS_VERSION}-MSMPI paceplugin.nsis + DEPENDS paceplugin + BYPRODUCTS LAMMPS-ML-PACE-plugin-${LAMMPS_VERSION}-MSMPI.exe) + else() + add_custom_target(package ${MAKENSIS_PATH} -V1 -DVERSION=${LAMMPS_VERSION}-MPI paceplugin.nsis + DEPENDS paceplugin + BYPRODUCTS LAMMPS-ML-PACE-plugin-${LAMMPS_VERSION}-MPI.exe) + endif() else() add_custom_target(package ${MAKENSIS_PATH} -V1 -DVERSION=${LAMMPS_VERSION} paceplugin.nsis COMMAND ${CMAKE_COMMAND} -E echo ${PWD} diff --git a/examples/plugins/LAMMPSInterfaceCXX.cmake b/examples/plugins/LAMMPSInterfaceCXX.cmake index a3313215d6..7eef5bd6e4 100644 --- a/examples/plugins/LAMMPSInterfaceCXX.cmake +++ b/examples/plugins/LAMMPSInterfaceCXX.cmake @@ -45,45 +45,80 @@ if(BUILD_MPI) set(MPI_CXX_SKIP_MPICXX TRUE) # We use a non-standard procedure to cross-compile with MPI on Windows if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND CMAKE_CROSSCOMPILING) - # Download and configure custom MPICH files for Windows - message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") - set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") - set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") - set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") - set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") - mark_as_advanced(MPICH2_WIN64_DEVEL_URL) - mark_as_advanced(MPICH2_WIN32_DEVEL_URL) - mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) - mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) + # Download and configure MinGW compatible MPICH development files for Windows + option(USE_MSMPI "Use Microsoft's MS-MPI SDK instead of MPICH2-1.4.1" OFF) + if(USE_MSMPI) + message(STATUS "Downloading and configuring MS-MPI 10.1 for Windows cross-compilation") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/msmpi-win64-devel.tar.gz" CACHE STRING "URL for MS-MPI (win64) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "86314daf1bffb809f1fcbefb8a547490" CACHE STRING "MD5 checksum of MS-MPI (win64) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) - include(ExternalProject) - if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN64_DEVEL_URL} - URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmsmpi.a) + else() + message(FATAL_ERROR "Only x86 64-bit builds are supported with MS-MPI") + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmsmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmsmpi.a") else() - ExternalProject_Add(mpi4win_build - URL ${MPICH2_WIN32_DEVEL_URL} - URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - BUILD_BYPRODUCTS /lib/libmpi.a) + # Download and configure custom MPICH files for Windows + message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN32_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) + mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) + + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + else() + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN32_DEVEL_URL} + URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") endif() - - ExternalProject_get_property(mpi4win_build SOURCE_DIR) - file(MAKE_DIRECTORY "${SOURCE_DIR}/include") - add_library(MPI::MPI_CXX UNKNOWN IMPORTED) - set_target_properties(MPI::MPI_CXX PROPERTIES - IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" - INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" - INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") - add_dependencies(MPI::MPI_CXX mpi4win_build) - - # set variables for status reporting at the end of CMake run - set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") - set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") - set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") else() find_package(MPI REQUIRED) option(LAMMPS_LONGLONG_TO_LONG "Workaround if your system or MPI version does not recognize 'long long' data types" OFF) From d5929a5cf80eff406b5ed66f43998671640322fb Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 4 Nov 2022 23:52:13 -0400 Subject: [PATCH 29/58] 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 30/58] 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 31/58] 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 32/58] 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 89f896ea7357d0380c833ad24361d458ef8551ae Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Nov 2022 23:57:29 -0400 Subject: [PATCH 33/58] Include MS-MPI in Windows build and test through GitHub Actions --- .github/workflows/compile-msvc.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/compile-msvc.yml b/.github/workflows/compile-msvc.yml index 3ec5fc4cb0..ee4044b2cd 100644 --- a/.github/workflows/compile-msvc.yml +++ b/.github/workflows/compile-msvc.yml @@ -33,6 +33,8 @@ jobs: run: | python3 -m pip install numpy python3 -m pip install pyyaml + nuget install MSMPIsdk + nuget install MSMPIDIST cmake -C cmake/presets/windows.cmake \ -D PKG_PYTHON=on \ -S cmake -B build \ From 09e490db404c98c7aa487a563d216a156dc695ad Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Nov 2022 07:24:17 -0500 Subject: [PATCH 34/58] add support for building/using the ADIOS package without MPI This needs the ADIOS2 installation being configured accordingly. --- cmake/CMakeLists.txt | 9 +++++++++ src/ADIOS/dump_atom_adios.cpp | 12 ++++++++++++ src/ADIOS/dump_custom_adios.cpp | 12 ++++++++++++ src/ADIOS/reader_adios.cpp | 8 ++++++++ 4 files changed, 41 insertions(+) diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index efdbfa2840..24522b6480 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -316,6 +316,15 @@ if(PKG_ADIOS) # script that defines the MPI::MPI_C target enable_language(C) find_package(ADIOS2 REQUIRED) + if(BUILD_MPI) + if(NOT ADIOS2_HAVE_MPI) + message(FATAL_ERROR "ADIOS2 must be built with MPI support when LAMMPS has MPI enabled") + endif() + else() + if(ADIOS2_HAVE_MPI) + message(FATAL_ERROR "ADIOS2 must be built without MPI support when LAMMPS has MPI disabled") + endif() + endif() target_link_libraries(lammps PRIVATE adios2::adios2) endif() diff --git a/src/ADIOS/dump_atom_adios.cpp b/src/ADIOS/dump_atom_adios.cpp index a31f3f4f25..deee23ba41 100644 --- a/src/ADIOS/dump_atom_adios.cpp +++ b/src/ADIOS/dump_atom_adios.cpp @@ -62,7 +62,11 @@ DumpAtomADIOS::DumpAtomADIOS(LAMMPS *lmp, int narg, char **arg) : DumpAtom(lmp, internal = new DumpAtomADIOSInternal(); try { +#if defined(MPI_STUBS) + internal->ad = new adios2::ADIOS("adios2_config.xml", adios2::DebugON); +#else internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON); +#endif } catch (std::ios_base::failure &e) { error->all(FLERR, "ADIOS initialization failed with error: {}", e.what()); } @@ -84,11 +88,19 @@ void DumpAtomADIOS::openfile() if (multifile) { // if one file per timestep, replace '*' with current timestep auto filecurrent = utils::star_subst(filename, update->ntimestep, padflag); +#if defined(MPI_STUBS) + internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write); +#else internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write, world); +#endif if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filecurrent); } else { if (!singlefile_opened) { +#if defined(MPI_STUBS) + internal->fh = internal->io.Open(filename, adios2::Mode::Write); +#else internal->fh = internal->io.Open(filename, adios2::Mode::Write, world); +#endif if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filename); singlefile_opened = 1; } diff --git a/src/ADIOS/dump_custom_adios.cpp b/src/ADIOS/dump_custom_adios.cpp index 194141c15a..be827f507d 100644 --- a/src/ADIOS/dump_custom_adios.cpp +++ b/src/ADIOS/dump_custom_adios.cpp @@ -70,7 +70,11 @@ DumpCustomADIOS::DumpCustomADIOS(LAMMPS *lmp, int narg, char **arg) : DumpCustom internal = new DumpCustomADIOSInternal(); try { +#if defined(MPI_STUBS) + internal->ad = new adios2::ADIOS("adios2_config.xml", adios2::DebugON); +#else internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON); +#endif } catch (std::ios_base::failure &e) { error->all(FLERR, "ADIOS initialization failed with error: {}", e.what()); } @@ -96,11 +100,19 @@ void DumpCustomADIOS::openfile() if (multifile) { // if one file per timestep, replace '*' with current timestep auto filecurrent = utils::star_subst(filename, update->ntimestep, padflag); +#if defined(MPI_STUBS) + internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write); +#else internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write, world); +#endif if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filecurrent); } else { if (!singlefile_opened) { +#if defined(MPI_STUBS) + internal->fh = internal->io.Open(filename, adios2::Mode::Write); +#else internal->fh = internal->io.Open(filename, adios2::Mode::Write, world); +#endif if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filename); singlefile_opened = 1; } diff --git a/src/ADIOS/reader_adios.cpp b/src/ADIOS/reader_adios.cpp index dd389be702..655df98d26 100644 --- a/src/ADIOS/reader_adios.cpp +++ b/src/ADIOS/reader_adios.cpp @@ -82,7 +82,11 @@ ReaderADIOS::ReaderADIOS(LAMMPS *lmp) : Reader(lmp) internal = new ReadADIOSInternal(); try { +#if defined(MPI_STUBS) + internal->ad = new adios2::ADIOS("adios2_config.xml", adios2::DebugON); +#else internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON); +#endif } catch (std::ios_base::failure &e) { error->one(FLERR, "ADIOS initialization failed with error: {}", e.what()); } @@ -134,7 +138,11 @@ void ReaderADIOS::open_file(const std::string &file) if (internal->fh) internal->fh.Close(); try { +#if defined(MPI_STUBS) + internal->fh = internal->io.Open(file, adios2::Mode::Read); +#else internal->fh = internal->io.Open(file, adios2::Mode::Read, world); +#endif } catch (std::ios_base::failure &e) { error->one(FLERR, "Error opening file {}: {}", file, e.what()); } From 401b5cee6d82bd07dc04a7c47bedaeac41b9ae4e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Nov 2022 07:32:08 -0500 Subject: [PATCH 35/58] add -y flag to add-repository commands to avoid issues with GPG support changes --- tools/singularity/ubuntu18.04.def | 2 +- tools/singularity/ubuntu18.04_amd_rocm.def | 2 +- tools/singularity/ubuntu18.04_gpu.def | 2 +- tools/singularity/ubuntu18.04_intel_opencl.def | 2 +- tools/singularity/ubuntu18.04_nvidia.def | 2 +- tools/singularity/ubuntu20.04.def | 2 +- tools/singularity/ubuntu20.04_amd_rocm.def | 2 +- tools/singularity/ubuntu20.04_gpu.def | 2 +- tools/singularity/ubuntu20.04_intel_opencl.def | 2 +- tools/singularity/ubuntu20.04_nvidia.def | 2 +- tools/singularity/ubuntu20.04_oneapi.def | 2 +- tools/singularity/ubuntu22.04.def | 2 +- 12 files changed, 12 insertions(+), 12 deletions(-) diff --git a/tools/singularity/ubuntu18.04.def b/tools/singularity/ubuntu18.04.def index 4f8dd22c43..02351d9ecb 100644 --- a/tools/singularity/ubuntu18.04.def +++ b/tools/singularity/ubuntu18.04.def @@ -5,7 +5,7 @@ From: ubuntu:18.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu18.04_amd_rocm.def b/tools/singularity/ubuntu18.04_amd_rocm.def index 721591c589..38eaa6e322 100644 --- a/tools/singularity/ubuntu18.04_amd_rocm.def +++ b/tools/singularity/ubuntu18.04_amd_rocm.def @@ -45,7 +45,7 @@ From: ubuntu:18.04 ########################################################################### # Common Software ########################################################################### - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get install --no-install-recommends -y \ bc \ diff --git a/tools/singularity/ubuntu18.04_gpu.def b/tools/singularity/ubuntu18.04_gpu.def index d6947bb742..cab62c623f 100644 --- a/tools/singularity/ubuntu18.04_gpu.def +++ b/tools/singularity/ubuntu18.04_gpu.def @@ -48,7 +48,7 @@ From: ubuntu:18.04 ########################################################################### # Common Software ########################################################################### - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get install --no-install-recommends -y \ bc \ diff --git a/tools/singularity/ubuntu18.04_intel_opencl.def b/tools/singularity/ubuntu18.04_intel_opencl.def index a1df05fe0e..2d562771bb 100644 --- a/tools/singularity/ubuntu18.04_intel_opencl.def +++ b/tools/singularity/ubuntu18.04_intel_opencl.def @@ -5,7 +5,7 @@ From: ubuntu:18.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu18.04_nvidia.def b/tools/singularity/ubuntu18.04_nvidia.def index 7c0819a7ab..2a3a34b1b2 100644 --- a/tools/singularity/ubuntu18.04_nvidia.def +++ b/tools/singularity/ubuntu18.04_nvidia.def @@ -5,7 +5,7 @@ From: nvidia/cuda:11.6.2-devel-ubuntu18.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu20.04.def b/tools/singularity/ubuntu20.04.def index 8b3c92ada0..2a2a1dd660 100644 --- a/tools/singularity/ubuntu20.04.def +++ b/tools/singularity/ubuntu20.04.def @@ -5,7 +5,7 @@ From: ubuntu:20.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu20.04_amd_rocm.def b/tools/singularity/ubuntu20.04_amd_rocm.def index 3e018363b2..f947de9ee9 100644 --- a/tools/singularity/ubuntu20.04_amd_rocm.def +++ b/tools/singularity/ubuntu20.04_amd_rocm.def @@ -36,7 +36,7 @@ From: ubuntu:20.04 ########################################################################### # Common Software ########################################################################### - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get install --no-install-recommends -y \ bc \ diff --git a/tools/singularity/ubuntu20.04_gpu.def b/tools/singularity/ubuntu20.04_gpu.def index 9c9dea571d..f84c1f8926 100644 --- a/tools/singularity/ubuntu20.04_gpu.def +++ b/tools/singularity/ubuntu20.04_gpu.def @@ -39,7 +39,7 @@ From: ubuntu:20.04 # Common Software ########################################################################### apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get install --no-install-recommends -y \ bc \ diff --git a/tools/singularity/ubuntu20.04_intel_opencl.def b/tools/singularity/ubuntu20.04_intel_opencl.def index 96769516d4..cc547fef29 100644 --- a/tools/singularity/ubuntu20.04_intel_opencl.def +++ b/tools/singularity/ubuntu20.04_intel_opencl.def @@ -5,7 +5,7 @@ From: ubuntu:20.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu20.04_nvidia.def b/tools/singularity/ubuntu20.04_nvidia.def index 654b69f0b9..8ec334ad8b 100644 --- a/tools/singularity/ubuntu20.04_nvidia.def +++ b/tools/singularity/ubuntu20.04_nvidia.def @@ -5,7 +5,7 @@ From: nvidia/cuda:11.6.2-devel-ubuntu20.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu20.04_oneapi.def b/tools/singularity/ubuntu20.04_oneapi.def index 8b66d8d441..25c913f392 100644 --- a/tools/singularity/ubuntu20.04_oneapi.def +++ b/tools/singularity/ubuntu20.04_oneapi.def @@ -5,7 +5,7 @@ From: ubuntu:20.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu22.04.def b/tools/singularity/ubuntu22.04.def index a1a2a45b49..8a4f262f75 100644 --- a/tools/singularity/ubuntu22.04.def +++ b/tools/singularity/ubuntu22.04.def @@ -5,7 +5,7 @@ From: ubuntu:22.04 export DEBIAN_FRONTEND=noninteractive apt-get update apt-get install --no-install-recommends -y software-properties-common gpg gpg-agent - add-apt-repository ppa:openkim/latest + add-apt-repository ppa:openkim/latest -y apt-get update apt-get upgrade --no-install-recommends -y apt-get install --no-install-recommends -y \ From dc4301dfa8d8d31f00f02c661dd1dab8e3df915e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Nov 2022 08:57:12 -0500 Subject: [PATCH 36/58] initialize ADIOS dumps only the first time when used in multiple runs --- src/ADIOS/dump_atom_adios.cpp | 96 ++++++++++++++++++----------------- 1 file changed, 49 insertions(+), 47 deletions(-) diff --git a/src/ADIOS/dump_atom_adios.cpp b/src/ADIOS/dump_atom_adios.cpp index deee23ba41..8e92300591 100644 --- a/src/ADIOS/dump_atom_adios.cpp +++ b/src/ADIOS/dump_atom_adios.cpp @@ -268,61 +268,63 @@ void DumpAtomADIOS::init_style() else if (scale_flag == 0 && image_flag == 1) pack_choice = &DumpAtomADIOS::pack_noscale_image; - /* Define the group of variables for the atom style here since it's a fixed - * set */ - internal->io = internal->ad->DeclareIO(internal->ioName); - if (!internal->io.InConfigFile()) { - // if not defined by user, we can change the default settings - // BPFile is the default writer - internal->io.SetEngine("BPFile"); - int num_aggregators = multiproc; - if (num_aggregators == 0) num_aggregators = 1; - auto nstreams = std::to_string(num_aggregators); - internal->io.SetParameters({{"substreams", nstreams}}); - if (me == 0) - utils::logmesg(lmp, "ADIOS method for {} is n-to-m (aggregation with {} writers)\n", filename, - nstreams); - } + /* Define the group of variables for the atom style here since it's a fixed set */ - internal->io.DefineVariable("ntimestep"); - internal->io.DefineVariable("natoms"); + if (!internal->io) { + internal->io = internal->ad->DeclareIO(internal->ioName); + if (!internal->io.InConfigFile()) { + // if not defined by user, we can change the default settings + // BPFile is the default writer + internal->io.SetEngine("BPFile"); + int num_aggregators = multiproc; + if (num_aggregators == 0) num_aggregators = 1; + auto nstreams = std::to_string(num_aggregators); + internal->io.SetParameters({{"substreams", nstreams}}); + if (me == 0) + utils::logmesg(lmp, "ADIOS method for {} is n-to-m (aggregation with {} writers)\n", filename, + nstreams); + } - internal->io.DefineVariable("nprocs"); - internal->io.DefineVariable("ncolumns"); + internal->io.DefineVariable("ntimestep"); + internal->io.DefineVariable("natoms"); - internal->io.DefineVariable("boxxlo"); - internal->io.DefineVariable("boxxhi"); - internal->io.DefineVariable("boxylo"); - internal->io.DefineVariable("boxyhi"); - internal->io.DefineVariable("boxzlo"); - internal->io.DefineVariable("boxzhi"); + internal->io.DefineVariable("nprocs"); + internal->io.DefineVariable("ncolumns"); - internal->io.DefineVariable("boxxy"); - internal->io.DefineVariable("boxxz"); - internal->io.DefineVariable("boxyz"); + internal->io.DefineVariable("boxxlo"); + internal->io.DefineVariable("boxxhi"); + internal->io.DefineVariable("boxylo"); + internal->io.DefineVariable("boxyhi"); + internal->io.DefineVariable("boxzlo"); + internal->io.DefineVariable("boxzhi"); - internal->io.DefineAttribute("triclinic", domain->triclinic); - internal->io.DefineAttribute("scaled", scale_flag); - internal->io.DefineAttribute("image", image_flag); + internal->io.DefineVariable("boxxy"); + internal->io.DefineVariable("boxxz"); + internal->io.DefineVariable("boxyz"); - int *boundaryptr = reinterpret_cast(domain->boundary); - internal->io.DefineAttribute("boundary", boundaryptr, 6); + internal->io.DefineAttribute("triclinic", domain->triclinic); + internal->io.DefineAttribute("scaled", scale_flag); + internal->io.DefineAttribute("image", image_flag); - auto nColumns = static_cast(size_one); - internal->io.DefineAttribute("columns", columnNames.data(), nColumns); - internal->io.DefineAttribute("columnstr", columns); - internal->io.DefineAttribute("boundarystr", boundstr); - internal->io.DefineAttribute("LAMMPS/dump_style", "atom"); - internal->io.DefineAttribute("LAMMPS/version", lmp->version); - internal->io.DefineAttribute("LAMMPS/num_ver", std::to_string(lmp->num_ver)); + int *boundaryptr = reinterpret_cast(domain->boundary); + internal->io.DefineAttribute("boundary", boundaryptr, 6); - // local dimension variables - internal->io.DefineVariable("nme", {adios2::LocalValueDim}); - internal->io.DefineVariable("offset", {adios2::LocalValueDim}); + auto nColumns = static_cast(size_one); + internal->io.DefineAttribute("columns", columnNames.data(), nColumns); + internal->io.DefineAttribute("columnstr", columns); + internal->io.DefineAttribute("boundarystr", boundstr); + internal->io.DefineAttribute("LAMMPS/dump_style", "atom"); + internal->io.DefineAttribute("LAMMPS/version", lmp->version); + internal->io.DefineAttribute("LAMMPS/num_ver", std::to_string(lmp->num_ver)); - // atom table size is not known at the moment - // it will be correctly defined at the moment of write - size_t UnknownSizeYet = 1; - internal->varAtoms = internal->io.DefineVariable( + // local dimension variables + internal->io.DefineVariable("nme", {adios2::LocalValueDim}); + internal->io.DefineVariable("offset", {adios2::LocalValueDim}); + + // atom table size is not known at the moment + // it will be correctly defined at the moment of write + size_t UnknownSizeYet = 1; + internal->varAtoms = internal->io.DefineVariable( "atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0}, {UnknownSizeYet, nColumns}); + } } From aa46f5560aa2760e6b5343cec8d215805a6d03b9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Nov 2022 22:57:32 -0500 Subject: [PATCH 37/58] add .. versiondded:: tag to BPM package commands and restrictions text consistent --- doc/src/bond_bpm_rotational.rst | 8 +++++--- doc/src/bond_bpm_spring.rst | 8 +++++--- doc/src/compute_nbond_atom.rst | 7 +++++-- doc/src/fix_nve_bpm_sphere.rst | 6 ++++++ doc/src/pair_bpm_spring.rst | 8 +++++++- 5 files changed, 28 insertions(+), 9 deletions(-) diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index 5c43071274..fb1db04b27 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -45,6 +45,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + The *bpm/rotational* bond style computes forces and torques based on deviations from the initial reference state of the two atoms. The reference state is stored by each bond when it is first computed in @@ -211,9 +213,9 @@ command, as *b1*, *b2*, ..., *b7*\ . Restrictions """""""""""" -This bond style can only be used if LAMMPS was built with the BPM -package. See the :doc:`Build package ` doc page for -more info. +This bond style is part of the BPM package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. By default if pair interactions are to be disabled, this bond style requires setting diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index de31357afe..46c300d364 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -45,6 +45,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + The *bpm/spring* bond style computes forces and torques based on deviations from the initial reference state of the two atoms. The reference state is stored by each bond when it is first computed in @@ -167,9 +169,9 @@ extra quantity can be accessed by the Restrictions """""""""""" -This bond style can only be used if LAMMPS was built with the BPM -package. See the :doc:`Build package ` doc page for -more info. +This bond style is part of the BPM package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. By default if pair interactions are to be disabled, this bond style requires setting diff --git a/doc/src/compute_nbond_atom.rst b/doc/src/compute_nbond_atom.rst index 0e53de3e49..f438836534 100644 --- a/doc/src/compute_nbond_atom.rst +++ b/doc/src/compute_nbond_atom.rst @@ -23,6 +23,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + Define a computation that computes the number of bonds each atom is part of. Bonds which are broken are not counted in the tally. See the :doc:`Howto broken bonds ` page for more information. @@ -40,8 +42,9 @@ LAMMPS output options. Restrictions """""""""""" -This fix can only be used if LAMMPS was built with the BPM package. -See the :doc:`Build package ` doc page for more info. +This compute is part of the BPM package. It is only enabled if LAMMPS was +built with that package. See the :doc:`Build package ` +page for more info. Related commands """""""""""""""" diff --git a/doc/src/fix_nve_bpm_sphere.rst b/doc/src/fix_nve_bpm_sphere.rst index 861586ab2a..ef170605a4 100644 --- a/doc/src/fix_nve_bpm_sphere.rst +++ b/doc/src/fix_nve_bpm_sphere.rst @@ -30,6 +30,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + Perform constant NVE integration to update position, velocity, angular velocity, and quaternion orientation for finite-size spherical particles in the group each timestep. V is volume; E is energy. This @@ -65,6 +67,10 @@ the :doc:`run ` command. This fix is not invoked during Restrictions """""""""""" +This fix is part of the BPM package. It is only enabled if LAMMPS was +built with that package. See the :doc:`Build package ` +page for more info. + This fix requires that atoms store torque, angular velocity (omega), a radius, and a quaternion as defined by the :doc:`atom_style bpm/sphere ` command. diff --git a/doc/src/pair_bpm_spring.rst b/doc/src/pair_bpm_spring.rst index 72e0083bd8..7dfa9bc12c 100644 --- a/doc/src/pair_bpm_spring.rst +++ b/doc/src/pair_bpm_spring.rst @@ -22,6 +22,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + Style *bpm/spring* computes pairwise forces with the formula .. math:: @@ -101,7 +103,11 @@ This pair style can only be used via the *pair* keyword of the Restrictions """""""""""" - none + +This pair style is part of the BPM package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + Related commands """""""""""""""" From d53ce7aba990183a3de802ff72f1837ad90c83a2 Mon Sep 17 00:00:00 2001 From: mehdibghk Date: Wed, 9 Nov 2022 11:15:18 +0800 Subject: [PATCH 38/58] initialize pointer to null --- src/ASPHERE/pair_ylz.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ASPHERE/pair_ylz.cpp b/src/ASPHERE/pair_ylz.cpp index 2d035b6f5d..afedf635d5 100644 --- a/src/ASPHERE/pair_ylz.cpp +++ b/src/ASPHERE/pair_ylz.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Mehdi Baghaee (SUSTech) + Contributing author: Hongyan Yuan (SUSTech) ------------------------------------------------------------------------- */ #include "pair_ylz.h" From f9d07d8932277edceea4eb97d879d9cc9d41c2c3 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 9 Nov 2022 17:40:57 -0500 Subject: [PATCH 39/58] 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 40/58] 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 4392b9c8cbc0c2b09827074509e9cbc9419d86a8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 02:24:35 -0500 Subject: [PATCH 41/58] store LAMMPS version of restart, if initialized from restart file --- src/lammps.cpp | 3 +++ src/lammps.h | 2 ++ src/read_restart.cpp | 2 ++ unittest/formats/test_file_operations.cpp | 5 +++++ 4 files changed, 12 insertions(+) diff --git a/src/lammps.cpp b/src/lammps.cpp index 9b9dda4ac8..5e9fed61b7 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -130,6 +130,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) : version = (const char *) LAMMPS_VERSION; num_ver = utils::date2num(version); + restart_ver = -1; external_comm = 0; mdicomm = nullptr; @@ -993,6 +994,8 @@ void LAMMPS::destroy() delete python; python = nullptr; + + restart_ver = -1; // reset last restart version id } /* ---------------------------------------------------------------------- diff --git a/src/lammps.h b/src/lammps.h index abe08b45f8..03b019ed43 100644 --- a/src/lammps.h +++ b/src/lammps.h @@ -49,6 +49,8 @@ class LAMMPS { // that is constructed so that will be greater // for newer versions in numeric or string // value comparisons + int restart_ver; // -1 or numeric version id of LAMMPS version in restart + // file, in case LAMMPS was initialized from a restart // MPI_Comm world; // MPI communicator FILE *infile; // infile diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 35d4629291..d4b455aaf9 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -599,6 +599,8 @@ void ReadRestart::header() if (flag == VERSION) { char *version = read_string(); + lmp->restart_ver = utils::date2num(version); + if (me == 0) utils::logmesg(lmp," restart file = {}, LAMMPS = {}\n", version, lmp->version); delete[] version; diff --git a/unittest/formats/test_file_operations.cpp b/unittest/formats/test_file_operations.cpp index 6149b0c99c..9c75029130 100644 --- a/unittest/formats/test_file_operations.cpp +++ b/unittest/formats/test_file_operations.cpp @@ -303,6 +303,7 @@ TEST_F(FileOperationsTest, error_all_one) TEST_F(FileOperationsTest, write_restart) { + ASSERT_EQ(lmp->restart_ver, -1); BEGIN_HIDE_OUTPUT(); command("echo none"); END_HIDE_OUTPUT(); @@ -356,6 +357,7 @@ TEST_F(FileOperationsTest, write_restart) BEGIN_HIDE_OUTPUT(); command("clear"); END_HIDE_OUTPUT(); + ASSERT_EQ(lmp->restart_ver, -1); ASSERT_EQ(lmp->atom->natoms, 0); ASSERT_EQ(lmp->update->ntimestep, 0); ASSERT_EQ(lmp->domain->triclinic, 0); @@ -369,18 +371,21 @@ TEST_F(FileOperationsTest, write_restart) command("change_box all triclinic"); command("write_restart triclinic.restart"); END_HIDE_OUTPUT(); + ASSERT_EQ(lmp->restart_ver, lmp->num_ver); ASSERT_EQ(lmp->atom->natoms, 1); ASSERT_EQ(lmp->update->ntimestep, 333); ASSERT_EQ(lmp->domain->triclinic, 1); BEGIN_HIDE_OUTPUT(); command("clear"); END_HIDE_OUTPUT(); + ASSERT_EQ(lmp->restart_ver, -1); ASSERT_EQ(lmp->atom->natoms, 0); ASSERT_EQ(lmp->update->ntimestep, 0); ASSERT_EQ(lmp->domain->triclinic, 0); BEGIN_HIDE_OUTPUT(); command("read_restart triclinic.restart"); END_HIDE_OUTPUT(); + ASSERT_EQ(lmp->restart_ver, lmp->num_ver); ASSERT_EQ(lmp->atom->natoms, 1); ASSERT_EQ(lmp->update->ntimestep, 333); ASSERT_EQ(lmp->domain->triclinic, 1); From c00a5d52d21cc3eb6af7dee8786a1658c3456f9b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 02:25:00 -0500 Subject: [PATCH 42/58] silence compiler warnings --- src/AMOEBA/amoeba_polar.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp index a0f4fce301..9342b6213d 100644 --- a/src/AMOEBA/amoeba_polar.cpp +++ b/src/AMOEBA/amoeba_polar.cpp @@ -250,10 +250,12 @@ void PairAmoeba::polar_real() double drc3[3],drc5[3],drc7[3]; double urc3[3],urc5[3],tep[3]; double fix[3],fiy[3],fiz[3]; +#if 0 // for poltyp TCG which is currently not supported double uax[3],uay[3],uaz[3]; double ubx[3],uby[3],ubz[3]; double uaxp[3],uayp[3],uazp[3]; double ubxp[3],ubyp[3],ubzp[3]; +#endif double dmpi[9],dmpk[9]; double dmpik[9]; double bn[5]; @@ -320,6 +322,7 @@ void PairAmoeba::polar_real() uixp = uinp[i][0]; uiyp = uinp[i][1]; uizp = uinp[i][2]; +#if 0 // for poltyp TCG which is currently not supported for (m = 0; m < tcgnab; m++) { uax[m] = uad[m][i][0]; uay[m] = uad[m][i][1]; @@ -334,7 +337,7 @@ void PairAmoeba::polar_real() ubyp[m] = ubp[m][i][1]; ubzp[m] = ubp[m][i][2]; } - +#endif if (amoeba) { pdi = pdamp[itype]; pti = thole[itype]; @@ -1030,7 +1033,8 @@ void PairAmoeba::polar_real() // get the dtau/dr terms used for TCG polarization force } else if (poltyp == TCG) { - /* +#if 0 + // poltyp TCG not yet supported for AMOEBA/HIPPO for (m = 0; m < tcgnab; m++) { ukx = ubd[m][j][0]; uky = ubd[m][j][1]; @@ -1131,8 +1135,7 @@ void PairAmoeba::polar_real() frcx += depx; frcy += depy; frcz += depz; - } - */ +#endif } // increment force-based gradient on the interaction sites From 1932b6390a04b982b9a19e8c59beef661986b26a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 06:51:43 -0500 Subject: [PATCH 43/58] update and sort codeowners lists --- .github/CODEOWNERS | 60 ++++++++++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 26 deletions(-) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 4e8803fc5a..b856b3a8dd 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -13,40 +13,44 @@ lib/kim/* @ellio167 lib/mesont/* @iafoss # whole packages +src/ADIOS/* @pnorbert src/AMOEBA/* @sjplimp -src/COMPRESS/* @rbberger -src/GPU/* @ndtrung81 -src/KOKKOS/* @stanmoore1 -src/KIM/* @ellio167 -src/LATTE/* @cnegre -src/MLIAP/* @athomps -src/SNAP/* @athomps -src/SPIN/* @julient31 +src/BPM/* @jtclemm src/BROWNIAN/* @samueljmcameron src/CG-DNA/* @ohenrich src/CG-SPICA/* @yskmiyazaki src/COLVARS/* @giacomofiorin +src/COMPRESS/* @rbberger src/DIELECTRIC/* @ndtrung81 src/ELECTRODE/* @ludwig-ahrens src/FEP/* @agiliopadua -src/ML-HDNNP/* @singraber +src/GPU/* @ndtrung81 +src/GRANULAR/* @jtclemm @dsbolin src/INTEL/* @wmbrownintel +src/KIM/* @ellio167 +src/KOKKOS/* @stanmoore1 +src/LATTE/* @cnegre src/MANIFOLD/* @Pakketeretet2 -src/MDI/* @taylor-a-barnes +src/MDI/* @taylor-a-barnes @sjplimp src/MEAM/* @martok src/MESONT/* @iafoss +src/ML-HDNNP/* @singraber +src/ML-IAP/* @athomps +src/ML-PACE/* @yury-lysogorskiy src/MOFFF/* @hheenen src/MOLFILE/* @akohlmey src/NETCDF/* @pastewka -src/ML-PACE/* @yury-lysogorskiy -src/PLUMED/* @gtribello -src/PHONON/* @lingtikong -src/PTM/* @pmla src/OPENMP/* @akohlmey +src/PHONON/* @lingtikong +src/PLUGIN/* @akohlmey +src/PLUMED/* @gtribello +src/PTM/* @pmla src/QMMM/* @akohlmey -src/REAXFF/* @hasanmetin @stanmoore1 src/REACTION/* @jrgissing +src/REAXFF/* @hasanmetin @stanmoore1 src/SCAFACOS/* @rhalver +src/SNAP/* @athomps +src/SPIN/* @julient31 src/TALLY/* @akohlmey src/UEF/* @danicholson src/VTK/* @rbberger @@ -120,27 +124,32 @@ src/dump_movie.* @akohlmey src/exceptions.h @rbberger src/fix_nh.* @athomps src/info.* @akohlmey @rbberger -src/timer.* @akohlmey src/min* @sjplimp @stanmoore1 +src/platform.* @akohlmey +src/timer.* @akohlmey src/utils.* @akohlmey @rbberger src/verlet.* @sjplimp @stanmoore1 src/math_eigen_impl.h @jewettaij # tools -tools/msi2lmp/* @akohlmey +tools/coding_standard/* @akohlmey @rbberger tools/emacs/* @HaoZeke -tools/singularity/* @akohlmey @rbberger -tools/coding_standard/* @rbberger -tools/valgrind/* @akohlmey -tools/swig/* @akohlmey +tools/lammps-shell/* @akohlmey +tools/msi2lmp/* @akohlmey tools/offline/* @rbberger +tools/singularity/* @akohlmey @rbberger +tools/swig/* @akohlmey +tools/valgrind/* @akohlmey tools/vim/* @hammondkd # tests -unittest/* @akohlmey @rbberger +unittest/* @akohlmey # cmake cmake/* @junghans @rbberger +cmake/Modules/LAMMPSInterfacePlugin.cmake @akohlmey +cmake/Modules/MPI4WIN.cmake @akohlmey +cmake/Modules/OpenCLLoader.cmake @akohlmey cmake/Modules/Packages/COLVARS.cmake @junghans @rbberger @giacomofiorin cmake/Modules/Packages/KIM.cmake @junghans @rbberger @ellio167 cmake/presets/*.cmake @akohlmey @@ -149,13 +158,12 @@ cmake/presets/*.cmake @akohlmey python/* @rbberger # fortran -fortran/* @akohlmey +fortran/* @akohlmey @hammondkd # docs -doc/utils/*/* @rbberger -doc/Makefile @rbberger -doc/README @rbberger +doc/* @akohlmey examples/plugin/* @akohlmey +examples/PACKAGES/pace/plugin/* @akohlmey # for releases src/version.h @sjplimp From 01b4600ba5f47c553bfda5ca85f22d88a7a50ce5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 07:04:57 -0500 Subject: [PATCH 44/58] remove obsolete file --- cmake/Modules/FindZMQ.cmake | 19 ------------------- 1 file changed, 19 deletions(-) delete mode 100644 cmake/Modules/FindZMQ.cmake diff --git a/cmake/Modules/FindZMQ.cmake b/cmake/Modules/FindZMQ.cmake deleted file mode 100644 index 7d612c2eb3..0000000000 --- a/cmake/Modules/FindZMQ.cmake +++ /dev/null @@ -1,19 +0,0 @@ -find_path(ZMQ_INCLUDE_DIR zmq.h) -find_library(ZMQ_LIBRARY NAMES zmq) - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(ZMQ DEFAULT_MSG ZMQ_LIBRARY ZMQ_INCLUDE_DIR) - -# Copy the results to the output variables and target. -if(ZMQ_FOUND) - set(ZMQ_LIBRARIES ${ZMQ_LIBRARY}) - set(ZMQ_INCLUDE_DIRS ${ZMQ_INCLUDE_DIR}) - - if(NOT TARGET ZMQ::ZMQ) - add_library(ZMQ::ZMQ UNKNOWN IMPORTED) - set_target_properties(ZMQ::ZMQ PROPERTIES - IMPORTED_LINK_INTERFACE_LANGUAGES "C" - IMPORTED_LOCATION "${ZMQ_LIBRARY}" - INTERFACE_INCLUDE_DIRECTORIES "${ZMQ_INCLUDE_DIR}") - endif() -endif() From 34a509322930d4db0fab3415affd3a8905231468 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 16:03:06 -0500 Subject: [PATCH 45/58] refactor handling of the python source command. document it and more limits. --- doc/src/python.rst | 236 +++++++++++++----------- src/PYTHON/python_impl.cpp | 18 +- unittest/python/test_python_package.cpp | 4 +- 3 files changed, 141 insertions(+), 117 deletions(-) diff --git a/doc/src/python.rst b/doc/src/python.rst index aad2f636d3..d777719181 100644 --- a/doc/src/python.rst +++ b/doc/src/python.rst @@ -8,14 +8,24 @@ Syntax .. parsed-literal:: - python func keyword args ... + python function keyword args ... -* func = name of Python function -* one or more keyword/args pairs must be appended +* function = *source* or name of Python function + + if function is *source*: .. parsed-literal:: - keyword = *invoke* or *input* or *return* or *format* or *length* or *file* or *here* or *exists* or *source* + keyword = *inline* or name of a Python file + inline = one or more lines of Python code which will be executed immediately + must be a single argument, typically enclosed in triple quotes + Python file = name of a file with Python code which will be executed immediately + +* if function is the name of a Python function, one or more keyword/args pairs must be appended + + .. parsed-literal:: + + keyword = *invoke* or *input* or *return* or *format* or *length* or *file* or *here* or *exists* *invoke* arg = none = invoke the previously defined Python function *input* args = N i1 i2 ... iN N = # of inputs to function @@ -38,10 +48,6 @@ Syntax inline = one or more lines of Python code which defines func must be a single argument, typically enclosed between triple quotes *exists* arg = none = Python code has been loaded by previous python command - *source* arg = *filename* or *inline* - filename = file of Python code which will be executed immediately - inline = one or more lines of Python code which will be executed immediately - must be a single argument, typically enclosed between triple quotes Examples """""""" @@ -70,80 +76,89 @@ Examples lmp.command("pair_style lj/cut ${cut}") # LAMMPS commands lmp.command("pair_coeff * * 1.0 1.0") lmp.command("run 100") - """ + """ + + python source funcdef.py + + python source inline "from lammps import lammps" + Description """"""""""" Define a Python function or execute a previously defined function or -execute some arbitrary python code. -Arguments, including LAMMPS variables, can be passed to the function -from the LAMMPS input script and a value returned by the Python -function to a LAMMPS variable. The Python code for the function can -be included directly in the input script or in a separate Python file. -The function can be standard Python code or it can make "callbacks" to -LAMMPS through its library interface to query or set internal values -within LAMMPS. This is a powerful mechanism for performing complex -operations in a LAMMPS input script that are not possible with the -simple input script and variable syntax which LAMMPS defines. Thus -your input script can operate more like a true programming language. +execute some arbitrary python code. Arguments, including LAMMPS +variables, can be passed to the function from the LAMMPS input script +and a value returned by the Python function to a LAMMPS variable. The +Python code for the function can be included directly in the input +script or in a separate Python file. The function can be standard +Python code or it can make "callbacks" to LAMMPS through its library +interface to query or set internal values within LAMMPS. This is a +powerful mechanism for performing complex operations in a LAMMPS input +script that are not possible with the simple input script and variable +syntax which LAMMPS defines. Thus your input script can operate more +like a true programming language. Use of this command requires building LAMMPS with the PYTHON package which links to the Python library so that the Python interpreter is embedded in LAMMPS. More details about this process are given below. -There are two ways to invoke a Python function once it has been -defined. One is using the *invoke* keyword. The other is to assign -the function to a :doc:`python-style variable ` defined in -your input script. Whenever the variable is evaluated, it will -execute the Python function to assign a value to the variable. Note -that variables can be evaluated in many different ways within LAMMPS. -They can be substituted for directly in an input script. Or they can -be passed to various commands as arguments, so that the variable is -evaluated during a simulation run. +There are two ways to invoke a Python function once it has been defined. +One is using the *invoke* keyword. The other is to assign the function +to a :doc:`python-style variable ` defined in your input +script. Whenever the variable is evaluated, it will execute the Python +function to assign a value to the variable. Note that variables can be +evaluated in many different ways within LAMMPS. They can be substituted +with their result directly in an input script, or they can be passed to +various commands as arguments, so that the variable is evaluated during +a simulation run. -A broader overview of how Python can be used with LAMMPS is given on -the :doc:`Python ` doc page. There is an examples/python -directory which illustrates use of the python command. +A broader overview of how Python can be used with LAMMPS is given in the +:doc:`Use Python with LAMMPS ` section of the +documentation. There is an ``examples/python`` directory which +illustrates use of the python command. ---------- -The *func* setting specifies the name of the Python function. The -code for the function is defined using the *file* or *here* keywords -as explained below. In case of the *source* keyword, the name of -the function is ignored. +The first argument of the *python* command is either the *source* +keyword or the name of a Python function. + +If the *source* keyword is used, no other keywords can be used. The +argument either can be a filename or the keyword *inline* followed by a +string with python commands, either on a single line enclosed in quotes, +or as multiple lines enclosed in triple quotes. These python commands +will be passed to the python interpreter and executed immediately +without registering a python function for future execution. + +In all other cases, the first argument is the name of a Python function +that will be registered with LAMMPS for future execution. The function +may already be defined (see *exists* keyword) or must be defined using +the *file* or *here* keywords as explained below. If the *invoke* keyword is used, no other keywords can be used, and a -previous python command must have defined the Python function +previous python command must have registered the Python function referenced by this command. This invokes the Python function with the -previously defined arguments and return value processed as explained -below. You can invoke the function as many times as you wish in your -input script. - -If the *source* keyword is used, no other keywords can be used. -The argument can be a filename or a string with python commands, -either on a single line enclosed in quotes, or as multiple lines -enclosed in triple quotes. These python commands will be passed -to the python interpreter and executed immediately without registering -a python function for future execution. +previously defined arguments and the return value is processed as +explained below. You can invoke the function as many times as you wish +in your input script. The *input* keyword defines how many arguments *N* the Python function -expects. If it takes no arguments, then the *input* keyword should -not be used. Each argument can be specified directly as a value, -e.g. 6 or 3.14159 or abc (a string of characters). The type of each -argument is specified by the *format* keyword as explained below, so -that Python will know how to interpret the value. If the word SELF is -used for an argument it has a special meaning. A pointer is passed to -the Python function which it converts into a reference to LAMMPS -itself. This enables the function to call back to LAMMPS through its -library interface as explained below. This allows the Python function -to query or set values internal to LAMMPS which can affect the -subsequent execution of the input script. A LAMMPS variable can also -be used as an argument, specified as v_name, where "name" is the name -of the variable. Any style of LAMMPS variable can be used, as defined -by the :doc:`variable ` command. Each time the Python -function is invoked, the LAMMPS variable is evaluated and its value is -passed to the Python function. +expects. If it takes no arguments, then the *input* keyword should not +be used. Each argument can be specified directly as a value, e.g. 6 or +3.14159 or abc (a string of characters). The type of each argument is +specified by the *format* keyword as explained below, so that Python +will know how to interpret the value. If the word SELF is used for an +argument it has a special meaning. A pointer is passed to the Python +function which it converts into a reference to LAMMPS itself. This +enables the function to call back to LAMMPS through its library +interface as explained below. This allows the Python function to query +or set values internal to LAMMPS which can affect the subsequent +execution of the input script. A LAMMPS variable can also be used as an +argument, specified as v_name, where "name" is the name of the variable. +Any style of LAMMPS variable can be used, as defined by the +:doc:`variable ` command. Each time the Python function is +invoked, the LAMMPS variable is evaluated and its value is passed to the +Python function. The *return* keyword is only needed if the Python function returns a value. The specified *varReturn* must be of the form v_name, where @@ -165,19 +180,19 @@ The two commands can appear in either order in the input script so long as both are specified before the Python function is invoked for the first time. -The *format* keyword must be used if the *input* or *return* keyword -is used. It defines an *fstring* with M characters, where M = sum of +The *format* keyword must be used if the *input* or *return* keyword is +used. It defines an *fstring* with M characters, where M = sum of number of inputs and outputs. The order of characters corresponds to the N inputs, followed by the return value (if it exists). Each character must be one of the following: "i" for integer, "f" for -floating point, "s" for string, or "p" for SELF. Each character -defines the type of the corresponding input or output value of the -Python function and affects the type conversion that is performed -internally as data is passed back and forth between LAMMPS and Python. -Note that it is permissible to use a :doc:`python-style variable ` in a LAMMPS command that allows for an -equal-style variable as an argument, but only if the output of the -Python function is flagged as a numeric value ("i" or "f") via the -*format* keyword. +floating point, "s" for string, or "p" for SELF. Each character defines +the type of the corresponding input or output value of the Python +function and affects the type conversion that is performed internally as +data is passed back and forth between LAMMPS and Python. Note that it +is permissible to use a :doc:`python-style variable ` in a +LAMMPS command that allows for an equal-style variable as an argument, +but only if the output of the Python function is flagged as a numeric +value ("i" or "f") via the *format* keyword. If the *return* keyword is used and the *format* keyword specifies the output as a string, then the default maximum length of that string is @@ -192,12 +207,12 @@ truncated. Either the *file*, *here*, or *exists* keyword must be used, but only one of them. These keywords specify what Python code to load into the -Python interpreter. The *file* keyword gives the name of a file, -which should end with a ".py" suffix, which contains Python code. The -code will be immediately loaded into and run in the "main" module of -the Python interpreter. Note that Python code which contains a -function definition does not "execute" the function when it is run; it -simply defines the function so that it can be invoked later. +Python interpreter. The *file* keyword gives the name of a file +containing Python code, which should end with a ".py" suffix. The code +will be immediately loaded into and run in the "main" module of the +Python interpreter. Note that Python code which contains a function +definition does not "execute" the function when it is run; it simply +defines the function so that it can be invoked later. The *here* keyword does the same thing, except that the Python code follows as a single argument to the *here* keyword. This can be done @@ -208,14 +223,15 @@ proper indentation, blank lines, and comments, as desired. See the how triple quotes can be used as part of input script syntax. The *exists* keyword takes no argument. It means that Python code -containing the required Python function defined by the *func* setting, -is assumed to have been previously loaded by another python command. +containing the required Python function with the given name has already +been executed, for example by a *python source* command or in the same +file that was used previously with the *file* keyword. -Note that the Python code that is loaded and run must contain a -function with the specified *func* name. To operate properly when -later invoked, the function code must match the *input* and -*return* and *format* keywords specified by the python command. -Otherwise Python will generate an error. +Note that the Python code that is loaded and run must contain a function +with the specified function name. To operate properly when later +invoked, the function code must match the *input* and *return* and +*format* keywords specified by the python command. Otherwise Python +will generate an error. ---------- @@ -225,19 +241,18 @@ LAMMPS. Whether you load Python code from a file or directly from your input script, via the *file* and *here* keywords, the code can be identical. It must be indented properly as Python requires. It can contain -comments or blank lines. If the code is in your input script, it -cannot however contain triple-quoted Python strings, since that will -conflict with the triple-quote parsing that the LAMMPS input script -performs. +comments or blank lines. If the code is in your input script, it cannot +however contain triple-quoted Python strings, since that will conflict +with the triple-quote parsing that the LAMMPS input script performs. All the Python code you specify via one or more python commands is loaded into the Python "main" module, i.e. __main__. The code can define global variables or statements that are outside of function definitions. It can contain multiple functions, only one of which matches the *func* setting in the python command. This means you can -use the *file* keyword once to load several functions, and the -*exists* keyword thereafter in subsequent python commands to access -the other functions previously loaded. +use the *file* keyword once to load several functions, and the *exists* +keyword thereafter in subsequent python commands to access the other +functions previously loaded. A Python function you define (or more generally, the code you load) can import other Python modules or classes, it can make calls to other @@ -495,24 +510,35 @@ Restrictions """""""""""" This command is part of the PYTHON package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. -Building LAMMPS with the PYTHON package will link LAMMPS with the -Python library on your system. Settings to enable this are in the +Building LAMMPS with the PYTHON package will link LAMMPS with the Python +library on your system. Settings to enable this are in the lib/python/Makefile.lammps file. See the lib/python/README file for information on those settings. -If you use Python code which calls back to LAMMPS, via the SELF input argument -explained above, there is an extra step required when building LAMMPS. LAMMPS -must also be built as a shared library and your Python function must be able to -load the :doc:`"lammps" Python module ` that wraps the LAMMPS -library interface. These are the same steps required to use Python by itself -to wrap LAMMPS. Details on these steps are explained on the :doc:`Python -` doc page. Note that it is important that the stand-alone LAMMPS -executable and the LAMMPS shared library be consistent (built from the same -source code files) in order for this to work. If the two have been built at +If you use Python code which calls back to LAMMPS, via the SELF input +argument explained above, there is an extra step required when building +LAMMPS. LAMMPS must also be built as a shared library and your Python +function must be able to load the :doc:`"lammps" Python module +` that wraps the LAMMPS library interface. These are the +same steps required to use Python by itself to wrap LAMMPS. Details on +these steps are explained on the :doc:`Python ` doc page. +Note that it is important that the stand-alone LAMMPS executable and the +LAMMPS shared library be consistent (built from the same source code +files) in order for this to work. If the two have been built at different times using different source files, problems may occur. +Another limitation of calling back to Python from the LAMMPS module +using the *python* command in a LAMMPS input is that both, the Python +interpreter and LAMMPS, must be linked to the same Python runtime as a +shared library. If the Python interpreter is linked to Python +statically (which seems to happen with Conda) then loading the shared +LAMMPS library will create a second python "main" module that hides the +one from the Python interpreter and all previous defined function and +global variables will become invisible. + Related commands """""""""""""""" diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index 92668674d1..f7c9684615 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -130,19 +130,17 @@ void PythonImpl::command(int narg, char **arg) return; } - // if source is only keyword, execute the python code + // if source is only keyword, execute the python code in file - if (narg == 3 && strcmp(arg[1], "source") == 0) { - int err; + if ((narg > 1) && (strcmp(arg[0], "source") == 0)) { + int err = -1; - FILE *fp = fopen(arg[2], "r"); - if (fp == nullptr) + if ((narg > 2) && (strcmp(arg[1], "inline") == 0)) { err = execute_string(arg[2]); - else - err = execute_file(arg[2]); - - if (fp) fclose(fp); - if (err) error->all(FLERR, "Could not process Python source command"); + } else { + if (platform::file_is_readable(arg[1])) err = execute_file(arg[1]); + } + if (err) error->warning(FLERR, "Could not process Python source command. Error code: {}", err); return; } diff --git a/unittest/python/test_python_package.cpp b/unittest/python/test_python_package.cpp index 4e5aa53b0c..8b4f998bfc 100644 --- a/unittest/python/test_python_package.cpp +++ b/unittest/python/test_python_package.cpp @@ -276,7 +276,7 @@ TEST_F(PythonPackageTest, RunSource) { // execute python script from file auto output = CAPTURE_OUTPUT([&] { - command("python xyz source ${input_dir}/run.py"); + command("python source ${input_dir}/run.py"); }); ASSERT_THAT(output, HasSubstr(LOREM_IPSUM)); @@ -286,7 +286,7 @@ TEST_F(PythonPackageTest, RunSourceInline) { // execute inline python script auto output = CAPTURE_OUTPUT([&] { - command("python xyz source \"\"\"\n" + command("python source inline \"\"\"\n" "from __future__ import print_function\n" "print(2+2)\n" "\"\"\""); From b6b81a951aeee743d3886177e171f1d3cd7f5471 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 17:17:23 -0500 Subject: [PATCH 46/58] improve error reporting for python style variables --- src/PYTHON/python_impl.cpp | 12 ++++++------ src/variable.cpp | 17 +++++++++++++++-- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index f7c9684615..af16d2866a 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -112,18 +112,18 @@ PythonImpl::~PythonImpl() void PythonImpl::command(int narg, char **arg) { - if (narg < 2) error->all(FLERR, "Invalid python command"); + if (narg < 2) utils::missing_cmd_args(FLERR, "python", error); // if invoke is only keyword, invoke the previously defined function if (narg == 2 && strcmp(arg[1], "invoke") == 0) { int ifunc = find(arg[0]); - if (ifunc < 0) error->all(FLERR, "Python invoke of undefined function"); + if (ifunc < 0) error->all(FLERR, "Python invoke of undefined function: {}", arg[0]); char *str = nullptr; if (pfuncs[ifunc].noutput) { str = input->variable->pythonstyle(pfuncs[ifunc].ovarname, pfuncs[ifunc].name); - if (!str) error->all(FLERR, "Python variable does not match Python function"); + if (!str) error->all(FLERR, "Python variable {} does not match variable {} registered with Python function {}",arg[0], pfuncs[ifunc].ovarname, pfuncs[ifunc].name); } invoke_function(ifunc, str); @@ -363,9 +363,9 @@ int PythonImpl::variable_match(const char *name, const char *varname, int numeri { int ifunc = find(name); if (ifunc < 0) return -1; - if (pfuncs[ifunc].noutput == 0) return -1; - if (strcmp(pfuncs[ifunc].ovarname, varname) != 0) return -1; - if (numeric && pfuncs[ifunc].otype == STRING) return -1; + if (pfuncs[ifunc].noutput == 0) return -2; + if (strcmp(pfuncs[ifunc].ovarname, varname) != 0) return -3; + if (numeric && pfuncs[ifunc].otype == STRING) return -4; return ifunc; } diff --git a/src/variable.cpp b/src/variable.cpp index 819f130a02..bb45649208 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -985,8 +985,21 @@ char *Variable::retrieve(const char *name) str = data[ivar][1] = utils::strdup(result); } else if (style[ivar] == PYTHON) { int ifunc = python->variable_match(data[ivar][0],name,0); - if (ifunc < 0) - error->all(FLERR,"Python variable {} does not match Python function {}", name, data[ivar][0]); + if (ifunc < 0) { + if (ifunc == -1) { + error->all(FLERR, "Could not find Python function {} linked to variable {}", + data[ivar][0], name); + } else if (ifunc == -2) { + error->all(FLERR, "Python function {} for variable {} does not have a return value", + data[ivar][0], name); + } else if (ifunc == -3) { + error->all(FLERR,"Python variable {} does not match variable name registered with " + "Python function {}", name, data[ivar][0]); + } else { + error->all(FLERR, "Unknown error verifying function {} linked to python style variable {}", + data[ivar][0],name); + } + } python->invoke_function(ifunc,data[ivar][1]); str = data[ivar][1]; // if Python func returns a string longer than VALUELENGTH From 148df8589bf61b4e6bf019e9ceb9c1f58c5197bc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Nov 2022 23:30:03 -0500 Subject: [PATCH 47/58] revise error reporting in the python command --- src/PYTHON/python_impl.cpp | 87 +++++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 34 deletions(-) diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index af16d2866a..cacf1ca237 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -118,12 +118,16 @@ void PythonImpl::command(int narg, char **arg) if (narg == 2 && strcmp(arg[1], "invoke") == 0) { int ifunc = find(arg[0]); - if (ifunc < 0) error->all(FLERR, "Python invoke of undefined function: {}", arg[0]); + if (ifunc < 0) error->all(FLERR, "Python invoke of unknown function: {}", arg[0]); char *str = nullptr; if (pfuncs[ifunc].noutput) { str = input->variable->pythonstyle(pfuncs[ifunc].ovarname, pfuncs[ifunc].name); - if (!str) error->all(FLERR, "Python variable {} does not match variable {} registered with Python function {}",arg[0], pfuncs[ifunc].ovarname, pfuncs[ifunc].name); + if (!str) + error->all(FLERR, + "Python variable {} does not match variable {} " + "registered with Python function {}", + arg[0], pfuncs[ifunc].ovarname, pfuncs[ifunc].name); } invoke_function(ifunc, str); @@ -138,9 +142,12 @@ void PythonImpl::command(int narg, char **arg) if ((narg > 2) && (strcmp(arg[1], "inline") == 0)) { err = execute_string(arg[2]); } else { - if (platform::file_is_readable(arg[1])) err = execute_file(arg[1]); + if (platform::file_is_readable(arg[1])) + err = execute_file(arg[1]); + else + error->all(FLERR, "Could not open python source file {} for processing", arg[1]); } - if (err) error->warning(FLERR, "Could not process Python source command. Error code: {}", err); + if (err) error->all(FLERR, "Failure in python source command"); return; } @@ -160,48 +167,51 @@ void PythonImpl::command(int narg, char **arg) int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg], "input") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python input", error); ninput = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); - if (ninput < 0) error->all(FLERR, "Invalid python command"); + if (ninput < 0) error->all(FLERR, "Invalid number of python input arguments: {}", ninput); iarg += 2; delete[] istr; istr = new char *[ninput]; - if (iarg + ninput > narg) error->all(FLERR, "Invalid python command"); + if (iarg + ninput > narg) utils::missing_cmd_args(FLERR, "python input", error); for (int i = 0; i < ninput; i++) istr[i] = arg[iarg + i]; iarg += ninput; } else if (strcmp(arg[iarg], "return") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python return", error); noutput = 1; ostr = arg[iarg + 1]; iarg += 2; } else if (strcmp(arg[iarg], "format") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python format", error); format = utils::strdup(arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "length") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python length", error); length_longstr = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); - if (length_longstr <= 0) error->all(FLERR, "Invalid python command"); + if (length_longstr <= 0) error->all(FLERR, "Invalid python return value length"); iarg += 2; } else if (strcmp(arg[iarg], "file") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python file", error); delete[] pyfile; pyfile = utils::strdup(arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "here") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Invalid python command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "python here", error); herestr = arg[iarg + 1]; iarg += 2; } else if (strcmp(arg[iarg], "exists") == 0) { existflag = 1; iarg++; } else - error->all(FLERR, "Invalid python command"); + error->all(FLERR, "Unknown python command keyword: {}", arg[iarg]); } - if (pyfile && herestr) error->all(FLERR, "Invalid python command"); - if (pyfile && existflag) error->all(FLERR, "Invalid python command"); - if (herestr && existflag) error->all(FLERR, "Invalid python command"); + if (pyfile && herestr) + error->all(FLERR, "Must not use python 'file' and 'here' keywords at the same time"); + if (pyfile && existflag) + error->all(FLERR, "Must not use python 'file' and 'exists' keywords at the same time"); + if (herestr && existflag) + error->all(FLERR, "Must not use python 'here' and 'exists' keywords at the same time"); // create or overwrite entry in pfuncs vector with name = arg[0] @@ -219,23 +229,21 @@ void PythonImpl::command(int narg, char **arg) if (fp == nullptr) { PyUtils::Print_Errors(); - error->all(FLERR, "Could not open Python file"); + error->all(FLERR, "Could not open Python file: {}", pyfile); } int err = PyRun_SimpleFile(fp, pyfile); - if (err) { PyUtils::Print_Errors(); - error->all(FLERR, "Could not process Python file"); + error->all(FLERR, "Could not process Python file: {}", pyfile); } - fclose(fp); + } else if (herestr) { int err = PyRun_SimpleString(herestr); - if (err) { PyUtils::Print_Errors(); - error->all(FLERR, "Could not process Python string"); + error->all(FLERR, "Could not process Python string: {}", herestr); } } @@ -278,14 +286,17 @@ void PythonImpl::invoke_function(int ifunc, char *result) int ninput = pfuncs[ifunc].ninput; PyObject *pArgs = PyTuple_New(ninput); - if (!pArgs) { error->all(FLERR, "Could not create Python function arguments"); } + if (!pArgs) + error->all(FLERR, "Could not prepare arguments for Python function {}", pfuncs[ifunc].name); for (int i = 0; i < ninput; i++) { int itype = pfuncs[ifunc].itype[i]; if (itype == INT) { if (pfuncs[ifunc].ivarflag[i]) { str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); - if (!str) { error->all(FLERR, "Could not evaluate Python function input variable"); } + if (!str) + error->all(FLERR, "Could not evaluate Python function {} input variable: {}", + pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); pValue = PY_INT_FROM_LONG(atoi(str)); } else { pValue = PY_INT_FROM_LONG(pfuncs[ifunc].ivalue[i]); @@ -293,7 +304,9 @@ void PythonImpl::invoke_function(int ifunc, char *result) } else if (itype == DOUBLE) { if (pfuncs[ifunc].ivarflag[i]) { str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); - if (!str) { error->all(FLERR, "Could not evaluate Python function input variable"); } + if (!str) + error->all(FLERR, "Could not evaluate Python function {} input variable: {}", + pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); pValue = PyFloat_FromDouble(atof(str)); } else { pValue = PyFloat_FromDouble(pfuncs[ifunc].dvalue[i]); @@ -301,7 +314,9 @@ void PythonImpl::invoke_function(int ifunc, char *result) } else if (itype == STRING) { if (pfuncs[ifunc].ivarflag[i]) { str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); - if (!str) { error->all(FLERR, "Could not evaluate Python function input variable"); } + if (!str) + error->all(FLERR, "Could not evaluate Python function {} input variable: {}", + pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); pValue = PY_STRING_FROM_STRING(str); } else { pValue = PY_STRING_FROM_STRING(pfuncs[ifunc].svalue[i]); @@ -309,7 +324,7 @@ void PythonImpl::invoke_function(int ifunc, char *result) } else if (itype == PTR) { pValue = PY_VOID_POINTER(lmp); } else { - error->all(FLERR, "Unsupported variable type"); + error->all(FLERR, "Unsupported variable type: {}", itype); } PyTuple_SetItem(pArgs, i, pValue); } @@ -322,7 +337,7 @@ void PythonImpl::invoke_function(int ifunc, char *result) if (!pValue) { PyUtils::Print_Errors(); - error->one(FLERR, "Python function evaluation failed"); + error->one(FLERR, "Python evaluation of function {} failed", pfuncs[ifunc].name); } // function returned a value @@ -398,9 +413,10 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, int length_lon pfuncs[ifunc].noutput = noutput; if (!format && ninput + noutput) - error->all(FLERR, "Invalid python command"); + error->all(FLERR, "Missing python format keyword"); else if (format && ((int) strlen(format) != ninput + noutput)) - error->all(FLERR, "Invalid python command"); + error->all(FLERR, "Input/output arguments ({}) and format characters ({}) are inconsistent", + (ninput + noutput), strlen(format)); // process inputs as values or variables @@ -446,7 +462,7 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, int length_lon if (strcmp(istr[i], "SELF") != 0) error->all(FLERR, "Invalid python command"); } else - error->all(FLERR, "Invalid python command"); + error->all(FLERR, "Invalid python format character: {}", type); } // process output as value or variable @@ -463,7 +479,7 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, int length_lon else if (type == 's') pfuncs[ifunc].otype = STRING; else - error->all(FLERR, "Invalid python command"); + error->all(FLERR, "Invalid python return format character: {}", type); if (length_longstr) { if (pfuncs[ifunc].otype != STRING) @@ -484,7 +500,9 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, int length_lon int PythonImpl::execute_string(char *cmd) { PyUtils::GIL lock; - return PyRun_SimpleString(cmd); + int err = PyRun_SimpleString(cmd); + if (err) PyUtils::Print_Errors(); + return err; } /* ---------------------------------------------------------------------- */ @@ -496,6 +514,7 @@ int PythonImpl::execute_file(char *fname) PyUtils::GIL lock; int err = PyRun_SimpleFile(fp, fname); + if (err) PyUtils::Print_Errors(); if (fp) fclose(fp); return err; From 0890bc026e9f79650d3c05fccb426b0cb69b45d0 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 11 Nov 2022 15:28:36 -0500 Subject: [PATCH 48/58] 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 49/58] 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 50/58] 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 5e832aa36021451398e09e56543e382dfceefd2b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 11 Nov 2022 22:05:09 -0500 Subject: [PATCH 51/58] revise documentation --- doc/src/python.rst | 224 ++++++++++++++++++++++++++------------------- 1 file changed, 130 insertions(+), 94 deletions(-) diff --git a/doc/src/python.rst b/doc/src/python.rst index d777719181..1a79f679b5 100644 --- a/doc/src/python.rst +++ b/doc/src/python.rst @@ -8,20 +8,21 @@ Syntax .. parsed-literal:: - python function keyword args ... + python mode keyword args ... -* function = *source* or name of Python function +* mode = *source* or name of Python function - if function is *source*: + if mode is *source*: .. parsed-literal:: - keyword = *inline* or name of a Python file - inline = one or more lines of Python code which will be executed immediately - must be a single argument, typically enclosed in triple quotes - Python file = name of a file with Python code which will be executed immediately + keyword = *here* or name of a *Python file* + *here* arg = inline + inline = one or more lines of Python code which defines func + must be a single argument, typically enclosed between triple quotes + *Python file* = name of a file with Python code which will be executed immediately -* if function is the name of a Python function, one or more keyword/args pairs must be appended +* if *mode* is the name of a Python function, one or more keywords with/without arguments must be appended .. parsed-literal:: @@ -34,7 +35,7 @@ Syntax SELF = reference to LAMMPS itself which can be accessed by Python function variable = v_name, where name = name of LAMMPS variable, e.g. v_abc *return* arg = varReturn - varReturn = v_name = LAMMPS variable name which return value of function will be assigned to + varReturn = v_name = LAMMPS variable name which the return value of the Python function will be assigned to *format* arg = fstring with M characters M = N if no return value, where N = # of inputs M = N+1 if there is a return value @@ -80,55 +81,65 @@ Examples python source funcdef.py - python source inline "from lammps import lammps" + python source here "from lammps import lammps" Description """"""""""" -Define a Python function or execute a previously defined function or -execute some arbitrary python code. Arguments, including LAMMPS -variables, can be passed to the function from the LAMMPS input script -and a value returned by the Python function to a LAMMPS variable. The -Python code for the function can be included directly in the input -script or in a separate Python file. The function can be standard -Python code or it can make "callbacks" to LAMMPS through its library -interface to query or set internal values within LAMMPS. This is a -powerful mechanism for performing complex operations in a LAMMPS input -script that are not possible with the simple input script and variable -syntax which LAMMPS defines. Thus your input script can operate more -like a true programming language. +The *python* command allows to interface LAMMPS with an embedded Python +interpreter and enabled to either execute arbitrary python code in that +interpreter, register a Python function for future execution (as a python +style variable, from a fix interfaced with python or for direct +invocation), or invoke such a previously registered function. + +Arguments, including LAMMPS variables, can be passed to the function +from the LAMMPS input script and a value returned by the Python function +to a LAMMPS variable. The Python code for the function can be included +directly in the input script or in a separate Python file. The function +can be standard Python code or it can make "callbacks" to LAMMPS through +its library interface to query or set internal values within LAMMPS. +This is a powerful mechanism for performing complex operations in a +LAMMPS input script that are not possible with the simple input script +and variable syntax which LAMMPS defines. Thus your input script can +operate more like a true programming language. Use of this command requires building LAMMPS with the PYTHON package which links to the Python library so that the Python interpreter is embedded in LAMMPS. More details about this process are given below. -There are two ways to invoke a Python function once it has been defined. -One is using the *invoke* keyword. The other is to assign the function -to a :doc:`python-style variable ` defined in your input -script. Whenever the variable is evaluated, it will execute the Python -function to assign a value to the variable. Note that variables can be -evaluated in many different ways within LAMMPS. They can be substituted -with their result directly in an input script, or they can be passed to -various commands as arguments, so that the variable is evaluated during -a simulation run. +There are two ways to invoke a Python function once it has been +registered. One is using the *invoke* keyword. The other is to assign +the function to a :doc:`python-style variable ` defined in +your input script. Whenever the variable is evaluated, it will execute +the Python function to assign a value to the variable. Note that +variables can be evaluated in many different ways within LAMMPS. They +can be substituted with their result directly in an input script, or +they can be passed to various commands as arguments, so that the +variable is evaluated during a simulation run. A broader overview of how Python can be used with LAMMPS is given in the :doc:`Use Python with LAMMPS ` section of the -documentation. There is an ``examples/python`` directory which +documentation. There also is an ``examples/python`` directory which illustrates use of the python command. ---------- The first argument of the *python* command is either the *source* -keyword or the name of a Python function. +keyword or the name of a Python function. This defines the mode +of the python command. -If the *source* keyword is used, no other keywords can be used. The -argument either can be a filename or the keyword *inline* followed by a -string with python commands, either on a single line enclosed in quotes, -or as multiple lines enclosed in triple quotes. These python commands -will be passed to the python interpreter and executed immediately -without registering a python function for future execution. +If the *source* keyword is used, it is followed by either a file name or +the *here* keyword. No other keywords can be used. The *here* keyword +is followed by a string with python commands, either on a single line +enclosed in quotes, or as multiple lines enclosed in triple +quotes. These Python commands will be passed to the python interpreter +and executed immediately without registering a Python function for +future execution. This allows to run arbitrary Python code at any +time while processing the LAMMPS input file. This can be used to pre-load +Python modules, initialize global variables or perform operations +using the python programming language. The Python code will be executed +in parallel on all MPI processes. No arguments can be passed. In all other cases, the first argument is the name of a Python function that will be registered with LAMMPS for future execution. The function @@ -136,7 +147,7 @@ may already be defined (see *exists* keyword) or must be defined using the *file* or *here* keywords as explained below. If the *invoke* keyword is used, no other keywords can be used, and a -previous python command must have registered the Python function +previous *python* command must have registered the Python function referenced by this command. This invokes the Python function with the previously defined arguments and the return value is processed as explained below. You can invoke the function as many times as you wish @@ -144,21 +155,23 @@ in your input script. The *input* keyword defines how many arguments *N* the Python function expects. If it takes no arguments, then the *input* keyword should not -be used. Each argument can be specified directly as a value, e.g. 6 or -3.14159 or abc (a string of characters). The type of each argument is -specified by the *format* keyword as explained below, so that Python -will know how to interpret the value. If the word SELF is used for an -argument it has a special meaning. A pointer is passed to the Python -function which it converts into a reference to LAMMPS itself. This +be used. Each argument can be specified directly as a value, e.g. '6' +or '3.14159' or 'abc' (a string of characters). The type of each +argument is specified by the *format* keyword as explained below, so +that Python will know how to interpret the value. If the word SELF is +used for an argument it has a special meaning. A pointer is passed to +the Python function which it can convert into a reference to LAMMPS +itself using the :doc:`LAMMPS Python module `. This enables the function to call back to LAMMPS through its library interface as explained below. This allows the Python function to query or set values internal to LAMMPS which can affect the subsequent execution of the input script. A LAMMPS variable can also be used as an argument, specified as v_name, where "name" is the name of the variable. -Any style of LAMMPS variable can be used, as defined by the -:doc:`variable ` command. Each time the Python function is -invoked, the LAMMPS variable is evaluated and its value is passed to the -Python function. +Any style of LAMMPS variable returning a scalar or a string can be used, +as defined by the :doc:`variable ` command. The *format* +keyword must be used to set the type of data that is passed to Python. +Each time the Python function is invoked, the LAMMPS variable is +evaluated and its value is passed to the Python function. The *return* keyword is only needed if the Python function returns a value. The specified *varReturn* must be of the form v_name, where @@ -168,8 +181,9 @@ numeric or string value, as specified by the *format* keyword. As explained on the :doc:`variable ` doc page, the definition of a python-style variable associates a Python function name with the -variable. This must match the *func* setting for this command. For -example these two commands would be self-consistent: +variable. This must match the *Python function name* first argument of +the *python* command. For example these two commands would be +consistent: .. code-block:: LAMMPS @@ -178,10 +192,11 @@ example these two commands would be self-consistent: The two commands can appear in either order in the input script so long as both are specified before the Python function is invoked for -the first time. +the first time. Afterwards, the variable 'foo' is associated with +the Python function 'myMultiply'. -The *format* keyword must be used if the *input* or *return* keyword is -used. It defines an *fstring* with M characters, where M = sum of +The *format* keyword must be used if the *input* or *return* keywords +are used. It defines an *fstring* with M characters, where M = sum of number of inputs and outputs. The order of characters corresponds to the N inputs, followed by the return value (if it exists). Each character must be one of the following: "i" for integer, "f" for @@ -246,13 +261,14 @@ however contain triple-quoted Python strings, since that will conflict with the triple-quote parsing that the LAMMPS input script performs. All the Python code you specify via one or more python commands is -loaded into the Python "main" module, i.e. __main__. The code can -define global variables or statements that are outside of function -definitions. It can contain multiple functions, only one of which -matches the *func* setting in the python command. This means you can -use the *file* keyword once to load several functions, and the *exists* -keyword thereafter in subsequent python commands to access the other -functions previously loaded. +loaded into the Python "main" module, i.e. ``__name__ == '__main__'``. +The code can define global variables, define global functions, define +classes or execute statements that are outside of function definitions. +It can contain multiple functions, only one of which matches the *func* +setting in the python command. This means you can use the *file* +keyword once to load several functions, and the *exists* keyword +thereafter in subsequent python commands to register the other functions +that were previously loaded with LAMMPS. A Python function you define (or more generally, the code you load) can import other Python modules or classes, it can make calls to other @@ -279,12 +295,13 @@ outside the function: nvaluelast = nvalue return nvalue -Nsteplast stores the previous timestep the function was invoked -(passed as an argument to the function). Nvaluelast stores the return -value computed on the last function invocation. If the function is -invoked again on the same timestep, the previous value is simply -returned, without re-computing it. The "global" statement inside the -Python function allows it to overwrite the global variables. +The variable 'nsteplast' stores the previous timestep the function was +invoked (passed as an argument to the function). The variable +'nvaluelast' stores the return value computed on the last function +invocation. If the function is invoked again on the same timestep, the +previous value is simply returned, without re-computing it. The +"global" statement inside the Python function allows it to overwrite the +global variables from within the local context of the function. Note that if you load Python code multiple times (via multiple python commands), you can overwrite previously loaded variables and functions @@ -300,19 +317,39 @@ copy of the Python function(s) you define. There is no connection between the Python interpreters running on different processors. This implies three important things. -First, if you put a print statement in your Python function, you will -see P copies of the output, when running on P processors. If the -prints occur at (nearly) the same time, the P copies of the output may -be mixed together. Welcome to the world of parallel programming and -debugging. +First, if you put a print or other statement creating output to the +screen in your Python function, you will see P copies of the output, +when running on P processors. If the prints occur at (nearly) the same +time, the P copies of the output may be mixed together. When loading +the LAMMPS Python module into the embedded Python interpreter, it is +possible to pass the pointer to the current LAMMPS class instance and +via the Python interface to the LAMMPS library interface, it is possible +to determine the MPI rank of the current process and thus adapt the +Python code so that output will only appear on MPI rank 0. The +following LAMMPS input demonstrates how this could be done. The text +'Hello, LAMMPS!' should be printed only once, even when running LAMMPS +in parallel. -Second, if your Python code loads modules that are not pre-loaded by -the Python library, then it will load the module from disk. This may -be a bottleneck if 1000s of processors try to load a module at the -same time. On some large supercomputers, loading of modules from disk -by Python may be disabled. In this case you would need to pre-build a -Python library that has the required modules pre-loaded and link -LAMMPS with that library. +.. code-block:: LAMMPS + + python python_hello input 1 SELF format p here """ + def python_hello(handle): + from lammps import lammps + lmp = lammps(ptr=handle) + me = lmp.extract_setting('world_rank') + if me == 0: + print('Hello, LAMMPS!') + """ + + python python_hello invoke + +If your Python code loads Python modules that are not pre-loaded by the +Python library, then it will load the module from disk. This may be a +bottleneck if 1000s of processors try to load a module at the same time. +On some large supercomputers, loading of modules from disk by Python may +be disabled. In this case you would need to pre-build a Python library +that has the required modules pre-loaded and link LAMMPS with that +library. Third, if your Python code calls back to LAMMPS (discussed in the next section) and causes LAMMPS to perform an MPI operation requires @@ -330,22 +367,21 @@ Python function is as follows: .. code-block:: python - def foo(lmpptr,...): + def foo(handle,...): from lammps import lammps - lmp = lammps(ptr=lmpptr) + lmp = lammps(ptr=handle) lmp.command('print "Hello from inside Python"') ... -The function definition must include a variable (lmpptr in this case) -which corresponds to SELF in the python command. The first line of the -function imports the :doc:`"lammps" Python module `. -The second line creates a Python object ``lmp`` which -wraps the instance of LAMMPS that called the function. The "ptr=lmpptr" -argument is what makes that happen. The third line invokes the -command() function in the LAMMPS library interface. It takes a single -string argument which is a LAMMPS input script command for LAMMPS to -execute, the same as if it appeared in your input script. In this case, -LAMMPS should output +The function definition must include a variable ('handle' in this case) +which corresponds to SELF in the *python* command. The first line of +the function imports the :doc:`"lammps" Python module `. +The second line creates a Python object ``lmp`` which wraps the instance +of LAMMPS that called the function. The 'ptr=handle' argument is what +makes that happen. The third line invokes the command() function in the +LAMMPS library interface. It takes a single string argument which is a +LAMMPS input script command for LAMMPS to execute, the same as if it +appeared in your input script. In this case, LAMMPS should output .. parsed-literal:: @@ -359,8 +395,8 @@ The :doc:`Python_head` page describes the syntax for how Python wraps the various functions included in the LAMMPS library interface. -A more interesting example is in the examples/python/in.python script -which loads and runs the following function from examples/python/funcs.py: +A more interesting example is in the ``examples/python/in.python`` script +which loads and runs the following function from ``examples/python/funcs.py``: .. code-block:: python From d4bd0a74a7e013147c791e0b73da49b504d76b7d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 11 Nov 2022 22:05:35 -0500 Subject: [PATCH 52/58] change "inline" keyword to "here" for consistency with the other uses --- src/PYTHON/python_impl.cpp | 2 +- unittest/python/test_python_package.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index cacf1ca237..9177381a35 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -139,7 +139,7 @@ void PythonImpl::command(int narg, char **arg) if ((narg > 1) && (strcmp(arg[0], "source") == 0)) { int err = -1; - if ((narg > 2) && (strcmp(arg[1], "inline") == 0)) { + if ((narg > 2) && (strcmp(arg[1], "here") == 0)) { err = execute_string(arg[2]); } else { if (platform::file_is_readable(arg[1])) diff --git a/unittest/python/test_python_package.cpp b/unittest/python/test_python_package.cpp index 8b4f998bfc..f998ebbcd8 100644 --- a/unittest/python/test_python_package.cpp +++ b/unittest/python/test_python_package.cpp @@ -286,7 +286,7 @@ TEST_F(PythonPackageTest, RunSourceInline) { // execute inline python script auto output = CAPTURE_OUTPUT([&] { - command("python source inline \"\"\"\n" + command("python source here \"\"\"\n" "from __future__ import print_function\n" "print(2+2)\n" "\"\"\""); From cc4c59649cdd510ce4a8df69324df7cf34606583 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 11 Nov 2022 22:06:53 -0500 Subject: [PATCH 53/58] spelling --- doc/utils/sphinx-config/false_positives.txt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 61c1d45ba7..13a3682870 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2248,6 +2248,7 @@ MxN myCompute myIndex mylammps +myMultiply MyPool mysocket mySpin @@ -2473,7 +2474,7 @@ nsq Nstart nstats Nstep -Nsteplast +nsteplast Nstop nsub Nsw @@ -2503,7 +2504,7 @@ numpy Numpy Nurdin Nvalue -Nvaluelast +nvaluelast Nvalues nvc nvcc From eee6862f2529ff63edfce30dc0642811c3496e77 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 12 Nov 2022 02:13:39 -0500 Subject: [PATCH 54/58] 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; From 175b2b045a930b8507ce6cff8c0667f9ee1a1d5a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Nov 2022 14:10:43 -0500 Subject: [PATCH 55/58] tweak grammar --- doc/src/python.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/python.rst b/doc/src/python.rst index 1a79f679b5..0efa27555b 100644 --- a/doc/src/python.rst +++ b/doc/src/python.rst @@ -88,10 +88,10 @@ Description """"""""""" The *python* command allows to interface LAMMPS with an embedded Python -interpreter and enabled to either execute arbitrary python code in that -interpreter, register a Python function for future execution (as a python -style variable, from a fix interfaced with python or for direct -invocation), or invoke such a previously registered function. +interpreter and enables either executing arbitrary python code in that +interpreter, registering a Python function for future execution (as a +python style variable, from a fix interfaced with python or for direct +invocation), or invoking such a previously registered function. Arguments, including LAMMPS variables, can be passed to the function from the LAMMPS input script and a value returned by the Python function From c191086812670186f3afc7c980c1b744359ef2a0 Mon Sep 17 00:00:00 2001 From: ssande7 <1731652+ssande7@users.noreply.github.com> Date: Tue, 15 Nov 2022 14:52:24 +1000 Subject: [PATCH 56/58] Fix segfault when using dynamic groups with r-RESPA --- src/respa.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/respa.cpp b/src/respa.cpp index 650b04e6d6..530c35b3fb 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -372,7 +372,7 @@ void Respa::setup(int flag) mesg += fmt::format(" {}:{}", ilevel + 1, step[ilevel]); mesg += "\n r-RESPA fixes :"; - for (int l = 0; l < modify->n_post_force_respa_any; ++l) { + for (int l = 0; l < modify->n_post_force_respa; ++l) { Fix *f = modify->get_fix_by_index(modify->list_post_force_respa[l]); if (f->respa_level >= 0) mesg += fmt::format(" {}:{}[{}]", MIN(f->respa_level + 1, nlevels), f->style, f->id); From 742265bfdbae7e513e9607e6bf843d22951ac0c0 Mon Sep 17 00:00:00 2001 From: davidfir3 <491197586@qq.com> Date: Tue, 15 Nov 2022 16:39:38 +0800 Subject: [PATCH 57/58] append suffix to pair style --- src/FEP/compute_fep.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/FEP/compute_fep.cpp b/src/FEP/compute_fep.cpp index bb11a64f54..13fcb98850 100644 --- a/src/FEP/compute_fep.cpp +++ b/src/FEP/compute_fep.cpp @@ -198,8 +198,12 @@ void ComputeFEP::init() if (pert->which == PAIR) { pairflag = 1; + Pair *pair = nullptr; - Pair *pair = force->pair_match(pert->pstyle, 1); + if (lmp->suffix_enable) + pair = force->pair_match(std::string(pert->pstyle)+"/"+lmp->suffix,1); + + if (pair == nullptr) pair = force->pair_match(pert->pstyle,1); if (pair == nullptr) error->all(FLERR, "compute fep pair style " From 3fd122311d1bfca88337de5bfdb8700d0fe0dc17 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 15 Nov 2022 03:45:04 -0500 Subject: [PATCH 58/58] some more grammar updates and clarifications --- doc/src/python.rst | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/doc/src/python.rst b/doc/src/python.rst index 0efa27555b..e5854faf46 100644 --- a/doc/src/python.rst +++ b/doc/src/python.rst @@ -87,15 +87,15 @@ Examples Description """"""""""" -The *python* command allows to interface LAMMPS with an embedded Python +The *python* command allows interfacing LAMMPS with an embedded Python interpreter and enables either executing arbitrary python code in that interpreter, registering a Python function for future execution (as a -python style variable, from a fix interfaced with python or for direct +python style variable, from a fix interfaced with python, or for direct invocation), or invoking such a previously registered function. Arguments, including LAMMPS variables, can be passed to the function from the LAMMPS input script and a value returned by the Python function -to a LAMMPS variable. The Python code for the function can be included +assigned to a LAMMPS variable. The Python code for the function can be included directly in the input script or in a separate Python file. The function can be standard Python code or it can make "callbacks" to LAMMPS through its library interface to query or set internal values within LAMMPS. @@ -132,14 +132,16 @@ of the python command. If the *source* keyword is used, it is followed by either a file name or the *here* keyword. No other keywords can be used. The *here* keyword is followed by a string with python commands, either on a single line -enclosed in quotes, or as multiple lines enclosed in triple -quotes. These Python commands will be passed to the python interpreter -and executed immediately without registering a Python function for -future execution. This allows to run arbitrary Python code at any -time while processing the LAMMPS input file. This can be used to pre-load -Python modules, initialize global variables or perform operations -using the python programming language. The Python code will be executed -in parallel on all MPI processes. No arguments can be passed. +enclosed in quotes, or as multiple lines enclosed in triple quotes. +These Python commands will be passed to the python interpreter and +executed immediately without registering a Python function for future +execution. The code will be loaded into and run in the "main" module of +the Python interpreter. This allows running arbitrary Python code at +any time while processing the LAMMPS input file. This can be used to +pre-load Python modules, initialize global variables, define functions +or classes, or perform operations using the python programming language. +The Python code will be executed in parallel on all MPI processes. No +arguments can be passed. In all other cases, the first argument is the name of a Python function that will be registered with LAMMPS for future execution. The function @@ -225,7 +227,8 @@ one of them. These keywords specify what Python code to load into the Python interpreter. The *file* keyword gives the name of a file containing Python code, which should end with a ".py" suffix. The code will be immediately loaded into and run in the "main" module of the -Python interpreter. Note that Python code which contains a function +Python interpreter. The Python code will be executed in parallel on all +MPI processes. Note that Python code which contains a function definition does not "execute" the function when it is run; it simply defines the function so that it can be invoked later.