ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3298 from jrgissing/reaxff_species_delete_keyword

Browse Source 
Reaxff species 'delete' keyword
This commit is contained in:
Axel Kohlmeyer
2022-07-07 16:37:43 -04:00
committed by GitHub
parent 6280c81151 1599afc623
commit 7d06ee9aa7
3 changed files with 283 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
doc/src/fix_reaxff_species.rst
View File

@ -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

230
src/REAXFF/fix_reaxff_species.cpp
View File

@ -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;

16
src/REAXFF/fix_reaxff_species.h
View File

@ -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<std::string> 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 &);