lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3514 from jrgissing/per-bond_custom_constraint 

bond/react feature updates
This commit is contained in:
Axel Kohlmeyer 2022-11-15 21:03:12 -05:00 committed by GitHub
parent 7a42d2a54a eee6862f25
commit 4655039774
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 384 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
doc/src/fix_bond_react.rst
View File

@ -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 *rescale_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
@ -56,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
@ -514,28 +520,40 @@ 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::
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 <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 <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 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 +565,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) <variable>`
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 l0 equal 1.0 # characteristic length
.. code-block:: LAMMPS
# in Constraints section of map file
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,
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
@ -598,6 +637,15 @@ 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
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 <variable>`.
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
@ -616,6 +664,19 @@ 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
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

351
src/REACTION/fix_bond_react.cpp
View File

@ -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};
@ -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;
@ -224,8 +229,10 @@ 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(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");
@ -249,8 +256,11 @@ 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;
rescale_charges_flag[i] = 0;
create_atoms_flag[i] = 0;
modify_create_fragid[i] = -1;
overlapsq[i] = 0;
@ -286,55 +296,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++]);
@ -344,7 +330,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++;
@ -354,14 +340,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);
@ -380,6 +360,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");
@ -398,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");
@ -433,21 +428,24 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
}
max_natoms = 0; // the number of atoms in largest molecule template
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
@ -462,6 +460,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++) {
@ -471,7 +473,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");
@ -598,6 +600,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);
@ -617,8 +620,10 @@ 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(rescale_charges_flag);
memory->destroy(molecule_keyword);
memory->destroy(nconstraints);
memory->destroy(constraintstr);
@ -819,6 +824,16 @@ 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++) {
@ -873,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;
@ -911,8 +928,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;
@ -1385,9 +1416,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;
@ -1395,8 +1442,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++)
@ -1414,14 +1464,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;
@ -2118,7 +2169,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 +2191,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 +2313,67 @@ 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<tagint> aset;
std::string computeid = varid;
std::map<std::set<tagint>,int>::iterator it;
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");
if (atoms2bondflag == 0) {
atoms2bondflag = 1;
get_atoms2bond(cperbond->groupbit);
}
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");
if (cperbond->invoked_local != lmp->update->ntimestep)
cperbond->compute_local();
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];
return perbondval;
}
int ivar = -1;
for (int i = 0; i < ncustomvars; i++) {
if (varid == customvarstrs[i]) {
@ -2280,13 +2387,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;
@ -2313,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<tagint> 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
------------------------------------------------------------------------- */
@ -2985,6 +3118,8 @@ void FixBondReact::update_everything()
}
delete [] iskip;
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
// invoke map_init() b/c atom count has grown
@ -3012,6 +3147,33 @@ void FixBondReact::update_everything()
}
}
// get charge rescale delta
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];
@ -3024,7 +3186,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;
}
}
}
@ -3795,10 +3957,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;
@ -4299,34 +4475,71 @@ 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;
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';
}
if (me == 0) {
int size = nreacts*sizeof(Set);
fwrite(&size,sizeof(int),1,fp);
fwrite(set,sizeof(Set),nreacts,fp);
int rbufcount = max_rate_limit_steps*nreacts;
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+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);
}
if (rbufcount) memory->destroy(rbuf);
}
/* ----------------------------------------------------------------------
use selected state info from restart file to restart the Fix
bond/react restart revisions numbers added after LAMMPS version 3 Nov 2022
------------------------------------------------------------------------- */
void FixBondReact::restart(char *buf)
{
Set *set_restart = (Set *) buf;
for (int i = 0; i < set_restart[0].nreacts; i++) {
int n,revision,r_nreacts,r_max_rate_limit_steps,ibufcount,n2cpy;
int **ibuf;
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[sizeof(int)+r_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++) {
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];
}
}
}
if (revision > 0) memory->destroy(ibuf);
}
/* ----------------------------------------------------------------------

26
src/REACTION/fix_bond_react.h
View File

@ -25,6 +25,10 @@ FixStyle(bond/react,FixBondReact);
#define LMP_FIX_BOND_REACT_H
#include "fix.h"
#include "compute.h"
#include <map>
#include <set>
namespace LAMMPS_NS {
@ -64,8 +68,11 @@ 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 *rescale_charges_flag;
int *create_atoms_flag;
int *modify_create_fragid;
double *overlapsq;
@ -74,9 +81,11 @@ class FixBondReact : public Fix {
int *nconstraints;
char **constraintstr;
int nrxnfunction;
std::vector<std::string> rxnfunclist;
std::vector<std::string> rxnfunclist; // lists current special rxn function
std::vector<int> 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;
@ -86,6 +95,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;
@ -160,7 +170,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);
@ -184,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 *);
@ -198,16 +210,18 @@ 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;
void restart(char *buf) override;
// store restart data
struct Set {
int nreacts;
char rxn_name[MAXLINE];
int reaction_count_total;
int max_rate_limit_steps;
};
Set *set;
@ -221,7 +235,9 @@ class FixBondReact : public Fix {
int ncustomvars;
std::vector<std::string> 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::map<std::set<tagint>, int> atoms2bond; // maps atom pair to index of local bond array
std::vector<std::vector<Constraint>> constraints;
// DEBUG