diff --git a/doc/src/fix_latte.rst b/doc/src/fix_latte.rst index 6ce84f9d52..df23027238 100644 --- a/doc/src/fix_latte.rst +++ b/doc/src/fix_latte.rst @@ -19,8 +19,8 @@ Syntax keyword = *coulomb* or *exclude* *coulomb* value = peID peID = ID of compute used to calculate per-atom energy - *exclude* value = groupID - groupID = ID of group of atoms to exclude before calling LATTE + *exclude* value = fixID + fixID = ID of fix which potentially excludes atoms before calling LATTE Examples """""""" @@ -28,7 +28,7 @@ Examples .. code-block:: LAMMPS fix dftb all latte - fix dftb all exclude GCMCmol + fix dftb all exclude GCMC Description """"""""""" @@ -64,27 +64,17 @@ The *coulomb* argument is not yet supported by fix latte (as of Sept Coulomb potential as an alternative to LATTE performing the calculation. -NOTE: after intitial debugging, change the exclude arg to -be the ID of another fix (GCMC in this case), and extract() -the exclusion group ID from fix gcmc. +The *exclude argument allows this fix to work in tandem with another +fix which may decide to delete one or more atoms of molecules. The +specified fixID is the ID of the other fix. -The *exclude argument allows this fix to work in tandem with the -:doc:`fix gcmc ` command which may decide to delete an atom -or molecule as one of its Monte Carlo events. In this case, LAMMPS -needs to pass LATTE the atoms for the system with the atom/molecule -removed. Fix gcmc does not actually remove the atom/molecule until -after the new energy is computed (in this case by LATTE), and a Monte -Carlo accept/reject decision is made for the event. - -The specified groupID must match the group ID which the :doc:`fix gcmc -` command assigns to atoms flagged for possible deletion. -It should be either its default exclusion group ID or group ID used -with its "exclude" keyword option. - -.. note:: - - The fix gcmc command must appear in the input script prior - to the fix latte command for this to work. +The one current example of such a fix is the :doc:`fix gcmc +` command which performs Monte Carlo insertions and +deletions. If a trial deletion is performed, then LAMMPS needs to +only pass LATTE the atoms which remain. Fix gcmc does not actually +remove any atoms until after the new energy is computed (in this case +by LATTE), and a Monte Carlo accept/reject decision is made for the +trial deletion. ---------- diff --git a/src/LATTE/fix_latte.cpp b/src/LATTE/fix_latte.cpp index 1095b9ab9a..0ba16a9c61 100644 --- a/src/LATTE/fix_latte.cpp +++ b/src/LATTE/fix_latte.cpp @@ -76,6 +76,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : coulomb = 0; id_pe = nullptr; exclude = 0; + id_exclude = nullptr; int iarg = 3; while (iarg < narg) { @@ -85,7 +86,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : coulomb = 1; error->all(FLERR,"Fix latte does not yet support a " "LAMMPS calculation of a Coulomb potential"); - id_pe = utils::strdup(arg[3]); + id_pe = utils::strdup(arg[iarg+1]); c_pe = modify->get_compute_by_id(id_pe); if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); if (c_pe->peatomflag == 0) @@ -96,11 +97,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix latte exclude", error); exclude = 1; - int excludegroup = group->find(arg[iarg+1]); - if (excludegroup == -1) - error->all(FLERR, "Fix latte couldd not find exclude group ID: {}", - arg[iarg+1]); - excludebit = group->bitmask[excludegroup]; + id_exclude = utils::strdup(arg[iarg+1]); iarg += 2; } else @@ -121,6 +118,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : FixLatte::~FixLatte() { delete[] id_pe; + delete[] id_exclude; memory->destroy(qpotential); memory->destroy(flatte); } @@ -148,9 +146,9 @@ void FixLatte::init() if (coulomb) { if (atom->q_flag == 0 || force->pair == nullptr || force->kspace == nullptr) error->all(FLERR,"Fix latte cannot compute Coulomb potential"); - c_pe = modify->get_compute_by_id(id_pe); - if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); + if (!c_pe) + error->all(FLERR,"Fix latte could not find Coulomb compute ID {}",id_pe); } // must be fully periodic or fully non-periodic @@ -167,6 +165,20 @@ void FixLatte::init() memory->create(qpotential,atom->nlocal,"latte:qpotential"); memory->create(flatte,atom->nlocal,3,"latte:flatte"); } + + // extract pointer to exclusion_group variable from id_exclude + // exclusion_group is index of a group the Fix defines + + if (exclude) { + Fix *f_exclude = modify->get_fix_by_id(id_exclude); + if (!f_exclude) + error->all(FLERR,"Fix latte could not find exclude fix ID {}", id_exclude); + int exclude_group_index,dim; + exclusion_group_ptr = (int *) f_exclude->extract("exclusion_group",dim); + if (!exclusion_group_ptr || dim != 0) + error->all(FLERR,"Fix latte could not query exclude_group of fix ID {}", + id_exclude); + } } /* ---------------------------------------------------------------------- */ @@ -267,12 +279,18 @@ void FixLatte::post_force(int vflag) if (!exclude) latte_wrapper_all(); else { - int *mask = atom->mask; - int nlocal = atom->nlocal; - int anyexclude = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & excludebit) anyexclude = 1; + + int exclusion_group = *exclusion_group_ptr; + if (exclusion_group) { + int excludebit = group->bitmask[exclusion_group]; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & excludebit) anyexclude = 1; + } if (!anyexclude) latte_wrapper_all(); else latte_wrapper_exclude(); @@ -316,6 +334,9 @@ void FixLatte::latte_wrapper_all() &domain->xy,&domain->xz,&domain->yz,forces,&maxiter,&latte_energy, &atom->v[0][0],&update->dt,virial,&newsystem,&latteerror); + printf("LATTE ALL: step %ld, natoms %d, LATTE eng %g\n", + update->ntimestep, natoms, latte_energy); + if (latteerror) error->all(FLERR,"Internal LATTE problem"); } @@ -325,6 +346,11 @@ void FixLatte::latte_wrapper_all() void FixLatte::latte_wrapper_exclude() { + int m; + + int exclusion_group = *exclusion_group_ptr; + int excludebit = group->bitmask[exclusion_group]; + int *mask = atom->mask; int nlocal = atom->nlocal; @@ -342,26 +368,29 @@ void FixLatte::latte_wrapper_exclude() memory->create(xinclude,nlatte,3,"latte:xinclude"); memory->create(finclude,nlatte,3,"latte:finclude"); + double *coords = &xinclude[0][0]; + int *types = typeinclude; + double *forces = &finclude[0][0]; + double **x = atom->x; int *type = atom->type; nlatte = 0; + m = 0; for (int i = 0; i < nlocal; i++) { if (mask[i] & excludebit) continue; - typeinclude[nlatte] = type[i]; - x[nlatte][0] = x[i][0]; - x[nlatte][1] = x[i][1]; - x[nlatte][2] = x[i][2]; + types[nlatte] = type[i]; nlatte++; + coords[m+0] = x[i][0]; + coords[m+1] = x[i][1]; + coords[m+2] = x[i][2]; + m += 3; } - double *coords = &xinclude[0][0]; - int *types = typeinclude; int ntypes = atom->ntypes; double *mass = &atom->mass[1]; double *boxlo = domain->boxlo; double *boxhi = domain->boxhi; - double *forces = &finclude[0][0]; bool latteerror = false; int maxiter = -1; @@ -369,13 +398,16 @@ void FixLatte::latte_wrapper_exclude() &domain->xy,&domain->xz,&domain->yz,forces,&maxiter,&latte_energy, &atom->v[0][0],&update->dt,virial,&newsystem,&latteerror); + printf("LATTE EXCLUDE: step %ld, natoms %d, LATTE eng %g\n", + update->ntimestep, nlatte, latte_energy); + if (latteerror) error->all(FLERR,"Internal LATTE problem"); // expand compressed forces array double **f = atom->f; - int m = 0; + m = 0; for (int i = 0; i < nlocal; i++) { if (mask[i] & excludebit) continue; f[i][0] = forces[m+0]; diff --git a/src/LATTE/fix_latte.h b/src/LATTE/fix_latte.h index fdc3fdee40..0ce6e0f4ce 100644 --- a/src/LATTE/fix_latte.h +++ b/src/LATTE/fix_latte.h @@ -44,10 +44,11 @@ class FixLatte : public Fix { double memory_usage() override; protected: - char *id_pe; int coulomb, pbcflag, pe_peratom, virial_global, virial_peratom, neighflag; int exclude, excludebit; int eflag_caller; + char *id_pe,*id_exclude; + int *exclusion_group_ptr; int flags_latte[6]; diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 0cc6590d5a..ef1ad324b0 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -2573,3 +2573,17 @@ void FixGCMC::grow_molecule_arrays(int nmolatoms) { molq = memory->grow(molq,nmaxmolatoms,"gcmc:molq"); molimage = memory->grow(molimage,nmaxmolatoms,"gcmc:molimage"); } + + +/* ---------------------------------------------------------------------- + extract variable which stores index of exclusion group +------------------------------------------------------------------------- */ + +void *FixGCMC::extract(const char *name, int &dim) +{ + if (strcmp(name,"exclusion_group") == 0) { + dim = 0; + return (void *) &exclusion_group; + } + return nullptr; +} diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index d945b82346..2bdd9eb461 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -31,32 +31,11 @@ class FixGCMC : public Fix { int setmask() override; void init() override; void pre_exchange() override; - void attempt_atomic_translation(); - void attempt_atomic_deletion(); - void attempt_atomic_insertion(); - void attempt_molecule_translation(); - void attempt_molecule_rotation(); - void attempt_molecule_deletion(); - void attempt_molecule_insertion(); - void attempt_atomic_translation_full(); - void attempt_atomic_deletion_full(); - void attempt_atomic_insertion_full(); - void attempt_molecule_translation_full(); - void attempt_molecule_rotation_full(); - void attempt_molecule_deletion_full(); - void attempt_molecule_insertion_full(); - double energy(int, int, tagint, double *); - double molecule_energy(tagint); - double energy_full(); - int pick_random_gas_atom(); - tagint pick_random_gas_molecule(); - void toggle_intramolecular(int); - void update_gas_atoms_list(); double compute_vector(int) override; double memory_usage() override; void write_restart(FILE *) override; void restart(char *) override; - void grow_molecule_arrays(int); + void *extract(const char *, int &); private: int molecule_group, molecule_group_bit; @@ -139,7 +118,35 @@ class FixGCMC : public Fix { class Compute *c_pe; + // private methods + void options(int, char **); + + void attempt_atomic_translation(); + void attempt_atomic_deletion(); + void attempt_atomic_insertion(); + void attempt_molecule_translation(); + void attempt_molecule_rotation(); + void attempt_molecule_deletion(); + void attempt_molecule_insertion(); + void attempt_atomic_translation_full(); + void attempt_atomic_deletion_full(); + void attempt_atomic_insertion_full(); + void attempt_molecule_translation_full(); + void attempt_molecule_rotation_full(); + void attempt_molecule_deletion_full(); + void attempt_molecule_insertion_full(); + + double energy(int, int, tagint, double *); + double energy_full(); + double molecule_energy(tagint); + + int pick_random_gas_atom(); + tagint pick_random_gas_molecule(); + void toggle_intramolecular(int); + void update_gas_atoms_list(); + + void grow_molecule_arrays(int); }; } // namespace LAMMPS_NS