diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index 9bf1a82933..a6da15b161 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:: @@ -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 ReaxFF pair_coeff command. Position files are not written by default. diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 8fa06cafb3..2c0775b9d9 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,9 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : nvalid = -1; 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); @@ -156,8 +159,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,14 +193,9 @@ 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]); - } - eleflag = 1; + for (int i = 0; i < ntypes; i++) + eletype[i] = arg[iarg + 1 + i]; + GetUniqueElements(); iarg += ntypes + 1; // delete species @@ -285,14 +282,6 @@ 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 (delflag && specieslistflag && masslimitflag) error->all(FLERR, "Incompatible combination fix reaxff/species command options"); @@ -312,7 +301,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : FixReaxFFSpecies::~FixReaxFFSpecies() { - memory->destroy(ele); memory->destroy(BOCut); memory->destroy(clusterID); memory->destroy(x0); @@ -358,7 +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++) + eletype[i] = reaxff->eletype[i+1]; + GetUniqueElements(); + } post_integrate(); } @@ -648,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); @@ -664,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; } @@ -677,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; @@ -693,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++; } @@ -708,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; @@ -733,6 +727,50 @@ int FixReaxFFSpecies::CheckExistence(int id, int ntypes) /* ---------------------------------------------------------------------- */ +void FixReaxFFSpecies::GetUniqueElements() +{ + eleflag = 1; + + // get unique 'element' labels + + 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]) { + skipflag = 1; + break; + } + if (skipflag) continue; + ueletype[nutypes++] = eletype[i]; + } + + // 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; + } +} + +/* ---------------------------------------------------------------------- */ + void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec) { int i, j, itemp; @@ -742,17 +780,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) { - if (eletype) - molname += eletype[j]; - else - molname += ele[j]; + molname += ueletype[j]; if (itemp != 1) molname += std::to_string(itemp); } } @@ -810,20 +845,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]; @@ -850,17 +885,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) { - if (eletype) - fprintf(pos, "%s", eletype[n]); - else - fprintf(pos, "%c", ele[n]); + fprintf(pos, "%s", ueletype[n].c_str()); if (Name[n] != 1) fprintf(pos, "%d", Name[n]); } } @@ -912,7 +944,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) @@ -944,13 +976,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; @@ -961,18 +993,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) { - if (eletype) - species_str += eletype[j]; - else - species_str += ele[j]; + species_str += ueletype[j]; if (Name[j] != 1) species_str += fmt::format("{}", Name[j]); } } @@ -1034,13 +1063,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) { - if (eletype) - fprintf(fdel, "%s", eletype[j]); - else - fprintf(fdel, "%c", ele[j]); + 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 8fcbb131a9..b9afc5466a 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -43,7 +43,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; @@ -62,7 +62,10 @@ 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 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); @@ -72,6 +75,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; 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 594e7e502d..7bd2c6efe1 100644 --- a/src/REAXFF/pair_reaxff.h +++ b/src/REAXFF/pair_reaxff.h @@ -53,6 +53,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];