diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index 107695f0f6e..068d1de9954 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -200,8 +200,8 @@ The 2 values in the global vector are as follows: The per-atom vector stores the molecule ID for each atom as identified by the fix. If an atom is not in a molecule, its ID will be 0. For atoms in the same molecule, the molecule ID for all of them -will be the same and will be equal to the smallest atom ID of -any atom in the molecule. +will be the same, and molecule IDs will range from 1 to the number +of molecules. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 1916597a692..09c3884c342 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -26,6 +26,7 @@ #include "domain.h" #include "error.h" #include "fix_ave_atom.h" +#include "fix_property_atom.h" #include "force.h" #include "group.h" #include "input.h" @@ -141,13 +142,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : } x0 = nullptr; - clusterID = nullptr; - - int ntmp = atom->nmax; - memory->create(x0, ntmp, "reaxff/species:x0"); - memory->create(clusterID, ntmp, "reaxff/species:clusterID"); - memset(clusterID, 0, sizeof(double) * ntmp); - vector_atom = clusterID; nmax = 0; setupflag = 0; @@ -304,7 +298,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : FixReaxFFSpecies::~FixReaxFFSpecies() { memory->destroy(BOCut); - memory->destroy(clusterID); memory->destroy(x0); memory->destroy(nd); @@ -330,6 +323,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies() try { modify->delete_compute(fmt::format("SPECATOM_{}", id)); modify->delete_fix(fmt::format("SPECBOND_{}", id)); + modify->delete_fix(fmt::format("clusterID_{}", id)); } catch (std::exception &) { } } @@ -372,9 +366,6 @@ void FixReaxFFSpecies::init() reaxff->fixspecies_flag = 1; - // reset next output timestep if not yet set or timestep has been reset - if (nvalid != update->ntimestep) nvalid = update->ntimestep + nfreq; - if (!setupflag) { // create a compute to store properties modify->add_compute(fmt::format("SPECATOM_{} all SPEC/ATOM q x y z vx vy vz abo01 abo02 " @@ -387,11 +378,25 @@ void FixReaxFFSpecies::init() auto fixcmd = fmt::format("SPECBOND_{} all ave/atom {} {} {}", id, nevery, nrepeat, nfreq); for (int i = 1; i < 32; ++i) fixcmd += fmt::format(" c_SPECATOM_{}[{}]", id, i); f_SPECBOND = dynamic_cast(modify->add_fix(fixcmd)); + + // create a fix to point to fix_property_atom for storing clusterID + fixcmd = fmt::format("clusterID_{} all property/atom d_clusterID ghost yes", id); + f_clusterID = dynamic_cast(modify->add_fix(fixcmd)); + + // per-atom property for clusterID + int flag,cols; + int index1 = atom->find_custom("clusterID",flag,cols); + clusterID = atom->dvector[index1]; + vector_atom = clusterID; + + int ntmp = atom->nmax; + memory->create(x0, ntmp, "reaxff/species:x0"); + setupflag = 1; } // check for valid variable name for delete Nlimit keyword - if (delete_Nsteps > 0) { + if (delete_Nsteps > 0 && delete_Nlimit_varid > -1) { delete_Nlimit_varid = input->variable->find(delete_Nlimit_varname.c_str()); if (delete_Nlimit_varid < 0) error->all(FLERR, "Fix reaxff/species: Variable name {} does not exist", @@ -424,10 +429,16 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) { int Nmole, Nspec; + // per-atom property for clusterID + int flag,cols; + int index1 = atom->find_custom("clusterID",flag,cols); + clusterID = atom->dvector[index1]; + vector_atom = clusterID; + // point to fix_ave_atom f_SPECBOND->end_of_step(); - if (ntimestep != nvalid) { + if (ntimestep != nvalid && nvalid != -1) { // push back delete_Tcount on every step if (delete_Nsteps > 0) for (int i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; @@ -439,11 +450,7 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) if (atom->nmax > nmax) { nmax = atom->nmax; memory->destroy(x0); - memory->destroy(clusterID); memory->create(x0, nmax, "reaxff/species:x0"); - memory->create(clusterID, nmax, "reaxff/species:clusterID"); - memset(clusterID, 0, sizeof(double) * nmax); - vector_atom = clusterID; } for (int i = 0; i < nmax; i++) { x0[i].x = x0[i].y = x0[i].z = 0.0; } @@ -464,9 +471,14 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) if (comm->me == 0) fflush(pos); } - if (delflag) DeleteSpecies(Nmole, Nspec); + if (delflag && nvalid != -1) { + DeleteSpecies(Nmole, Nspec); + + // reset molecule ID to index from 1 + SortMolecule(Nmole); + } - nvalid += nfreq; + nvalid = ntimestep + nfreq; } /* ---------------------------------------------------------------------- */ @@ -826,7 +838,8 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) int count, count_tmp, m, n, k; int *Nameall; int *mask = atom->mask; - double avq, avq_tmp, avx[3], avx_tmp, box[3], halfbox[3]; + double *rmass = atom->rmass; + double totq, totq_tmp, com[3], com_tmp, thism, totm, box[3], halfbox[3]; double **spec_atom = f_SPECBOND->array_atom; if (multipos) OpenPos(); @@ -844,7 +857,7 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) update->ntimestep, Nmole, Nspec, domain->boxlo[0], domain->boxhi[0], domain->boxlo[1], domain->boxhi[1], domain->boxlo[2], domain->boxhi[2]); - fprintf(pos, "ID\tAtom_Count\tType\tAve_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n"); + fprintf(pos, "ID\tAtom_Count\tType\tTot_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n"); } Nameall = nullptr; @@ -853,8 +866,9 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) for (m = 1; m <= Nmole; m++) { count = 0; - avq = 0.0; - for (n = 0; n < 3; n++) avx[n] = 0.0; + totq = 0.0; + totm = 0.0; + for (n = 0; n < 3; n++) com[n] = 0.0; for (n = 0; n < nutypes; n++) Name[n] = 0; for (i = 0; i < nlocal; i++) { @@ -864,30 +878,37 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) itype = ele2uele[atom->type[i] - 1]; Name[itype]++; count++; - avq += spec_atom[i][0]; + totq += spec_atom[i][0]; if ((x0[i].x - spec_atom[i][1]) > halfbox[0]) spec_atom[i][1] += box[0]; if ((spec_atom[i][1] - x0[i].x) > halfbox[0]) spec_atom[i][1] -= box[0]; if ((x0[i].y - spec_atom[i][2]) > halfbox[1]) spec_atom[i][2] += box[1]; if ((spec_atom[i][2] - x0[i].y) > halfbox[1]) spec_atom[i][2] -= box[1]; if ((x0[i].z - spec_atom[i][3]) > halfbox[2]) spec_atom[i][3] += box[2]; if ((spec_atom[i][3] - x0[i].z) > halfbox[2]) spec_atom[i][3] -= box[2]; - for (n = 0; n < 3; n++) avx[n] += spec_atom[i][n + 1]; + if (rmass) thism = rmass[i]; + else thism = atom->mass[atom->type[i]]; + for (n = 0; n < 3; n++) com[n] += spec_atom[i][n+1]*thism; + totm += thism; } } - avq_tmp = 0.0; - MPI_Allreduce(&avq, &avq_tmp, 1, MPI_DOUBLE, MPI_SUM, world); - avq = avq_tmp; + totq_tmp = 0.0; + MPI_Allreduce(&totq, &totq_tmp, 1, MPI_DOUBLE, MPI_SUM, world); + totq = totq_tmp; for (n = 0; n < 3; n++) { - avx_tmp = 0.0; - MPI_Reduce(&avx[n], &avx_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); - avx[n] = avx_tmp; + com_tmp = 0.0; + MPI_Reduce(&com[n], &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); + com[n] = com_tmp; } MPI_Reduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, 0, world); count = count_tmp; + com_tmp = 0.0; + MPI_Reduce(&totm, &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); + totm = com_tmp; + MPI_Reduce(Name, Nameall, nutypes, MPI_INT, MPI_SUM, 0, world); for (n = 0; n < nutypes; n++) Name[n] = Nameall[n]; @@ -900,16 +921,15 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) } } if (count > 0) { - avq /= count; for (k = 0; k < 3; k++) { - avx[k] /= count; - if (avx[k] >= domain->boxhi[k]) avx[k] -= box[k]; - if (avx[k] < domain->boxlo[k]) avx[k] += box[k]; + com[k] /= totm; + if (com[k] >= domain->boxhi[k]) com[k] -= box[k]; + if (com[k] < domain->boxlo[k]) com[k] += box[k]; - avx[k] -= domain->boxlo[k]; - avx[k] /= box[k]; + com[k] -= domain->boxlo[k]; + com[k] /= box[k]; } - fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", avq, avx[0], avx[1], avx[2]); + fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", totq, com[0], com[1], com[2]); } fprintf(pos, "\n"); } @@ -922,21 +942,29 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) { - int ndeletions; + int i, ndeletions; int headroom = -1; if (delete_Nsteps > 0) { - if (delete_Tcount[delete_Nsteps - 1] == -1) return; + if (delete_Tcount[delete_Nsteps - 1] == -1) { + for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; + return; + } ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps - 1]; if (delete_Nlimit_varid > -1) delete_Nlimit = input->variable->compute_equal(delete_Nlimit_varid); headroom = MAX(0, delete_Nlimit - ndeletions); - if (headroom == 0) return; + if (headroom == 0) { + for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; + return; + } } - int i, j, m, n, itype, cid; + int j, m, n, itype, cid; int ndel, ndelone, count, count_tmp; int *Nameall; int *mask = atom->mask; + double *mass = atom->mass; + double *rmass = atom->rmass; double localmass, totalmass; std::string species_str; @@ -989,7 +1017,8 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) Name[itype]++; count++; marklist[nmarklist++] = i; - localmass += atom->mass[atom->type[i]]; + if (rmass) localmass += rmass[i]; + else localmass += atom->mass[atom->type[i]]; } } @@ -1177,7 +1206,7 @@ double FixReaxFFSpecies::memory_usage() { double bytes; - bytes = 4 * nmax * sizeof(double); // clusterID + x0 + bytes = 3 * nmax * sizeof(double); // x0 return bytes; } diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index b9afc5466a8..d378065a82e 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -88,6 +88,7 @@ class FixReaxFFSpecies : public Fix { class NeighList *list; class FixAveAtom *f_SPECBOND; + class FixPropertyAtom *f_clusterID; class PairReaxFF *reaxff; }; } // namespace LAMMPS_NS