diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index c78c05a35e..383b8212f9 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -39,6 +39,9 @@ Syntax *masslimit* value = massmin massmax massmin = minimum molecular weight of species to delete massmax = maximum molecular weight of species to delete + *delete_rate_limit* value = Nlimit Nsteps + Nlimit = maximum number of deletions allowed to occur within interval + Nsteps = the interval (number of timesteps) over which to count deletions Examples """""""" @@ -142,7 +145,13 @@ When using the *masslimit* keyword, each line of the *filedel* file contains the timestep on which deletions occurs, followed by how many of each species are deleted (with quantities preceding chemical formulae). The *specieslist* and *masslimit* keywords cannot both be -used in the same *reaxff/species* fix. +used in the same *reaxff/species* fix. The *delete_rate_limit* +keyword can enforce an upper limit on the overall rate of molecule +deletion. The number of deletion occurrences is limited to Nlimit +within an interval of Nsteps timesteps. When using the +*delete_rate_limit* keyword, no deletions are permitted to occur +within the first Nsteps timesteps of the first run (after reading a +either a data or restart file). ---------- diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 65a2e6d8ce..ce04be2cc8 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -21,6 +21,7 @@ #include "atom.h" #include "atom_vec.h" +#include "citeme.h" #include "comm.h" #include "domain.h" #include "error.h" @@ -36,12 +37,25 @@ #include "pair_reaxff.h" #include "reaxff_defs.h" +#include #include #include +#include using namespace LAMMPS_NS; using namespace FixConst; +static const char cite_reaxff_species_delete[] = + "fix reaxff/species, 'delete' keyword: https://doi.org/10.1016/j.carbon.2022.11.002\n\n" + "@Article{Gissinger23,\n" + " author = {J. R. Gissinger, S. R. Zavada, J. G. Smith, J. Kemppainen, I. Gallegos, G. M. Odegard, E. J. Siochi, K. E. Wise},\n" + " title = {Predicting char yield of high-temperature resins},\n" + " journal = {Carbon},\n" + " year = 2023,\n" + " volume = 202,\n" + " pages = {336-347}\n" + "}\n\n"; + /* ---------------------------------------------------------------------- */ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : @@ -145,6 +159,7 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : ele = filepos = filedel = nullptr; eleflag = posflag = padflag = 0; delflag = specieslistflag = masslimitflag = 0; + delete_Nlimit = delete_Nsteps = 0; singlepos_opened = multipos_opened = del_opened = 0; multipos = 0; @@ -221,7 +236,12 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR, "Unknown fix reaxff/species delete option: {}", arg[iarg]); - + // rate limit when deleting molecules + } else if (strcmp(arg[iarg], "delete_rate_limit") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix reaxff/species delete_rate_limit", error); + delete_Nlimit = utils::numeric(FLERR, arg[iarg+1], false, lmp); + delete_Nsteps = utils::numeric(FLERR, arg[iarg+2], false, lmp); + iarg += 3; // position of molecules } else if (strcmp(arg[iarg], "position") == 0) { if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix reaxff/species position", error); @@ -260,6 +280,15 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : if (delflag && specieslistflag && masslimitflag) error->all(FLERR, "Incompatible combination fix reaxff/species command options"); + if (delete_Nlimit > 0) { + if (lmp->citeme) lmp->citeme->add(cite_reaxff_species_delete); + memory->create(delete_Tcount,delete_Nsteps,"reaxff/species:delete_Tcount"); + + for (int i = 0; i < delete_Nsteps; i++) + delete_Tcount[i] = -1; + delete_Tcount[0] = 0; + } + vector_nmole = 0; vector_nspec = 0; } @@ -279,6 +308,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies() memory->destroy(Mol2Spec); memory->destroy(MolType); memory->destroy(MolName); + memory->destroy(delete_Tcount); delete[] filepos; delete[] filedel; @@ -375,7 +405,13 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) // point to fix_ave_atom f_SPECBOND->end_of_step(); - if (ntimestep != nvalid) return; + if (ntimestep != nvalid) { + // push back delete_Tcount on every step + if (delete_Nlimit > 0) + for (int i = delete_Nsteps-1; i > 0; i--) + delete_Tcount[i] = delete_Tcount[i-1]; + return; + } nlocal = atom->nlocal; @@ -826,6 +862,15 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) { + int ndeletions; + int headroom = -1; + if (delete_Nlimit > 0) { + if (delete_Tcount[delete_Nsteps-1] == -1) return; + ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps-1]; + headroom = MAX(0, delete_Nlimit - ndeletions); + if (headroom == 0) return; + } + int i, j, m, n, itype, cid; int ndel, ndelone, count, count_tmp; int *Nameall; @@ -856,7 +901,23 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) int *marklist; memory->create(marklist, nlocal, "reaxff/species:marklist"); - for (m = 1; m <= Nmole; m++) { + std::random_device rnd; + std::minstd_rand park_rng(rnd()); + int *molrange; + memory->create(molrange,Nmole,"reaxff/species:molrange"); + for (m = 0; m < Nmole; m++) + molrange[m] = m + 1; + if (delete_Nlimit > 0) { + // shuffle index when using rate_limit, in case order is biased + if (comm->me == 0) + std::shuffle(&molrange[0],&molrange[Nmole], park_rng); + MPI_Bcast(&molrange[0], Nmole, MPI_INT, 0, world); + } + + int this_delete_Tcount = 0; + for (int mm = 0; mm < Nmole; mm++) { + if (this_delete_Tcount == headroom) break; + m = molrange[mm]; localmass = totalmass = count = nmarklist = 0; for (n = 0; n < ntypes; n++) Name[n] = 0; @@ -896,6 +957,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) // find corresponding moltype if (totalmass > massmin && totalmass < massmax) { + this_delete_Tcount++; for (j = 0; j < nmarklist; j++) { mark[marklist[j]] = 1; deletecount[Mol2Spec[m - 1]] += 1.0 / (double) count; @@ -905,6 +967,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) if (count > 0) { for (i = 0; i < ndelspec; i++) { if (del_species[i] == species_str) { + this_delete_Tcount++; for (j = 0; j < nmarklist; j++) { mark[marklist[j]] = 1; deletecount[i] += 1.0 / (double) count; @@ -976,6 +1039,14 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) } } + + // push back delete_Tcount on every step + if (delete_Nlimit > 0) { + for (i = delete_Nsteps-1; i > 0; i--) + delete_Tcount[i] = delete_Tcount[i-1]; + delete_Tcount[0] += this_delete_Tcount; + } + if (ndel && (atom->map_style != Atom::MAP_NONE)) { atom->nghost = 0; atom->map_init(); @@ -988,6 +1059,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) memory->destroy(marklist); memory->destroy(mark); memory->destroy(deletecount); + memory->destroy(molrange); } /* ---------------------------------------------------------------------- */ diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index 65eeae4c60..329e17145b 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -60,6 +60,7 @@ class FixReaxFFSpecies : public Fix { FILE *fp, *pos, *fdel; int eleflag, posflag, multipos, padflag, setupflag; int delflag, specieslistflag, masslimitflag; + int delete_Nlimit, delete_Nsteps, *delete_Tcount; double massmin, massmax; int singlepos_opened, multipos_opened, del_opened; char *ele, **eletype, *filepos, *filedel;