diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index c557d2617d..15098761cc 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 *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 ` 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 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) ` -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 `. + 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 diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp index 5dc5967c29..8857590c4f 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}; @@ -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 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_'"); + 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 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); } /* ---------------------------------------------------------------------- diff --git a/src/REACTION/fix_bond_react.h b/src/REACTION/fix_bond_react.h index e9ca7ec362..b55608b07e 100644 --- a/src/REACTION/fix_bond_react.h +++ b/src/REACTION/fix_bond_react.h @@ -25,6 +25,10 @@ FixStyle(bond/react,FixBondReact); #define LMP_FIX_BOND_REACT_H #include "fix.h" +#include "compute.h" + +#include +#include 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 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 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 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, int> atoms2bond; // maps atom pair to index of local bond array std::vector> constraints; // DEBUG