diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index f12cfb7336..d7b1609619 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -15,7 +15,9 @@ General commands ================ -An alphabetic list of general LAMMPS commands. +An alphabetic list of general LAMMPS commands. Note that style +commands with many variants, can be more easily accessed via the small +table above. .. table_from_list:: :columns: 5 diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index cb61cbe17b..c5c2fb7cba 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -165,6 +165,7 @@ OPT. * :doc:`orient/fcc ` * :doc:`orient/eco ` * :doc:`pafi ` + * :doc:`pair ` * :doc:`phonon ` * :doc:`pimd ` * :doc:`planeforce ` diff --git a/doc/src/compute_pair.rst b/doc/src/compute_pair.rst index 2ec7fc0421..84ea6a803c 100644 --- a/doc/src/compute_pair.rst +++ b/doc/src/compute_pair.rst @@ -93,7 +93,9 @@ Restrictions Related commands """""""""""""""" -:doc:`compute pe `, :doc:`compute bond ` +:doc:`compute pe `, :doc:`compute bond `, +:doc:`fix pair ` + Default """"""" diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 3d79ae0101..7b16a9ce4d 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -17,7 +17,7 @@ Syntax * one or more keyword/value pairs may be appended * these keywords apply to various dump styles -* keyword = *append* or *at* or *balance* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* +* keyword = *append* or *at* or *balance* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* .. parsed-literal:: @@ -65,6 +65,8 @@ Syntax *refresh* arg = c_ID = compute ID that supports a refresh operation *scale* arg = *yes* or *no* *sfactor* arg = coordinate scaling factor (> 0.0) + *skip* arg = v_name + v_name = variable with name which evaluates to non-zero (skip) or 0 *sort* arg = *off* or *id* or N or -N off = no sorting of per-atom lines within a snapshot id = sort per-atom lines by atom ID @@ -694,14 +696,33 @@ most effective when the typical magnitude of position data is between ---------- +.. versionadded:: 15Sep2022 + +The *skip* keyword can be used with all dump styles. It allows a dump +snapshot to be skipped (not written to the dump file), if a condition +is met. The condition is computed by an :doc:`equal-style variable +`, which should be specified as v_name, where name is the +variable name. If the variable evaluation returns a non-zero value, +then the dump snapshot is skipped. If it returns zero, the dump +proceeds as usual. Note that :doc:`equal-style variable ` +can contain Boolean operators which effectively evaluate as a true +(non-zero) or false (zero) result. + +The *skip* keyword can be useful for debugging purposes, e.g. to dump +only on a particular timestep. Or to limit output to conditions of +interest, e.g. only when the force on some atom exceeds a threshold +value. + +---------- + The *sort* keyword determines whether lines of per-atom output in a snapshot are sorted or not. A sort value of *off* means they will typically be written in indeterminate order, either in serial or -parallel. This is the case even in serial if the -:doc:`atom_modify sort ` option is turned on, which it is by -default, to improve performance. A sort value of *id* means sort the output by -atom ID. A sort value of :math:`N` or :math:`-N` means sort the output by the -value in the :math:`N`\ th column of per-atom info in either ascending or +parallel. This is the case even in serial if the :doc:`atom_modify sort +` option is turned on, which it is by default, to improve +performance. A sort value of *id* means sort the output by atom ID. A +sort value of :math:`N` or :math:`-N` means sort the output by the value +in the :math:`N`\ th column of per-atom info in either ascending or descending order. The dump *local* style cannot be sorted by atom ID, since there are @@ -745,8 +766,8 @@ attributes that can be tested for are the same as those that can be specified in the :doc:`dump custom ` command, with the exception of the *element* attribute, since it is not a numeric value. Note that a different attributes can be used than those output by the -:doc:`dump custom ` command. For example, you can output the coordinates -and stress of atoms whose energy is above some threshold. +:doc:`dump custom ` command. For example, you can output the +coordinates and stress of atoms whose energy is above some threshold. If an atom-style variable is used as the attribute, then it can produce continuous numeric values or effective Boolean 0/1 values, diff --git a/doc/src/fix.rst b/doc/src/fix.rst index db08e64b14..333e920bd9 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -312,6 +312,7 @@ accelerated styles exist. * :doc:`orient/fcc ` - add grain boundary migration force for FCC * :doc:`orient/eco ` - add generalized grain boundary migration force * :doc:`pafi ` - constrained force averages on hyper-planes to compute free energies (PAFI) +* :doc:`pair ` - access per-atom info from pair styles * :doc:`phonon ` - calculate dynamical matrix from MD simulations * :doc:`pimd ` - Feynman path integral molecular dynamics * :doc:`planeforce ` - constrain atoms to move in a plane diff --git a/doc/src/fix_pair.rst b/doc/src/fix_pair.rst new file mode 100644 index 0000000000..abb44718cd --- /dev/null +++ b/doc/src/fix_pair.rst @@ -0,0 +1,110 @@ +.. index:: fix pair + +fix pair command +======================= + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID pair N pstyle name flag ... + +* ID, group-ID are documented in :doc:`fix ` command +* pair = style name of this fix command +* N = invoke this fix once every N timesteps +* pstyle = name of pair style to extract info from (e.g. eam) +* one or more name/flag pairs can be listed +* name = name of quantity the pair style allows extraction of +* flag = 1 if pair style needs to be triggered to produce data for name, 0 if not + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix request all pair 100 eam rho 0 + fix request all pair 100 amoeba uind 0 uinp 0 + + +Description +""""""""""" + +.. versionadded:: 15Sep2022 + +Extract per-atom quantities from a pair style and store them in this +fix so they can be accessed by other LAMMPS commands, e.g. by a +:doc:`dump ` command or by another :doc:`fix `, +:doc:`compute `, or :doc:`variable ` command. + +These are example use cases: + +* extract per-atom density from :doc:`pair_style eam ` to a dump file +* extract induced dipoles from :doc:`pair_style amoeba ` to a dump file +* extract accuracy metrics from a machine-learned potential to trigger output when + a condition is met (see the :doc:`dump_modify skip ` command) + +The *N* argument determines how often the fix is invoked. + +The *pstyle* argument is the name of the pair style. It can be a +sub-style used in a :doc:`pair_style hybrid ` command. + +One or more *name/flag* pairs of arguments follow. Each *name* is a +per-atom quantity which the pair style must recognize as an extraction +request. See the doc pages for individual :doc:`pair_styles +` to see what fix pair requests (if any) they support. + +The *flag* setting determines whether this fix will also trigger the +pair style to compute the named quantity so it can be extracted. If the +quantity is always computed by the pair style, no trigger is needed; +specify *flag* = 0. If the quantity is not always computed +(e.g. because it is expensive to calculate), then specify *flag* = 1. +This will trigger the quantity to be calculated only on timesteps it is +needed. Again, see the doc pages for individual :doc:`pair_styles +` to determine which fix pair requests (if any) need to be +triggered with a *flag* = 1 setting. + +The per-atom data extracted from the pair style is stored by this fix +as either a per-atom vector or array. If there is only one *name* +argument specified and the pair style computes a single value for each +atom, then this fix stores it as a per-atom vector. Otherwise a +per-atom array is created, with its data in the order of the *name* +arguments. + +For example, :doc:`pair_style amoeba ` allows extraction of +two named quantities: "uind" and "uinp", both of which are 3-vectors for +each atom, i.e. dipole moments. In the example below a 6-column per-atom +array will be created. Columns 1-3 will store the "uind" values; +columns 4-6 will store the "uinp" values. + +.. code-block:: LAMMPS + + pair_style amoeba + fix ex all pair amoeba 10 uind 0 uinp 0 + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. + +As explained above, this fix produces a per-atom vector or array which +can be accessed by various :doc:`output commands `. If +an array is produced, the number of columns is the sum of the number +of per-atom quantities produced by each *name* argument requested from +the pair style. + +Restrictions +"""""""""""" +none + +Related commands +"""""""""""""""" + +:doc:`compute pair ` + +Default +""""""" + +none diff --git a/doc/src/pair_amoeba.rst b/doc/src/pair_amoeba.rst index a58eb2fabc..44a14f0f45 100644 --- a/doc/src/pair_amoeba.rst +++ b/doc/src/pair_amoeba.rst @@ -156,6 +156,20 @@ settings. ---------- +.. versionadded:: TBD + +The *amoeba* and *hippo* pair styles support extraction of two per-atom +quantities by the :doc:`fix pair ` command. This allows the +quantities to be output to files by the :doc:`dump ` or otherwise +processed by other LAMMPS commands. + +The names of the two quantities are "uind" and "uinp" for the induced +dipole moments for each atom. Neither quantity needs to be triggered by +the :doc:`fix pair ` command in order for these pair styles to +calculate it. + +---------- + Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/pair_eam.rst b/doc/src/pair_eam.rst index fd05d7189e..af936a6bc9 100644 --- a/doc/src/pair_eam.rst +++ b/doc/src/pair_eam.rst @@ -444,6 +444,20 @@ identical to the FS EAM files (see above). ---------- +.. versionadded:: TBD + +The *eam*, *eam/alloy*, *eam/fs*, and *eam/he* pair styles support +extraction of two per-atom quantities by the :doc:`fix pair ` +command. This allows the quantities to be output to files by the +:doc:`dump ` or otherwise processed by other LAMMPS commands. + +The names of the two quantities are "rho" and "fp" for the density and +derivative of the embedding energy for each atom. Neither quantity +needs to be triggered by the :doc:`fix pair ` command in order +for these pair styles to calculate it. + +---------- + .. include:: accel_styles.rst ---------- @@ -459,21 +473,26 @@ a pair_coeff command with I != J arguments for the eam styles. This pair style does not support the :doc:`pair_modify ` shift, table, and tail options. -The eam pair styles do not write their information to :doc:`binary restart files `, since it is stored in tabulated potential files. -Thus, you need to re-specify the pair_style and pair_coeff commands in -an input script that reads a restart file. +The eam pair styles do not write their information to :doc:`binary +restart files `, since it is stored in tabulated potential +files. Thus, you need to re-specify the pair_style and pair_coeff +commands in an input script that reads a restart file. The eam pair styles can only be used via the *pair* keyword of the :doc:`run_style respa ` command. They do not support the *inner*, *middle*, *outer* keywords. + + + ---------- Restrictions """""""""""" All of these styles are part of the MANYBODY package. They are only -enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +enabled if LAMMPS was built with that package. See the :doc:`Build +package ` page for more info. Related commands """""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c6c6a2de6a..15316b4f09 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -3581,7 +3581,9 @@ UEF ufm Uhlenbeck Ui +uind uInfParallel +uinp uk ul ulb diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 1f1baa0280..a717416524 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -2094,6 +2094,28 @@ void *PairAmoeba::extract(const char *str, int &dim) return nullptr; } +/* ---------------------------------------------------------------------- + peratom requests from FixPair + return ptr to requested data + also return ncol = # of quantites per atom + 0 = per-atom vector + 1 or more = # of columns in per-atom array + return NULL if str is not recognized +---------------------------------------------------------------------- */ + +void *PairAmoeba::extract_peratom(const char *str, int &ncol) +{ + if (strcmp(str,"uind") == 0) { + ncol = 3; + return (void *) uind; + } else if (strcmp(str,"uinp") == 0) { + ncol = 3; + return (void *) uinp; + } + + return nullptr; +} + /* ---------------------------------------------------------------------- grow local vectors and arrays if necessary keep them all atom->nmax in length even if ghost storage not needed diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index ad80a74879..7f6596c3bd 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -50,6 +50,7 @@ class PairAmoeba : public Pair { void unpack_reverse_grid(int, void *, int, int *) override; void *extract(const char *, int &) override; + void *extract_peratom(const char *, int &) override; double memory_usage() override; protected: diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 834c71a947..b14d47f7f0 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -933,5 +933,28 @@ void *PairEAM::extract(const char *str, int &dim) { dim = 2; if (strcmp(str,"scale") == 0) return (void *) scale; + + return nullptr; +} + +/* ---------------------------------------------------------------------- + peratom requests from FixPair + return ptr to requested data + also return ncol = # of quantites per atom + 0 = per-atom vector + 1 or more = # of columns in per-atom array + return NULL if str is not recognized +---------------------------------------------------------------------- */ + +void *PairEAM::extract_peratom(const char *str, int &ncol) +{ + if (strcmp(str,"rho") == 0) { + ncol = 0; + return (void *) rho; + } else if (strcmp(str,"fp") == 0) { + ncol = 0; + return (void *) fp; + } + return nullptr; } diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h index 3589ab4ab0..11ab969d18 100644 --- a/src/MANYBODY/pair_eam.h +++ b/src/MANYBODY/pair_eam.h @@ -53,6 +53,7 @@ class PairEAM : public Pair { double init_one(int, int) override; double single(int, int, int, int, double, double, double, double &) override; void *extract(const char *, int &) override; + void *extract_peratom(const char *, int &) override; int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; diff --git a/src/dump.cpp b/src/dump.cpp index e8134adc0f..42d7186bcc 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -20,11 +20,13 @@ #include "error.h" #include "fix.h" #include "group.h" +#include "input.h" #include "irregular.h" #include "memory.h" #include "modify.h" #include "output.h" #include "update.h" +#include "variable.h" #include @@ -87,6 +89,9 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp) delay_flag = 0; write_header_flag = 1; + skipflag = 0; + skipvar = nullptr; + maxfiles = -1; numfiles = 0; fileidx = 0; @@ -164,6 +169,7 @@ Dump::~Dump() delete[] format_bigint_user; delete[] refresh; + delete[] skipvar; // format_column_user is deallocated by child classes that use it @@ -298,6 +304,15 @@ void Dump::init() else error->all(FLERR,"Dump could not find refresh compute ID"); } + // if skipflag, check skip variable + + if (skipflag) { + skipindex = input->variable->find(skipvar); + if (skipindex < 0) error->all(FLERR,"Dump skip variable not found"); + if (!input->variable->equalstyle(skipindex)) + error->all(FLERR,"Variable for dump skip is invalid style"); + } + // preallocation for PBC copies if requested if (pbcflag && atom->nlocal > maxpbc) pbc_allocate(); @@ -325,15 +340,6 @@ void Dump::write() imageint *imagehold; double **xhold,**vhold; - // if timestep < delaystep, just return - - if (delay_flag && update->ntimestep < delaystep) return; - - // if file per timestep, open new file - - if (multifile) openfile(); - if (fp) clearerr(fp); - // simulation box bounds if (domain->triclinic == 0) { @@ -359,6 +365,24 @@ void Dump::write() nme = count(); + // if timestep < delaystep, just return + // if skip condition is defined and met, just return + // must do both these tests after count() b/c it invokes computes, + // this enables caller to trigger future invocation of needed computes + + if (delay_flag && update->ntimestep < delaystep) return; + + if (skipflag) { + double value = input->variable->compute_equal(skipindex); + if (value != 0.0) return; + } + + // if file per timestep, open new file + // do this after skip check, so no file is opened if skip occurs + + if (multifile) openfile(); + if (fp) clearerr(fp); + // ntotal = total # of dump lines in snapshot // nmax = max # of dump lines on any proc @@ -1264,6 +1288,15 @@ void Dump::modify_params(int narg, char **arg) pbcflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; + } else if (strcmp(arg[iarg],"skip") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + skipflag = 1; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + delete[] skipvar; + skipvar = utils::strdup(&arg[iarg+1][2]); + } else error->all(FLERR,"Illegal dump_modify command"); + iarg += 2; + } else if (strcmp(arg[iarg],"sort") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "dump_modify sort", error); if (strcmp(arg[iarg+1],"off") == 0) sort_flag = 0; diff --git a/src/dump.h b/src/dump.h index 7d17288710..ce2025fe79 100644 --- a/src/dump.h +++ b/src/dump.h @@ -97,6 +97,10 @@ class Dump : protected Pointers { char *refresh; // compute ID to invoke refresh() on int irefresh; // index of compute + int skipflag; // 1 if skip condition defined + char *skipvar; // name of variable to check for skip condition + int skipindex; // index of skip variable + char boundstr[9]; // encoding of boundary flags char *format; // format string for the file write diff --git a/src/fix_pair.cpp b/src/fix_pair.cpp new file mode 100644 index 0000000000..6f98ad2790 --- /dev/null +++ b/src/fix_pair.cpp @@ -0,0 +1,343 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_pair.h" + +#include "atom.h" +#include "dump.h" +#include "error.h" +#include "force.h" +#include "fix.h" +#include "input.h" +#include "memory.h" +#include "pair.h" +#include "output.h" +#include "variable.h" +#include "update.h" +#include "variable.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 7) utils::missing_cmd_args(FLERR, "fix pair", error); + + nevery = utils::inumeric(FLERR,arg[3],false,lmp); + if (nevery < 1) error->all(FLERR,"Illegal fix pair every value: {}", nevery); + + pairname = utils::strdup(arg[4]); + pstyle = force->pair_match(pairname,1,0); + if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname); + + nfield = (narg-5) / 2; + fieldname = new char*[nfield]; + trigger = new int[nfield]; + + nfield = 0; + int iarg = 5; + + while (iarg < narg) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix pair {}", arg[iarg]), error); + fieldname[nfield] = utils::strdup(arg[iarg]); + int flag = utils::inumeric(FLERR,arg[iarg+1],true,lmp); + if (flag == 0) trigger[nfield] = 0; + else if (flag == 1) trigger[nfield] = 1; + else error->all(FLERR,"Illegal fix pair {} command flag: {}", arg[iarg], arg[iarg+1]); + nfield++; + iarg += 2; + } + + // set trigger names = fieldname + "_flag" + + triggername = new char*[nfield]; + + for (int ifield = 0; ifield < nfield; ifield++) { + if (trigger[ifield] == 0) triggername[ifield] = nullptr; + else triggername[nfield] = utils::strdup(fmt::format("{}_flag", fieldname[ifield])); + } + + // extract all fields just to get number of per-atom values + // returned data ptr may be NULL, if pair style has not allocated field yet + // check for recognized field cannot be done until post_force() + // also check if triggername can be extracted as a scalar value + + triggerptr = new int*[nfield]; + + ncols = 0; + + for (int ifield = 0; ifield < nfield; ifield++) { + int columns = 0; // set in case fieldname not recognized by pstyle + void *pvoid = pstyle->extract_peratom(fieldname[ifield],columns); + if (columns) ncols += columns; + else ncols++; + if (trigger[ifield]) { + int dim; + triggerptr[ifield] = (int *) pstyle->extract(triggername[ifield],dim); + if (!triggerptr[ifield]) + error->all(FLERR,"Fix pair pair style cannot extract {}", triggername[ifield]); + if (dim) + error->all(FLERR,"Fix pair pair style {} trigger {} is not a scalar", + pairname, triggername[ifield]); + } + } + + // if set peratom_freq = Nevery, then cannot access the per-atom + // values as part of thermo output during minimiziation + // at different frequency or on last step of minimization + // instead set peratom_freq = 1 + // ok, since vector/array always have values + // but requires the vector/array be persisted between Nevery steps + // since it may be accessed + + peratom_flag = 1; + if (ncols == 1) size_peratom_cols = 0; + else size_peratom_cols = ncols; + peratom_freq = 1; + + // perform initial allocation of atom-based array + // register with Atom class + + vector = nullptr; + array = nullptr; + grow_arrays(atom->nmax); + atom->add_callback(Atom::GROW); + + // zero the vector/array since dump may access it on timestep 0 + // zero the vector/array since a variable may access it before first run + + int nlocal = atom->nlocal; + + if (ncols == 1) { + for (int i = 0; i < nlocal; i++) + vector[i] = 0.0; + } else { + for (int i = 0; i < nlocal; i++) + for (int m = 0; m < ncols; m++) + array[i][m] = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +FixPair::~FixPair() +{ + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,Atom::GROW); + + delete[] pairname; + for (int ifield = 0; ifield < nfield; ifield++) { + delete[] fieldname[ifield]; + delete[] triggername[ifield]; + } + + delete[] fieldname; + delete[] trigger; + delete[] triggername; + delete[] triggerptr; + + if (ncols == 1) memory->destroy(vector); + else memory->destroy(array); +} + +/* ---------------------------------------------------------------------- */ + +int FixPair::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + mask |= MIN_PRE_FORCE; + mask |= POST_FORCE; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixPair::init() +{ + // insure pair style still exists + + pstyle = force->pair_match(pairname,1,0); + if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname); +} + +/* ---------------------------------------------------------------------- */ + +void FixPair::setup(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixPair::setup_pre_force(int vflag) +{ + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- + trigger pair style computation on steps which are multiples of Nevery +------------------------------------------------------------------------- */ + +void FixPair::pre_force(int /*vflag*/) +{ + if (update->ntimestep % nevery) return; + + // set pair style triggers + + for (int ifield = 0; ifield < nfield; ifield++) + if (trigger[ifield]) *(triggerptr[ifield]) = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixPair::min_pre_force(int vflag) +{ + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- + extract results from pair style +------------------------------------------------------------------------- */ + +void FixPair::post_force(int /*vflag*/) +{ + if (update->ntimestep % nevery) return; + + // extract pair style fields one by one + // store their values in this fix + + int nlocal = atom->nlocal; + + int icol = 0; + int columns; + + for (int ifield = 0; ifield < nfield; ifield++) { + void *pvoid = pstyle->extract_peratom(fieldname[ifield],columns); + if (pvoid == nullptr) + error->all(FLERR,"Fix pair pair style cannot extract {}",fieldname[ifield]); + + if (columns == 0) { + double *pvector = (double *) pvoid; + if (ncols == 1) { + for (int i = 0; i < nlocal; i++) + vector[i] = pvector[i]; + } else { + for (int i = 0; i < nlocal; i++) + array[i][icol] = pvector[i]; + } + icol++; + + } else { + double **parray = (double **) pvoid; + for (int i = 0; i < nlocal; i++) + for (int m = 0; m < columns; m++) { + array[i][icol] = parray[i][m]; + icol++; + } + } + } + + // unset pair style triggers + + for (int ifield = 0; ifield < nfield; ifield++) + if (trigger[ifield]) *(triggerptr[ifield]) = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixPair::min_post_force(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- + allocate atom-based vector or array +------------------------------------------------------------------------- */ + +void FixPair::grow_arrays(int nmax) +{ + if (ncols == 1) { + memory->grow(vector,nmax,"store/state:vector"); + vector_atom = vector; + } else { + memory->grow(array,nmax,ncols,"store/state:array"); + array_atom = array; + } +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixPair::copy_arrays(int i, int j, int /*delflag*/) +{ + if (ncols == 1) { + vector[j] = vector[i]; + } else { + for (int m = 0; m < ncols; m++) + array[j][m] = array[i][m]; + } +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixPair::pack_exchange(int i, double *buf) +{ + if (ncols == 1) { + buf[0] = vector[i]; + } else { + for (int m = 0; m < ncols; m++) + buf[m] = array[i][m]; + } + + return ncols; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixPair::unpack_exchange(int nlocal, double *buf) +{ + if (ncols == 1) { + vector[nlocal] = buf[0]; + } else { + for (int m = 0; m < ncols; m++) + array[nlocal][m] = buf[m]; + } + + return ncols; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based vector or array +------------------------------------------------------------------------- */ + +double FixPair::memory_usage() +{ + double bytes = 0.0; + if (ncols == 1) bytes += (double)atom->nmax * sizeof(double); + else bytes += (double)atom->nmax*ncols * sizeof(double); + return bytes; +} diff --git a/src/fix_pair.h b/src/fix_pair.h new file mode 100644 index 0000000000..c4d0bb593c --- /dev/null +++ b/src/fix_pair.h @@ -0,0 +1,62 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(pair,FixPair); +// clang-format on +#else + +#ifndef LMP_FIX_PAIR_H +#define LMP_FIX_PAIR_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixPair : public Fix { + public: + FixPair(class LAMMPS *, int, char **); + ~FixPair() override; + int setmask() override; + void init() override; + void setup(int) override; + void setup_pre_force(int) override; + void pre_force(int) override; + void min_pre_force(int) override; + void post_force(int) override; + void min_post_force(int) override; + + void grow_arrays(int) override; + void copy_arrays(int, int, int) override; + int pack_exchange(int, double *) override; + int unpack_exchange(int, double *) override; + + double memory_usage() override; + + private: + int nevery,nfield,ncols; + char *pairname; + char **fieldname,**triggername; + int *trigger; + int **triggerptr; + + class Pair *pstyle; + double *vector; + double **array; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/pair.h b/src/pair.h index 2a68aea5fa..4f4d66ed97 100644 --- a/src/pair.h +++ b/src/pair.h @@ -215,6 +215,7 @@ class Pair : protected Pointers { // specific child-class methods for certain Pair styles virtual void *extract(const char *, int &) { return nullptr; } + virtual void *extract_peratom(const char *, int &) { return nullptr; } virtual void swap_eam(double *, double **) {} virtual void reset_dt() {} virtual void min_xf_pointers(int, double **, double **) {}