diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index a652777ba7..11f0e7b7e7 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -20,7 +20,7 @@ Syntax * Nfreq = calculate average bond-order every this many timesteps * filename = name of output file * zero or more keyword/value pairs may be appended -* keyword = *cutoff* or *element* or *position* +* keyword = *cutoff* or *element* or *position* or *delete* .. parsed-literal:: @@ -31,6 +31,14 @@ Syntax *position* value = posfreq filepos posfreq = write position files every this many timestep filepos = name of position output file + *delete* value = filedel keyword value + filedel = name of delete species output file + keyword = *specieslist* or *masslimit* + *specieslist* value = Nspecies Species1 Species2 ... + Nspecies = number of species in list + *masslimit* value = massmin massmax + massmin = minimum molecular weight of species to delete + massmax = maximum molecular weight of species to delete Examples """""""" @@ -40,6 +48,7 @@ Examples fix 1 all reaxff/species 10 10 100 species.out fix 1 all reaxff/species 1 2 20 species.out cutoff 1 1 0.40 cutoff 1 2 0.55 fix 1 all reaxff/species 1 100 100 species.out element Au O H position 1000 AuOH.pos + fix 1 all reaxff/species 1 100 100 species.out delete species.del masslimit 0 50 Description """"""""""" @@ -59,13 +68,18 @@ the first line. .. warning:: In order to compute averaged data, it is required that there are no - neighbor list rebuilds between the *Nfreq* steps. For that reason, fix - *reaxff/species* may change your neighbor list settings. There will - be a warning message showing the new settings. Having an *Nfreq* - setting that is larger than what is required for correct computation - of the ReaxFF force field interactions can thus lead to incorrect - results. For typical ReaxFF calculations a value of 100 is already - quite large. + neighbor list rebuilds for at least Nrepeat\*Nevery steps preceding + each *Nfreq* step. For that reason, fix *reaxff/species* may + change your neighbor list settings. Reneighboring will occur no + more frequently than every Nrepeat\*Nevery timesteps, and will + occur less frequently if *Nfreq* is not a multiple of + Nrepeat\*Nevery. There will be a warning message showing the new + settings. Having a *Nfreq* setting that is larger than what is + required for correct computation of the ReaxFF force field + interactions, in combination with certain *Nrepeat* and *Nevery* + settings, can thus lead to incorrect results. For typical ReaxFF + calculations, reneighboring only every 100 steps is already quite a + low frequency. If the filename ends with ".gz", the output file is written in gzipped format. A gzipped dump file will be about 3x smaller than the text version, @@ -104,6 +118,30 @@ character appears in *filepos*, then one file per snapshot is written at *posfreq* and the "\*" character is replaced with the timestep value. For example, AuO.pos.\* becomes AuO.pos.0, AuO.pos.1000, etc. +The optional keyword *delete* enables the periodic removal of +molecules from the system. Criteria for deletion can be either a list +of specific chemical formulae or a range of molecular weights. +Molecules are deleted every *Nfreq* timesteps, and bond connectivity +is determined using the *Nevery* and *Nrepeat* keywords. The +*filedel* argument is the name of the output file that records the +species that are removed from the system. The *specieslist* keyword +permits specific chemical species to be deleted. The *Nspecies* +argument specifies how many species are eligible for deletion and is +followed by a list of chemical formulae, whose strings are compared to +species identified by this fix. For example, "specieslist 2 CO CO2" +deletes molecules that are identified as "CO" and "CO2" in the species +output file. When using the *specieslist* keyword, the *filedel* file +has the following format: the first line lists the chemical formulae +eligible for deletion, and each additional line contains the timestep +on which a molecule deletion occurs and the number of each species +deleted on that timestep. The *masslimit* keyword permits deletion of +molecules with molecular weights between *massmin* and *massmax*. +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. + ---------- The *Nevery*, *Nrepeat*, and *Nfreq* arguments specify on what diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index a89cd6f0dc..e5298bdc4d 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -14,11 +14,13 @@ /* ---------------------------------------------------------------------- Contributing authors: Ray Shan (Sandia, tnshan@sandia.gov) Oleg Sergeev (VNIIA, sergeev@vniia.ru) + Jacob Gissinger (NASA, jacob.r.gissinger@gmail.com), 'delete' keyword ------------------------------------------------------------------------- */ #include "fix_reaxff_species.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" #include "domain.h" #include "error.h" @@ -45,7 +47,8 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, { if (narg < 7) error->all(FLERR, "Illegal fix reaxff/species command"); - force_reneighbor = 0; + force_reneighbor = 1; + next_reneighbor = -1; vector_flag = 1; size_vector = 2; @@ -125,8 +128,10 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, MolName = nullptr; MolType = nullptr; NMol = nullptr; + Mol2Spec = nullptr; nd = nullptr; molmap = nullptr; + mark = nullptr; nmax = 0; setupflag = 0; @@ -142,10 +147,11 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, // optional args eletype = nullptr; - ele = filepos = nullptr; + ele = filepos = filedel = nullptr; eleflag = posflag = padflag = 0; + delflag = specieslistflag = masslimitflag = 0; - singlepos_opened = multipos_opened = 0; + singlepos_opened = multipos_opened = del_opened = 0; multipos = 0; posfreq = 0; @@ -180,6 +186,44 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, eleflag = 1; iarg += ntypes + 1; + // delete species + } else if (strcmp(arg[iarg],"delete") == 0) { + delflag = 1; + filedel = new char[255]; + strcpy(filedel,arg[iarg+1]); + if (me == 0) { + fdel = fopen(filedel, "w"); + if (fdel == nullptr) error->one(FLERR,"Cannot open fix reaxff/species delete file"); + } + + del_opened = 1; + + if (strcmp(arg[iarg+2],"masslimit") == 0) { + if (iarg+5 > narg) error->all(FLERR,"Illegal fix reaxff/species command"); + masslimitflag = 1; + massmin = atof(arg[iarg+3]); + massmax = atof(arg[iarg+4]); + iarg += 5; + } else if (strcmp(arg[iarg+2],"specieslist") == 0) { + specieslistflag = 1; + ndelspec = atoi(arg[iarg+3]); + if (iarg+ndelspec+4 > narg) error->all(FLERR,"Illegal fix reaxff/species command"); + + del_species.resize(ndelspec); + for (int i = 0; i < ndelspec; i ++) + del_species[i] = arg[iarg+4+i]; + + if (me == 0) { + fprintf(fdel,"Timestep"); + for (i = 0; i < ndelspec; i++) + fprintf(fdel,"\t%s",del_species[i].c_str()); + fprintf(fdel,"\n"); + fflush(fdel); + } + + iarg += ndelspec + 4; + } else error->all(FLERR, "Illegal fix reaxff/species command"); + // position of molecules } else if (strcmp(arg[iarg], "position") == 0) { if (iarg + 3 > narg) error->all(FLERR, "Illegal fix reaxff/species command"); @@ -213,6 +257,9 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, if (ntypes > 3) ele[3] = 'N'; } + if (delflag && specieslistflag && masslimitflag) + error->all(FLERR, "Illegal fix reaxff/species command"); + vector_nmole = 0; vector_nspec = 0; } @@ -229,6 +276,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies() memory->destroy(nd); memory->destroy(Name); memory->destroy(NMol); + memory->destroy(Mol2Spec); memory->destroy(MolType); memory->destroy(MolName); @@ -358,6 +406,8 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) if (me == 0) fflush(pos); } + if (delflag) DeleteSpecies(Nmole, Nspec); + nvalid += nfreq; } @@ -540,6 +590,11 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) memory->create(Nameall, ntypes, "reaxff/species:Nameall"); memory->create(NMolall, Nmole, "reaxff/species:NMolall"); + memory->destroy(Mol2Spec); + Mol2Spec = nullptr; + memory->create(Mol2Spec, Nmole, "reaxff/species:Mol2Spec"); + for (m = 0; m < Nmole; m++) Mol2Spec[m] = -1; + for (m = 1, Nspec = 0; m <= Nmole; m++) { for (n = 0; n < ntypes; n++) Name[n] = 0; for (n = 0, flag_mol = 0; n < nlocal; n++) { @@ -563,11 +618,15 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) flag_spec = 0; for (l = 0; l < ntypes; l++) if (MolName[ntypes * k + l] != Name[l]) flag_spec = 1; - if (flag_spec == 0) NMol[k]++; + if (flag_spec == 0) { + NMol[k]++; + Mol2Spec[m-1] = k; + } flag_identity *= flag_spec; } if (Nspec == 0 || flag_identity == 1) { for (l = 0; l < ntypes; l++) MolName[ntypes * Nspec + l] = Name[l]; + Mol2Spec[m-1] = Nspec; Nspec++; } } @@ -758,6 +817,169 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) /* ---------------------------------------------------------------------- */ +void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) +{ + int i, j, m, n, k, itype, cid; + int ndel, ndelone, count, count_tmp; + int *Nameall; + int *mask = atom->mask; + double localmass, totalmass; + double **spec_atom = f_SPECBOND->array_atom; + std::string species_str; + + AtomVec *avec = atom->avec; + + mark = nullptr; + memory->create(mark, nlocal, "reaxff/species:mark"); + for (i = 0; i < nlocal; i++) mark[i] = 0; + + Nameall = nullptr; + memory->create(Nameall, ntypes, "reaxff/species:Nameall"); + + int ndelcomm; + if (masslimitflag) ndelcomm = Nspec; + else ndelcomm = ndelspec; + + double *deletecount; + memory->create(deletecount, ndelcomm, "reaxff/species:deletecount"); + for (i = 0; i < ndelcomm; i++) deletecount[i] = 0; + + int nmarklist; + int *marklist; + memory->create(marklist, nlocal, "reaxff/species:marklist"); + + for (m = 1; m <= Nmole; m++) { + localmass = totalmass = count = nmarklist = 0; + for (n = 0; n < ntypes; n++) Name[n] = 0; + + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + cid = nint(clusterID[i]); + if (cid == m) { + itype = atom->type[i]-1; + Name[itype]++; + count++; + marklist[nmarklist++] = i; + localmass += atom->mass[atom->type[i]]; + } + } + + MPI_Allreduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, world); + count = count_tmp; + + MPI_Allreduce(Name, Nameall, ntypes, MPI_INT, MPI_SUM, world); + for (n = 0; n < ntypes; n++) Name[n] = Nameall[n]; + + MPI_Allreduce(&localmass, &totalmass, 1 , MPI_DOUBLE, MPI_SUM, world); + + species_str = ""; + for (j = 0; j < ntypes; j++) { + if (Name[j] != 0) { + if (eletype) species_str += eletype[j]; + else species_str += ele[j]; + if (Name[j] != 1) species_str += fmt::format("{}",Name[j]); + } + } + + if (masslimitflag) { + + // find corresponding moltype + + if (totalmass > massmin && totalmass < massmax) { + for (j = 0; j < nmarklist; j++) { + mark[marklist[j]] = 1; + deletecount[Mol2Spec[m-1]] += 1.0 / (double) count; + } + } + } else { + if (count > 0) { + for (i = 0; i < ndelspec; i++) { + if (del_species[i] == species_str) { + for (j = 0; j < nmarklist; j++) { + mark[marklist[j]] = 1; + deletecount[i] += 1.0 / (double) count; + } + break; + } + } + } + } + } + + // delete atoms. loop in reverse order to avoid copying marked atoms + + ndel = ndelone = 0; + for (i = atom->nlocal-1; i >= 0; i--) { + if (mark[i] == 1) { + avec->copy(atom->nlocal-1,i,1); + atom->nlocal--; + ndelone++; + } + } + + MPI_Allreduce(&ndelone, &ndel, 1, MPI_INT, MPI_SUM, world); + + atom->natoms -= ndel; + + if (me == 0) MPI_Reduce(MPI_IN_PLACE, deletecount, ndelcomm, MPI_DOUBLE, MPI_SUM, 0, world); + else MPI_Reduce(deletecount, deletecount, ndelcomm, MPI_DOUBLE, MPI_SUM, 0, world); + + if (me == 0) { + if (masslimitflag) { + int printflag = 0; + for (int m = 0; m < Nspec; m++) { + if (deletecount[m] > 0) { + if (printflag == 0) { + fprintf(fdel, "Timestep %lld", update->ntimestep); + printflag = 1; + } + fprintf(fdel, " %g ", deletecount[m]); + for (j = 0; j < ntypes; j ++) { + int itemp = MolName[ntypes * m + j]; + if (itemp != 0) { + if (eletype) fprintf(fdel, "%s", eletype[j]); + else fprintf(fdel, "%c", ele[j]); + if (itemp != 1) fprintf(fdel, "%d", itemp); + } + } + } + } + if (printflag) { + fprintf(fdel, "\n"); + fflush(fdel); + } + } else { + int writeflag = 0; + for (i = 0; i < ndelspec; i++) + if (deletecount[i]) writeflag = 1; + + if (writeflag) { + fprintf(fdel, "%lld", update->ntimestep); + for (i = 0; i < ndelspec; i++) { + fprintf(fdel, "\t%g", deletecount[i]); + } + fprintf(fdel, "\n"); + fflush(fdel); + } + } + } + + if (ndel && (atom->map_style != Atom::MAP_NONE)) { + atom->nghost = 0; + atom->map_init(); + atom->map_set(); + } + + next_reneighbor = update->ntimestep; + + memory->destroy(Nameall); + memory->destroy(marklist); + memory->destroy(mark); + memory->destroy(deletecount); +} + +/* ---------------------------------------------------------------------- */ + double FixReaxFFSpecies::compute_vector(int n) { if (n == 0) return vector_nmole; diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index f441c3bc92..e272ff7d6c 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -45,19 +45,24 @@ class FixReaxFFSpecies : public Fix { protected: int me, nprocs, nmax, nlocal, ntypes, ntotal; - int nrepeat, nfreq, posfreq, compressed; + int nrepeat, nfreq, posfreq, compressed, ndelspec; int Nmoltype, vector_nmole, vector_nspec; - int *Name, *MolName, *NMol, *nd, *MolType, *molmap; + int *Name, *MolName, *NMol, *nd, *MolType, *molmap, *mark; + int *Mol2Spec; double *clusterID; AtomCoord *x0; double bg_cut; double **BOCut; - FILE *fp, *pos; + std::vector del_species; + + FILE *fp, *pos, *fdel; int eleflag, posflag, multipos, padflag, setupflag; - int singlepos_opened, multipos_opened; - char *ele, **eletype, *filepos; + int delflag, specieslistflag, masslimitflag; + double massmin, massmax; + int singlepos_opened, multipos_opened, del_opened; + char *ele, **eletype, *filepos, *filedel; void Output_ReaxFF_Bonds(bigint, FILE *); AtomCoord chAnchor(AtomCoord, AtomCoord); @@ -65,6 +70,7 @@ class FixReaxFFSpecies : public Fix { void SortMolecule(int &); void FindSpecies(int, int &); void WriteFormulas(int, int); + void DeleteSpecies(int, int); int CheckExistence(int, int); int nint(const double &);