From 492b0641f2335fbcafec80f4dcc45c9d4d61c76e Mon Sep 17 00:00:00 2001 From: jrgissing Date: Fri, 16 Jun 2023 20:50:22 -0400 Subject: [PATCH 1/9] record element list from pair_coeff --- src/REAXFF/pair_reaxff.cpp | 5 ++++- src/REAXFF/pair_reaxff.h | 1 + 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/REAXFF/pair_reaxff.cpp b/src/REAXFF/pair_reaxff.cpp index 1867aec81b..99f7510a49 100644 --- a/src/REAXFF/pair_reaxff.cpp +++ b/src/REAXFF/pair_reaxff.cpp @@ -297,14 +297,17 @@ void PairReaxFF::coeff(int nargs, char **args) } int n = atom->ntypes; + eletype.resize(n+1); // pair_coeff element map - for (int i = 3; i < nargs; i++) + for (int i = 3; i < nargs; i++) { + eletype[i-2] = args[i]; for (int j = 0; j < nreax_types; j++) if (utils::lowercase(args[i]) == utils::lowercase(api->system->reax_param.sbp[j].name)) { map[i-2] = j; itmp ++; } + } // error check if (itmp != n) diff --git a/src/REAXFF/pair_reaxff.h b/src/REAXFF/pair_reaxff.h index 926e0060c9..a1c7e1c223 100644 --- a/src/REAXFF/pair_reaxff.h +++ b/src/REAXFF/pair_reaxff.h @@ -54,6 +54,7 @@ class PairReaxFF : public Pair { int fixbond_flag, fixspecies_flag; int **tmpid; double **tmpbo, **tmpr; + std::vector eletype; ReaxFF::API *api; typedef double rvec[3]; From d18d7edad99babd9e3d3ae48e23c93f158532e0e Mon Sep 17 00:00:00 2001 From: jrgissing Date: Fri, 16 Jun 2023 21:23:46 -0400 Subject: [PATCH 2/9] reaxff/species: default elements from pair reaxff --- src/REAXFF/fix_reaxff_species.cpp | 46 ++++++++----------------------- src/REAXFF/fix_reaxff_species.h | 3 +- 2 files changed, 14 insertions(+), 35 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 8fa06cafb3..0569abe1a0 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -65,7 +65,7 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), Name(nullptr), MolName(nullptr), NMol(nullptr), nd(nullptr), MolType(nullptr), molmap(nullptr), mark(nullptr), Mol2Spec(nullptr), clusterID(nullptr), x0(nullptr), BOCut(nullptr), fp(nullptr), pos(nullptr), fdel(nullptr), delete_Tcount(nullptr), - ele(nullptr), eletype(nullptr), filepos(nullptr), filedel(nullptr) + filepos(nullptr), filedel(nullptr) { if (narg < 7) utils::missing_cmd_args(FLERR, "fix reaxff/species", error); @@ -84,6 +84,7 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : nvalid = -1; ntypes = atom->ntypes; + eletype.resize(ntypes); nevery = utils::inumeric(FLERR, arg[3], false, lmp); nrepeat = utils::inumeric(FLERR, arg[4], false, lmp); @@ -156,8 +157,7 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : for (int j = 1; j < np1; j++) BOCut[i][j] = bo_cut; // optional args - eletype = nullptr; - ele = filepos = filedel = nullptr; + filepos = filedel = nullptr; eleflag = posflag = padflag = 0; delflag = specieslistflag = masslimitflag = 0; delete_Nlimit = delete_Nsteps = 0; @@ -191,13 +191,8 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : if (iarg + ntypes + 1 > narg) utils::missing_cmd_args(FLERR, "fix reaxff/species element", error); - eletype = (char **) malloc(ntypes * sizeof(char *)); - int len; - for (int i = 0; i < ntypes; i++) { - len = strlen(arg[iarg + 1 + i]) + 1; - eletype[i] = (char *) malloc(len * sizeof(char)); - strcpy(eletype[i], arg[iarg + 1 + i]); - } + for (int i = 0; i < ntypes; i++) + eletype[i] = arg[iarg + 1 + i]; eleflag = 1; iarg += ntypes + 1; @@ -285,13 +280,9 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Unknown fix reaxff/species keyword: {}", arg[iarg]); } - if (!eleflag) { - memory->create(ele, ntypes + 1, "reaxff/species:ele"); - ele[0] = 'C'; - if (ntypes > 1) ele[1] = 'H'; - if (ntypes > 2) ele[2] = 'O'; - if (ntypes > 3) ele[3] = 'N'; - } + if (!eleflag) + for (int i = 0; i < ntypes; i++) + eletype[i] = reaxff->eletype[i+1]; if (delflag && specieslistflag && masslimitflag) error->all(FLERR, "Incompatible combination fix reaxff/species command options"); @@ -312,7 +303,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : FixReaxFFSpecies::~FixReaxFFSpecies() { - memory->destroy(ele); memory->destroy(BOCut); memory->destroy(clusterID); memory->destroy(x0); @@ -749,10 +739,7 @@ void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec) for (j = 0; j < ntypes; j++) { itemp = MolType[ntypes * i + j]; if (itemp != 0) { - if (eletype) - molname += eletype[j]; - else - molname += ele[j]; + molname += eletype[j]; if (itemp != 1) molname += std::to_string(itemp); } } @@ -857,10 +844,7 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) fprintf(pos, "%d\t%d\t", m, count); for (n = 0; n < ntypes; n++) { if (Name[n] != 0) { - if (eletype) - fprintf(pos, "%s", eletype[n]); - else - fprintf(pos, "%c", ele[n]); + fprintf(pos, "%s", eletype[n].c_str()); if (Name[n] != 1) fprintf(pos, "%d", Name[n]); } } @@ -969,10 +953,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) species_str = ""; for (j = 0; j < ntypes; j++) { if (Name[j] != 0) { - if (eletype) - species_str += eletype[j]; - else - species_str += ele[j]; + species_str += eletype[j]; if (Name[j] != 1) species_str += fmt::format("{}", Name[j]); } } @@ -1037,10 +1018,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) 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]); + fprintf(fdel, "%s", eletype[j].c_str()); if (itemp != 1) fprintf(fdel, "%d", itemp); } } diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index f711cdeb11..829ab61cd8 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -63,7 +63,8 @@ class FixReaxFFSpecies : public Fix { int delete_Nsteps, *delete_Tcount; double massmin, massmax; int singlepos_opened, multipos_opened, del_opened; - char *ele, **eletype, *filepos, *filedel; + char *filepos, *filedel; + std::vector eletype; void Output_ReaxFF_Bonds(bigint, FILE *); AtomCoord chAnchor(AtomCoord, AtomCoord); From 32e4aac9f73d72be28c21a7a3c0007f3ddf8d466 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Fri, 16 Jun 2023 21:34:06 -0400 Subject: [PATCH 3/9] Update fix_reaxff_species.rst --- doc/src/fix_reaxff_species.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index 6a171ede5e..9c3d1d01fd 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -110,10 +110,10 @@ all types from 1 to :math:`N`. A leading asterisk means all types from The optional keyword *element* can be used to specify the chemical symbol printed for each LAMMPS atom type. The number of symbols must match the number of LAMMPS atom types and each symbol must consist of -1 or 2 alphanumeric characters. Normally, these symbols should be -chosen to match the chemical identity of each LAMMPS atom type, as -specified using the :doc:`reaxff pair_coeff ` command and -the ReaxFF force field file. +1 or 2 alphanumeric characters. By default, these symbols are the same +as the chemical identity of each LAMMPS atom type, as specified by the +:doc:`reaxff pair_coeff ` command and the ReaxFF force +field file. The optional keyword *position* writes center-of-mass positions of each identified molecules to file *filepos* every *posfreq* timesteps. @@ -233,5 +233,5 @@ Default """"""" The default values for bond-order cutoffs are 0.3 for all I-J pairs. -The default element symbols are C, H, O, N. +The default element symbols are taken from the pair_coeff command. Position files are not written by default. From 33d82c30ca146b71ec238eb343ef175a9653de1e Mon Sep 17 00:00:00 2001 From: jrgissing Date: Thu, 22 Jun 2023 18:40:50 -0400 Subject: [PATCH 4/9] fix to allow reaxff/species before pair_coeff --- src/REAXFF/fix_reaxff_species.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 0569abe1a0..98c4280eef 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -280,10 +280,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Unknown fix reaxff/species keyword: {}", arg[iarg]); } - if (!eleflag) - for (int i = 0; i < ntypes; i++) - eletype[i] = reaxff->eletype[i+1]; - if (delflag && specieslistflag && masslimitflag) error->all(FLERR, "Incompatible combination fix reaxff/species command options"); @@ -350,6 +346,12 @@ void FixReaxFFSpecies::setup(int /*vflag*/) ntotal = static_cast(atom->natoms); if (Name == nullptr) memory->create(Name, ntypes, "reaxff/species:Name"); + if (!eleflag) + for (int i = 0; i < ntypes; i++) { + eletype[i] = reaxff->eletype[i+1]; + eleflag = 1; + } + post_integrate(); } From 6318b09a07d04c5cad74a06d3075eb266dcd8c8f Mon Sep 17 00:00:00 2001 From: jrgissing Date: Tue, 4 Jul 2023 20:29:31 -0400 Subject: [PATCH 5/9] report unique species when duplicate elements previously, duplicate species were reported when there were duplicate elements in the element-to-type mapping. for example, H2 and HH and multiple other H2s and HHs could appear in the 'species' and 'delete species' files --- src/REAXFF/fix_reaxff_species.cpp | 110 ++++++++++++++++++------------ src/REAXFF/fix_reaxff_species.h | 5 +- 2 files changed, 72 insertions(+), 43 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 98c4280eef..fade2c97d3 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -85,6 +85,8 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : ntypes = atom->ntypes; eletype.resize(ntypes); + ueletype.resize(ntypes); + ele2uele.resize(ntypes); nevery = utils::inumeric(FLERR, arg[3], false, lmp); nrepeat = utils::inumeric(FLERR, arg[4], false, lmp); @@ -193,7 +195,7 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < ntypes; i++) eletype[i] = arg[iarg + 1 + i]; - eleflag = 1; + GetUniqueElements(); iarg += ntypes + 1; // delete species @@ -344,13 +346,13 @@ int FixReaxFFSpecies::setmask() void FixReaxFFSpecies::setup(int /*vflag*/) { ntotal = static_cast(atom->natoms); - if (Name == nullptr) memory->create(Name, ntypes, "reaxff/species:Name"); + if (Name == nullptr) memory->create(Name, nutypes, "reaxff/species:Name"); - if (!eleflag) - for (int i = 0; i < ntypes; i++) { + if (!eleflag) { + for (int i = 0; i < ntypes; i++) eletype[i] = reaxff->eletype[i+1]; - eleflag = 1; - } + GetUniqueElements(); + } post_integrate(); } @@ -640,14 +642,14 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) memory->destroy(MolName); MolName = nullptr; - memory->create(MolName, Nmole * (ntypes + 1), "reaxff/species:MolName"); + memory->create(MolName, Nmole * (nutypes + 1), "reaxff/species:MolName"); memory->destroy(NMol); NMol = nullptr; memory->create(NMol, Nmole, "reaxff/species:NMol"); for (m = 0; m < Nmole; m++) NMol[m] = 1; - memory->create(Nameall, ntypes, "reaxff/species:Nameall"); + memory->create(Nameall, nutypes, "reaxff/species:Nameall"); memory->create(NMolall, Nmole, "reaxff/species:NMolall"); memory->destroy(Mol2Spec); @@ -656,12 +658,12 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) 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; n < nutypes; n++) Name[n] = 0; for (n = 0, flag_mol = 0; n < nlocal; n++) { if (!(mask[n] & groupbit)) continue; cid = nint(clusterID[n]); if (cid == m) { - itype = atom->type[n] - 1; + itype = ele2uele[atom->type[n] - 1]; Name[itype]++; flag_mol = 1; } @@ -669,15 +671,15 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) MPI_Allreduce(&flag_mol, &flag_tmp, 1, MPI_INT, MPI_MAX, world); flag_mol = flag_tmp; - MPI_Allreduce(Name, Nameall, ntypes, MPI_INT, MPI_SUM, world); - for (n = 0; n < ntypes; n++) Name[n] = Nameall[n]; + MPI_Allreduce(Name, Nameall, nutypes, MPI_INT, MPI_SUM, world); + for (n = 0; n < nutypes; n++) Name[n] = Nameall[n]; if (flag_mol == 1) { flag_identity = 1; for (k = 0; k < Nspec; k++) { flag_spec = 0; - for (l = 0; l < ntypes; l++) - if (MolName[ntypes * k + l] != Name[l]) flag_spec = 1; + for (l = 0; l < nutypes; l++) + if (MolName[nutypes * k + l] != Name[l]) flag_spec = 1; if (flag_spec == 0) { NMol[k]++; Mol2Spec[m - 1] = k; @@ -685,7 +687,7 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) flag_identity *= flag_spec; } if (Nspec == 0 || flag_identity == 1) { - for (l = 0; l < ntypes; l++) MolName[ntypes * Nspec + l] = Name[l]; + for (l = 0; l < nutypes; l++) MolName[nutypes * Nspec + l] = Name[l]; Mol2Spec[m - 1] = Nspec; Nspec++; } @@ -700,24 +702,24 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec) memory->destroy(MolType); MolType = nullptr; - memory->create(MolType, Nspec * (ntypes + 2), "reaxff/species:MolType"); + memory->create(MolType, Nspec * (nutypes + 2), "reaxff/species:MolType"); } /* ---------------------------------------------------------------------- */ -int FixReaxFFSpecies::CheckExistence(int id, int ntypes) +int FixReaxFFSpecies::CheckExistence(int id, int nutypes) { int i, j, molid, flag; for (i = 0; i < Nmoltype; i++) { flag = 0; - for (j = 0; j < ntypes; j++) { - molid = MolType[ntypes * i + j]; - if (molid != MolName[ntypes * id + j]) flag = 1; + for (j = 0; j < nutypes; j++) { + molid = MolType[nutypes * i + j]; + if (molid != MolName[nutypes * id + j]) flag = 1; } if (flag == 0) return i; } - for (i = 0; i < ntypes; i++) MolType[ntypes * Nmoltype + i] = MolName[ntypes * id + i]; + for (i = 0; i < nutypes; i++) MolType[nutypes * Nmoltype + i] = MolName[nutypes * id + i]; Nmoltype++; return Nmoltype - 1; @@ -725,6 +727,30 @@ int FixReaxFFSpecies::CheckExistence(int id, int ntypes) /* ---------------------------------------------------------------------- */ +void FixReaxFFSpecies::GetUniqueElements() +{ + // count unique 'element' labels + // map user input to unique list + + nutypes = 0; + int skipflag; + for (int i = 0; i < ntypes; i++) { + skipflag = 0; + for (int j = 0; j < nutypes; j++) + if (eletype[i] == ueletype[j]) { + ele2uele[i] = j; + skipflag = 1; + break; + } + if (skipflag) continue; + ele2uele[i] = nutypes; + ueletype[nutypes++] = eletype[i]; + } + eleflag = 1; +} + +/* ---------------------------------------------------------------------- */ + void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec) { int i, j, itemp; @@ -734,14 +760,14 @@ void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec) Nmoltype = 0; - for (i = 0; i < Nspec; i++) nd[i] = CheckExistence(i, ntypes); + for (i = 0; i < Nspec; i++) nd[i] = CheckExistence(i, nutypes); for (i = 0; i < Nmoltype; i++) { std::string molname; - for (j = 0; j < ntypes; j++) { - itemp = MolType[ntypes * i + j]; + for (j = 0; j < nutypes; j++) { + itemp = MolType[nutypes * i + j]; if (itemp != 0) { - molname += eletype[j]; + molname += ueletype[j]; if (itemp != 1) molname += std::to_string(itemp); } } @@ -799,20 +825,20 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) } Nameall = nullptr; - memory->create(Nameall, ntypes, "reaxff/species:Nameall"); + memory->create(Nameall, nutypes, "reaxff/species:Nameall"); for (m = 1; m <= Nmole; m++) { count = 0; avq = 0.0; for (n = 0; n < 3; n++) avx[n] = 0.0; - for (n = 0; n < ntypes; n++) Name[n] = 0; + for (n = 0; n < nutypes; 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; + itype = ele2uele[atom->type[i] - 1]; Name[itype]++; count++; avq += spec_atom[i][0]; @@ -839,14 +865,14 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) MPI_Reduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, 0, world); count = count_tmp; - MPI_Reduce(Name, Nameall, ntypes, MPI_INT, MPI_SUM, 0, world); - for (n = 0; n < ntypes; n++) Name[n] = Nameall[n]; + MPI_Reduce(Name, Nameall, nutypes, MPI_INT, MPI_SUM, 0, world); + for (n = 0; n < nutypes; n++) Name[n] = Nameall[n]; if (comm->me == 0) { fprintf(pos, "%d\t%d\t", m, count); - for (n = 0; n < ntypes; n++) { + for (n = 0; n < nutypes; n++) { if (Name[n] != 0) { - fprintf(pos, "%s", eletype[n].c_str()); + fprintf(pos, "%s", ueletype[n].c_str()); if (Name[n] != 1) fprintf(pos, "%d", Name[n]); } } @@ -898,7 +924,7 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) for (i = 0; i < nlocal; i++) mark[i] = 0; Nameall = nullptr; - memory->create(Nameall, ntypes, "reaxff/species:Nameall"); + memory->create(Nameall, nutypes, "reaxff/species:Nameall"); int ndelcomm; if (masslimitflag) @@ -930,13 +956,13 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) if (this_delete_Tcount == headroom) break; m = molrange[mm]; localmass = totalmass = count = nmarklist = 0; - for (n = 0; n < ntypes; n++) Name[n] = 0; + for (n = 0; n < nutypes; 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; + itype = ele2uele[atom->type[i] - 1]; Name[itype]++; count++; marklist[nmarklist++] = i; @@ -947,15 +973,15 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) 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(Name, Nameall, nutypes, MPI_INT, MPI_SUM, world); + for (n = 0; n < nutypes; n++) Name[n] = Nameall[n]; MPI_Allreduce(&localmass, &totalmass, 1, MPI_DOUBLE, MPI_SUM, world); species_str = ""; - for (j = 0; j < ntypes; j++) { + for (j = 0; j < nutypes; j++) { if (Name[j] != 0) { - species_str += eletype[j]; + species_str += ueletype[j]; if (Name[j] != 1) species_str += fmt::format("{}", Name[j]); } } @@ -1017,10 +1043,10 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) printflag = 1; } fprintf(fdel, " %g ", deletecount[m]); - for (j = 0; j < ntypes; j++) { - int itemp = MolName[ntypes * m + j]; + for (j = 0; j < nutypes; j++) { + int itemp = MolName[nutypes * m + j]; if (itemp != 0) { - fprintf(fdel, "%s", eletype[j].c_str()); + fprintf(fdel, "%s", ueletype[j].c_str()); if (itemp != 1) fprintf(fdel, "%d", itemp); } } diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index 829ab61cd8..91cdad8246 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -44,7 +44,7 @@ class FixReaxFFSpecies : public Fix { double compute_vector(int) override; protected: - int nmax, nlocal, ntypes, ntotal; + int nmax, nlocal, ntypes, nutypes, ntotal; int nrepeat, nfreq, posfreq, compressed, ndelspec; int Nmoltype, vector_nmole, vector_nspec; int *Name, *MolName, *NMol, *nd, *MolType, *molmap, *mark; @@ -64,7 +64,9 @@ class FixReaxFFSpecies : public Fix { double massmin, massmax; int singlepos_opened, multipos_opened, del_opened; char *filepos, *filedel; + std::vector ele2uele; std::vector eletype; + std::vector ueletype; void Output_ReaxFF_Bonds(bigint, FILE *); AtomCoord chAnchor(AtomCoord, AtomCoord); @@ -74,6 +76,7 @@ class FixReaxFFSpecies : public Fix { void WriteFormulas(int, int); void DeleteSpecies(int, int); int CheckExistence(int, int); + void GetUniqueElements(); int nint(const double &); int pack_forward_comm(int, int *, double *, int, int *) override; From 4eba3791f3b201e7c9c3f07fcf84955b3031be53 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 9 Jul 2023 15:30:56 -0400 Subject: [PATCH 6/9] enforce order for printing CHON --- src/REAXFF/fix_reaxff_species.cpp | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index fade2c97d3..2c0775b9d9 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -729,8 +729,9 @@ int FixReaxFFSpecies::CheckExistence(int id, int nutypes) void FixReaxFFSpecies::GetUniqueElements() { - // count unique 'element' labels - // map user input to unique list + eleflag = 1; + + // get unique 'element' labels nutypes = 0; int skipflag; @@ -738,15 +739,34 @@ void FixReaxFFSpecies::GetUniqueElements() skipflag = 0; for (int j = 0; j < nutypes; j++) if (eletype[i] == ueletype[j]) { - ele2uele[i] = j; skipflag = 1; break; } if (skipflag) continue; - ele2uele[i] = nutypes; ueletype[nutypes++] = eletype[i]; } - eleflag = 1; + + // reorder CHON, if necessary + + int incr = 0; + std::vector CHON = {"C", "H", "O", "N"}; + for (auto it = CHON.begin(); it != CHON.end(); ++it) + for (int j = incr; j < nutypes; j++) { + if (ueletype[j] == *it) { + ueletype.erase(ueletype.begin() + j); + ueletype.insert(ueletype.begin() + incr++, *it); + break; + } + } + + // map user input to unique list + + for (int i = 0; i < ntypes; i++) + for (int j = 0; j < nutypes; j++) + if (eletype[i] == ueletype[j]) { + ele2uele[i] = j; + break; + } } /* ---------------------------------------------------------------------- */ From ada61d96fe1fdf8f708c09e3ecadaeb4e4e6ecc7 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 8 May 2024 19:46:07 -0400 Subject: [PATCH 7/9] Update fix_reaxff_species.rst --- doc/src/fix_reaxff_species.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index bc8e9cf27f..76ecc934ff 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* or *delete* +* keyword = *cutoff* or *element* or *position* or *delete* or *delete_rate_limit* .. parsed-literal:: @@ -233,5 +233,5 @@ Default """"""" The default values for bond-order cutoffs are 0.3 for all I-J pairs. -The default element symbols are taken from the pair_coeff command. +The default element symbols are taken from the ReaxFF pair_coeff command. Position files are not written by default. From 8e6a232dffde3112401dd910a69f18b35a9229f2 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 8 May 2024 23:31:36 -0400 Subject: [PATCH 8/9] Update fix_reaxff_species.rst --- doc/src/fix_reaxff_species.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index 76ecc934ff..a6da15b161 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -112,7 +112,7 @@ symbol printed for each LAMMPS atom type. The number of symbols must match the number of LAMMPS atom types and each symbol must consist of 1 or 2 alphanumeric characters. By default, these symbols are the same as the chemical identity of each LAMMPS atom type, as specified by the -:doc:`reaxff pair_coeff ` command and the ReaxFF force +:doc:`ReaxFF pair_coeff ` command and the ReaxFF force field file. The optional keyword *position* writes center-of-mass positions of From b4f18700dc6dd20df812999c4f9398131b0c4c93 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Thu, 9 May 2024 00:06:10 -0400 Subject: [PATCH 9/9] Update fix_reaxff_species.h --- src/REAXFF/fix_reaxff_species.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index 1968d413ca..b9afc5466a 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -63,9 +63,9 @@ class FixReaxFFSpecies : public Fix { double massmin, massmax; int singlepos_opened, multipos_opened, del_opened; char *filepos, *filedel; - std::vector ele2uele; - std::vector eletype; - std::vector ueletype; + std::vector ele2uele; // for element eletype[i], ele2uele[i] stores index of unique element + std::vector eletype; // list of ReaxFF elements of length ntypes + std::vector ueletype; // list of unique elements, of quantity nutypes void Output_ReaxFF_Bonds(bigint, FILE *); AtomCoord chAnchor(AtomCoord, AtomCoord);