From 3eb8fd219afbb97303efd4fcc5dcd85958dd7954 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 16 Aug 2023 16:19:11 -0600 Subject: [PATCH 01/22] update of a Howto_output doc page --- doc/src/Howto_output.rst | 142 ++++++++++++++++++++++++++------------- 1 file changed, 94 insertions(+), 48 deletions(-) diff --git a/doc/src/Howto_output.rst b/doc/src/Howto_output.rst index 851b7703fd4..6fcd36ab56c 100644 --- a/doc/src/Howto_output.rst +++ b/doc/src/Howto_output.rst @@ -1,7 +1,7 @@ Output from LAMMPS (thermo, dumps, computes, fixes, variables) ============================================================== -There are four basic kinds of LAMMPS output: +There are four basic forms of LAMMPS output: * :doc:`Thermodynamic output `, which is a list of quantities printed every few timesteps to the screen and logfile. @@ -20,18 +20,17 @@ output files, depending on what :doc:`dump ` and :doc:`fix ` commands you specify. As discussed below, LAMMPS gives you a variety of ways to determine -what quantities are computed and printed when the thermodynamics, +what quantities are calculated and printed when the thermodynamics, dump, or fix commands listed above perform output. Throughout this discussion, note that users can also :doc:`add their own computes and -fixes to LAMMPS ` which can then generate values that can then -be output with these commands. +fixes to LAMMPS ` which can generate values that can then be +output with these commands. The following subsections discuss different LAMMPS commands related to output and the kind of data they operate on and produce: * :ref:`Global/per-atom/local/per-grid data ` * :ref:`Scalar/vector/array data ` -* :ref:`Per-grid data ` * :ref:`Disambiguation ` * :ref:`Thermodynamic output ` * :ref:`Dump file output ` @@ -48,34 +47,65 @@ to output and the kind of data they operate on and produce: Global/per-atom/local/per-grid data ----------------------------------- -Various output-related commands work with four different styles of +Various output-related commands work with four different "styles" of data: global, per-atom, local, and per-grid. A global datum is one or more system-wide values, e.g. the temperature of the system. A per-atom datum is one or more values per atom, e.g. the kinetic energy of each atom. Local datums are calculated by each processor based on -the atoms it owns, but there may be zero or more per atom, e.g. a list +the atoms it owns, and there may be zero or more per atom, e.g. a list of bond distances. A per-grid datum is one or more values per grid cell, for a grid which -overlays the simulation domain. The grid cells and the data they -store are distributed across processors; each processor owns the grid -cells whose center point falls within its subdomain. +overlays the simulation domain. Similar to atoms and per-atom data, +the grid cells and the data they store are distributed across +processors; each processor owns the grid cells whose center points +fall within its subdomain. .. _scalar: Scalar/vector/array data ------------------------ -Global, per-atom, and local datums can come in three kinds: a single -scalar value, a vector of values, or a 2d array of values. The doc -page for a "compute" or "fix" or "variable" that generates data will -specify both the style and kind of data it produces, e.g. a per-atom -vector. - -When a quantity is accessed, as in many of the output commands -discussed below, it can be referenced via the following bracket -notation, where ID in this case is the ID of a compute. The leading -"c\_" would be replaced by "f\_" for a fix, or "v\_" for a variable: +Global, per-atom, local, and per-grid datums can come in three +"kinds": a single scalar value, a vector of values, or a 2d array of +values. More specifically these are the valid kinds for each style: + +* global scalar +* global vector +* global array +* per-atom vector +* per-atom array +* local vector +* local array +* per-grid vector +* per-grid array + +A per-atom vector means a single value per atom; the "vector" is the +length of the number of atoms. A per-atom array means multiple values +per atom. Similarly a local vector or array means one or multiple +values per entity (e.g. per bond in the system). And a per-grid +vector or array means one or multiple values per grid cell. + +The doc page for a compute or fix or variable that generates data will +specify both the styles and kinds of data it produces, e.g. a per-atom +vector. Note that a compute or fix may generate multiple styles and +kinds of output. However, for per-atom data only a vector or array is +output, never both. Likewise for per-local and per-grid data. An +example of a fix which generates multiple styles and kinds of data is +the :doc:`fix mdi/qm ` command. It outputs a global +scalar, global vector, and per-atom array for the quantum mechanical +energy and virial of the system and forces on each atom. + +By contrast, different variable styles generate only a single kind of +data: a global scalar for an equal-style variable, global vector for a +vector-style variable, and a per-atom vector for an atom-style +variable. + +When data is accessed by another command, as in many of the output +commands discussed below, it can be referenced via the following +bracket notation, where ID in this case is the ID of a compute. The +leading "c\_" would be replaced by "f\_" for a fix, or "v\_" for a +variable (and ID would be the name of the variable): +-------------+--------------------------------------------+ | c_ID | entire scalar, vector, or array | @@ -85,40 +115,56 @@ notation, where ID in this case is the ID of a compute. The leading | c_ID[I][J] | one element of array | +-------------+--------------------------------------------+ -In other words, using one bracket reduces the dimension of the data -once (vector -> scalar, array -> vector). Using two brackets reduces -the dimension twice (array -> scalar). Thus a command that uses -scalar values as input can typically also process elements of a vector -or array. - -.. _grid: - -Per-grid data ------------------------- +Note that using one bracket reduces the dimension of the data once +(vector -> scalar, array -> vector). Using two brackets reduces the +dimension twice (array -> scalar). Thus a command that uses scalar +values as input can also conceptually operate on an element of a +vector or array. -Per-grid data can come in two kinds: a vector of values (one per grid -cekk), or a 2d array of values (multiple values per grid ckk). The -doc page for a "compute" or "fix" that generates data will specify -names for both the grid(s) and datum(s) it produces, e.g. per-grid -vectors or arrays, which can be referenced by other commands. See the -:doc:`Howto grid ` doc page for more details. +Per-grid vectors or arrays are accessed similarly, except that the ID +for the compute or fix includes a grid name and a data name. This is +because a fix or compute can create multiple grids (of different +sizes) and multiple sets of data (for each grid). The fix or compute +defines names for each grid and for each data set, so that all of them +can be accessed by other commands. See the :doc:`Howto grid +` doc page for more details. .. _disambiguation: Disambiguation -------------- -Some computes and fixes produce data in multiple styles, e.g. a global -scalar and a per-atom vector. Usually the context in which the input -script references the data determines which style is meant. Example: -if a compute provides both a global scalar and a per-atom vector, the -former will be accessed by using ``c_ID`` in an equal-style variable, -while the latter will be accessed by using ``c_ID`` in an atom-style -variable. Note that atom-style variable formulas can also access -global scalars, but in this case it is not possible to do this -directly because of the ambiguity. Instead, an equal-style variable -can be defined which accesses the global scalar, and that variable can -be used in the atom-style variable formula in place of ``c_ID``. +When a compute or fix produces data in multiple styles, e.g. global +and per-atom, a reference to the data can sometimes be ambiguous. +Usually the context in which the input script references the data +determines which style is meant. + +For example, if a compute outputs a global vector and a per-atom +array, an element of the global vector will be accessed by using +``c_ID[I]`` in :doc:`thermodynamic output `, while a +column of the per-atom array will be accessed by using ``c_ID[I]`` in +a :doc:`dump custom ` command. + +However, if a :doc:`atom-style variable ` references +``c_ID[I]``, then it could be intended to refer to a single element of +the global vector or a column of the per-atom array. The doc page for +any command that has a potential ambiguity (variables are the most +common) will explain how to resolve the ambiguity. + +In this case, an atom-style variables references per-atom data if it +exists. If access to an element of a global vector is needed (as in +this example), an equal-style variable which references the value can +be defined and used in the atom-style variable formula instead. + +Similarly, :doc:`thermodynamic output ` can only +reference global data from a compute or fix. But you can indirectly +access per-atom data as follows. The reference ``c_ID[245][2]`` for +the ID of a :doc:`compute displace/atom ` +command, refers to the y-component of displacement for the atom with +ID 245. While you cannot use that reference directly in the +:doc:`thermo_style ` command, you can use it an +equal-style variable formula, and then reference the variable in +thermodynamic output. .. _thermo: @@ -389,7 +435,7 @@ output and input data types must match, e.g. global/per-atom/local data and scalar/vector/array data. Also note that, as described above, when a command takes a scalar as -input, that could be an element of a vector or array. Likewise a +input, that could also be an element of a vector or array. Likewise a vector input could be a column of an array. +--------------------------------------------------------+----------------------------------------------+----------------------------------------------------+ From 0d739439c7cebf64716561022151921c04cd35c8 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 17 Aug 2023 12:47:48 -0600 Subject: [PATCH 02/22] changes to compute voronoi/atom --- doc/src/compute_voronoi_atom.rst | 126 ++++++++++++--------------- src/VORONOI/compute_voronoi_atom.cpp | 37 ++++---- 2 files changed, 75 insertions(+), 88 deletions(-) diff --git a/doc/src/compute_voronoi_atom.rst b/doc/src/compute_voronoi_atom.rst index 274be1b702b..3e67bb6cbfa 100644 --- a/doc/src/compute_voronoi_atom.rst +++ b/doc/src/compute_voronoi_atom.rst @@ -13,7 +13,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * voronoi/atom = style name of this compute command * zero or more keyword/value pairs may be appended -* keyword = *only_group* or *occupation* or *surface* or *radius* or *edge_histo* or *edge_threshold* or *face_threshold* or *neighbors* or *peratom* +* keyword = *only_group* or *occupation* or *surface* or *radius* or *edge_histo* or *edge_threshold* or *face_threshold* or *neighbors* .. parsed-literal:: @@ -31,7 +31,6 @@ Syntax *face_threshold* arg = minarea minarea = minimum area for a face to be counted *neighbors* value = *yes* or *no* = store list of all neighbors or no - *peratom* value = *yes* or *no* = per-atom quantities accessible or no Examples """""""" @@ -53,14 +52,12 @@ atoms in the simulation box. The tessellation is calculated using all atoms in the simulation, but non-zero values are only stored for atoms in the group. -By default two per-atom quantities are calculated by this compute. -The first is the volume of the Voronoi cell around each atom. Any -point in an atom's Voronoi cell is closer to that atom than any other. -The second is the number of faces of the Voronoi cell. This is -equal to the number of nearest neighbors of the central atom, -plus any exterior faces (see note below). If the *peratom* keyword -is set to "no", the per-atom quantities are still calculated, -but they are not accessible. +Two per-atom quantities are calculated by this compute. The first is +the volume of the Voronoi cell around each atom. Any point in an +atom's Voronoi cell is closer to that atom than any other. The second +is the number of faces of the Voronoi cell. This is equal to the +number of nearest neighbors of the central atom, plus any exterior +faces (see note below). ---------- @@ -97,13 +94,13 @@ present in atom_style sphere for granular models. The *edge_histo* keyword activates the compilation of a histogram of number of edges on the faces of the Voronoi cells in the compute -group. The argument *maxedge* of the this keyword is the largest number -of edges on a single Voronoi cell face expected to occur in the -sample. This keyword adds the generation of a global vector with -*maxedge*\ +1 entries. The last entry in the vector contains the number of -faces with more than *maxedge* edges. Since the polygon with the -smallest amount of edges is a triangle, entries 1 and 2 of the vector -will always be zero. +group. The argument *maxedge* of the this keyword is the largest +number of edges on a single Voronoi cell face expected to occur in the +sample. This keyword generates output of a global vector by this +compute with *maxedge*\ +1 entries. The last entry in the vector +contains the number of faces with more than *maxedge* edges. Since the +polygon with the smallest amount of edges is a triangle, entries 1 and +2 of the vector will always be zero. The *edge_threshold* and *face_threshold* keywords allow the suppression of edges below a given minimum length and faces below a @@ -127,8 +124,8 @@ to locate vacancies (the coordinates are given by the atom coordinates at the time step when the compute was first invoked), while column two data can be used to identify interstitial atoms. -If the *neighbors* value is set to yes, then this compute creates a -local array with 3 columns. There is one row for each face of each +If the *neighbors* value is set to yes, then this compute also creates +a local array with 3 columns. There is one row for each face of each Voronoi cell. The 3 columns are the atom ID of the atom that owns the cell, the atom ID of the atom in the neighboring cell (or zero if the face is external), and the area of the face. The array can be @@ -143,8 +140,8 @@ containing all the Voronoi neighbors in a system: compute 6 all voronoi/atom neighbors yes dump d2 all local 1 dump.neighbors index c_6[1] c_6[2] c_6[3] -If the *face_threshold* keyword is used, then only faces -with areas greater than the threshold are stored. +If the *face_threshold* keyword is used, then only faces with areas +greater than the threshold are stored. ---------- @@ -158,48 +155,48 @@ Voro++ software in the src/VORONOI/README file. .. note:: - The calculation of Voronoi volumes is performed by each processor for - the atoms it owns, and includes the effect of ghost atoms stored by - the processor. This assumes that the Voronoi cells of owned atoms - are not affected by atoms beyond the ghost atom cut-off distance. - This is usually a good assumption for liquid and solid systems, but - may lead to underestimation of Voronoi volumes in low density - systems. By default, the set of ghost atoms stored by each processor - is determined by the cutoff used for :doc:`pair_style ` - interactions. The cutoff can be set explicitly via the - :doc:`comm_modify cutoff ` command. The Voronoi cells - for atoms adjacent to empty regions will extend into those regions up - to the communication cutoff in :math:`x`, :math:`y`, or :math:`z`. - In that situation, an exterior face is created at the cutoff distance - normal to the :math:`x`, :math:`y`, or :math:`z` direction. For - triclinic systems, the exterior face is parallel to the corresponding - reciprocal lattice vector. + The calculation of Voronoi volumes is performed by each processor + for the atoms it owns, and includes the effect of ghost atoms + stored by the processor. This assumes that the Voronoi cells of + owned atoms are not affected by atoms beyond the ghost atom cut-off + distance. This is usually a good assumption for liquid and solid + systems, but may lead to underestimation of Voronoi volumes in low + density systems. By default, the set of ghost atoms stored by each + processor is determined by the cutoff used for :doc:`pair_style + ` interactions. The cutoff can be set explicitly via + the :doc:`comm_modify cutoff ` command. The Voronoi + cells for atoms adjacent to empty regions will extend into those + regions up to the communication cutoff in :math:`x`, :math:`y`, or + :math:`z`. In that situation, an exterior face is created at the + cutoff distance normal to the :math:`x`, :math:`y`, or :math:`z` + direction. For triclinic systems, the exterior face is parallel to + the corresponding reciprocal lattice vector. .. note:: - The Voro++ package performs its calculation in 3d. This will - still work for a 2d LAMMPS simulation, provided all the atoms have the - same :math:`z`-coordinate. The Voronoi cell of each atom will be a columnar - polyhedron with constant cross-sectional area along the :math:`z`-direction - and two exterior faces at the top and bottom of the simulation box. If - the atoms do not all have the same :math:`z`-coordinate, then the columnar - cells will be accordingly distorted. The cross-sectional area of each - Voronoi cell can be obtained by dividing its volume by the :math:`z` extent - of the simulation box. Note that you define the :math:`z` extent of the - simulation box for 2d simulations when using the - :doc:`create_box ` or :doc:`read_data ` commands. + The Voro++ package performs its calculation in 3d. This will still + work for a 2d LAMMPS simulation, provided all the atoms have the + same :math:`z`-coordinate. The Voronoi cell of each atom will be a + columnar polyhedron with constant cross-sectional area along the + :math:`z`-direction and two exterior faces at the top and bottom of + the simulation box. If the atoms do not all have the same + :math:`z`-coordinate, then the columnar cells will be accordingly + distorted. The cross-sectional area of each Voronoi cell can be + obtained by dividing its volume by the :math:`z` extent of the + simulation box. Note that you define the :math:`z` extent of the + simulation box for 2d simulations when using the :doc:`create_box + ` or :doc:`read_data ` commands. Output info """"""""""" -By default, this compute calculates a per-atom array with two -columns. In regular dynamic tessellation mode the first column is the -Voronoi volume, the second is the neighbor count, as described above -(read above for the output data in case the *occupation* keyword is -specified). These values can be accessed by any command that uses -per-atom values from a compute as input. See the :doc:`Howto output ` page for an overview of LAMMPS output -options. If the *peratom* keyword is set to "no", the per-atom array -is still created, but it is not accessible. +This compute calculates a per-atom array with two columns. In regular +dynamic tessellation mode the first column is the Voronoi volume, the +second is the neighbor count, as described above (read above for the +output data in case the *occupation* keyword is specified). These +values can be accessed by any command that uses per-atom values from a +compute as input. See the :doc:`Howto output ` page for +an overview of LAMMPS output options. If the *edge_histo* keyword is used, then this compute generates a global vector of length *maxedge*\ +1, containing a histogram of the @@ -209,17 +206,6 @@ If the *neighbors* value is set to *yes*, then this compute calculates a local array with three columns. There is one row for each face of each Voronoi cell. -.. note:: - - Some LAMMPS commands such as the :doc:`compute reduce ` - command can accept either a per-atom or local quantity. If this compute - produces both quantities, the command - may access the per-atom quantity, even if you want to access the local - quantity. This effect can be eliminated by using the *peratom* - keyword to turn off the production of the per-atom quantities. For - the default value *yes* both quantities are produced. For the value - *no*, only the local array is produced. - The Voronoi cell volume will be in distance :doc:`units ` cubed. The Voronoi face area will be in distance :doc:`units ` squared. @@ -227,7 +213,8 @@ Restrictions """""""""""" This compute is part of the VORONOI package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. It also requires you have a copy of the Voro++ library built and installed on your system. See instructions on obtaining and @@ -241,5 +228,4 @@ Related commands Default """"""" -*neighbors* no, *peratom* yes - +The default for the neighobrs keyword is no. diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp index 28bab271a28..eb4f53986f8 100644 --- a/src/VORONOI/compute_voronoi_atom.cpp +++ b/src/VORONOI/compute_voronoi_atom.cpp @@ -111,12 +111,7 @@ ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) : if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command"); faces_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; - } else if (strcmp(arg[iarg], "peratom") == 0) { - if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command"); - peratom_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); - iarg += 2; - } - else error->all(FLERR,"Illegal compute voronoi/atom command"); + } else error->all(FLERR,"Illegal compute voronoi/atom command"); } if (occupation && ( surface!=VOROSURF_NONE || maxedge>0 ) ) @@ -394,27 +389,29 @@ void ComputeVoronoi::checkOccupation() // clear occupation vector memset(occvec, 0, oldnatoms*sizeof(*occvec)); - int i, j, k, - nlocal = atom->nlocal, - nall = atom->nghost + nlocal; - double rx, ry, rz, - **x = atom->x; + int i, j, k; + double rx, ry, rz; + + int nlocal = atom->nlocal; + int nall = atom->nghost + nlocal; + double **x = atom->x; // prepare destination buffer for variable evaluation + if (atom->nmax > lmax) { memory->destroy(lnext); lmax = atom->nmax; memory->create(lnext,lmax,"voronoi/atom:lnext"); } - // clear lroot - for (i=0; ifind_voronoi_cell(x[i][0], x[i][1], x[i][2], rx, ry, rz, k)) || @@ -435,6 +432,7 @@ void ComputeVoronoi::checkOccupation() } // MPI sum occupation + #ifdef NOTINPLACE memcpy(sendocc, occvec, oldnatoms*sizeof(*occvec)); MPI_Allreduce(sendocc, occvec, oldnatoms, MPI_INT, MPI_SUM, world); @@ -443,6 +441,7 @@ void ComputeVoronoi::checkOccupation() #endif // determine the total number of atoms in this atom's currently occupied cell + int c; for (i=0; itag[i]; if (mytag > oldmaxtag) voro[i][0] = 0; @@ -479,6 +479,7 @@ void ComputeVoronoi::checkOccupation() void ComputeVoronoi::loopCells() { // invoke voro++ and fetch results for owned atoms in group + voronoicell_neighbor c; int i; if (faces_flag) nfaces = 0; From 299eda8ca36eb4f1eda63f31bd021dca479a4ba3 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 17 Aug 2023 16:12:14 -0600 Subject: [PATCH 03/22] have compute_reduce require either peratom or local inputs --- doc/src/compute_reduce.rst | 134 +++++++++++++++++++--------------- doc/src/fix_rigid.rst | 5 +- src/compute_reduce.cpp | 95 +++++++++++++++--------- src/compute_reduce.h | 3 +- src/compute_reduce_region.cpp | 18 +++-- 5 files changed, 149 insertions(+), 106 deletions(-) diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 204f1c090d7..31591d4419e 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -37,13 +37,16 @@ Syntax v_name = per-atom vector calculated by an atom-style variable with name * zero or more keyword/args pairs may be appended -* keyword = *replace* +* keyword = *replace* or *inputs* .. parsed-literal:: *replace* args = vec1 vec2 vec1 = reduced value from this input vector will be replaced vec2 = replace it with vec1[N] where N is index of max/min value from vec2 + *inputs* arg = peratom or local + peratom = all inputs are per-atom quantities (default) + local = all input are local quantities Examples """""""" @@ -61,26 +64,30 @@ Description Define a calculation that "reduces" one or more vector inputs into scalar values, one per listed input. The inputs can be per-atom or -local quantities; they cannot be global quantities. Atom attributes -are per-atom quantities, :doc:`computes ` and :doc:`fixes ` -may generate any of the three kinds of quantities, and :doc:`atom-style variables ` generate per-atom quantities. See the -:doc:`variable ` command and its special functions which can -perform the same operations as the compute reduce command on global -vectors. +local quantities and must all be the same kind (per-atom or local); +see discussion of the optional *inputs* keyword below. + +Atom attributes are per-atom quantities, :doc:`computes ` and +:doc:`fixes ` can generate either per-atom or local quantities, +and :doc:`atom-style variables ` generate per-atom +quantities. See the :doc:`variable ` command and its +special functions which can perform the same reduction operations as +the compute reduce command on global vectors. The reduction operation is specified by the *mode* setting. The *sum* option adds the values in the vector into a global total. The *min* or *max* options find the minimum or maximum value across all vector values. The *minabs* or *maxabs* options find the minimum or maximum value across all absolute vector values. The *ave* setting adds the -vector values into a global total, then divides by the number of values -in the vector. The *sumsq* option sums the square of the values in the -vector into a global total. The *avesq* setting does the same as *sumsq*, -then divides the sum of squares by the number of values. The last two options -can be useful for calculating the variance of some quantity (e.g., variance = -sumsq :math:`-` ave\ :math:`^2`). The *sumabs* option sums the absolute -values in the vector into a global total. The *aveabs* setting does the same -as *sumabs*, then divides the sum of absolute values by the number of +vector values into a global total, then divides by the number of +values in the vector. The *sumsq* option sums the square of the +values in the vector into a global total. The *avesq* setting does +the same as *sumsq*, then divides the sum of squares by the number of +values. The last two options can be useful for calculating the +variance of some quantity (e.g., variance = sumsq :math:`-` ave\ +:math:`^2`). The *sumabs* option sums the absolute values in the +vector into a global total. The *aveabs* setting does the same as +*sumabs*, then divides the sum of absolute values by the number of values. Each listed input is operated on independently. For per-atom inputs, @@ -123,52 +130,54 @@ array with six columns: ---------- -The atom attribute values (*x*, *y*, *z*, *vx*, *vy*, *vz*, *fx*, *fy*, and -*fz*) are self-explanatory. Note that other atom attributes can be used as -inputs to this fix by using the -:doc:`compute property/atom ` command and then specifying -an input value from that compute. +The atom attribute values (*x*, *y*, *z*, *vx*, *vy*, *vz*, *fx*, +*fy*, and *fz*) are self-explanatory. Note that other atom attributes +can be used as inputs to this fix by using the :doc:`compute +property/atom ` command and then specifying an +input value from that compute. If a value begins with "c\_", a compute ID must follow which has been -previously defined in the input script. Computes can generate -per-atom or local quantities. See the individual -:doc:`compute ` page for details. If no bracketed integer -is appended, the vector calculated by the compute is used. If a -bracketed integer is appended, the Ith column of the array calculated -by the compute is used. Users can also write code for their own -compute styles and :doc:`add them to LAMMPS `. See the -discussion above for how :math:`I` can be specified with a wildcard asterisk -to effectively specify multiple values. +previously defined in the input script. Valid computes can generate +per-atom or local quantities. See the individual :doc:`compute +` page for details. If no bracketed integer is appended, the +vector calculated by the compute is used. If a bracketed integer is +appended, the Ith column of the array calculated by the compute is +used. Users can also write code for their own compute styles and +:doc:`add them to LAMMPS `. See the discussion above for how +:math:`I` can be specified with a wildcard asterisk to effectively +specify multiple values. If a value begins with "f\_", a fix ID must follow which has been -previously defined in the input script. Fixes can generate per-atom -or local quantities. See the individual :doc:`fix ` page for -details. Note that some fixes only produce their values on certain -timesteps, which must be compatible with when compute reduce +previously defined in the input script. Valid fixes can generate +per-atom or local quantities. See the individual :doc:`fix ` +page for details. Note that some fixes only produce their values on +certain timesteps, which must be compatible with when compute reduce references the values, else an error results. If no bracketed integer is appended, the vector calculated by the fix is used. If a bracketed integer is appended, the Ith column of the array calculated by the fix is used. Users can also write code for their own fix style and :doc:`add them to LAMMPS `. See the discussion above for how -:math:`I` can be specified with a wildcard asterisk to effectively specify -multiple values. +:math:`I` can be specified with a wildcard asterisk to effectively +specify multiple values. If a value begins with "v\_", a variable name must follow which has been previously defined in the input script. It must be an :doc:`atom-style variable `. Atom-style variables can reference thermodynamic keywords and various per-atom attributes, or invoke other computes, fixes, or variables when they are evaluated, so -this is a very general means of generating per-atom quantities to reduce. +this is a very general means of generating per-atom quantities to +reduce. ---------- If the *replace* keyword is used, two indices *vec1* and *vec2* are -specified, where each index ranges from 1 to the number of input values. -The replace keyword can only be used if the *mode* is *min* or *max*\ . -It works as follows. A min/max is computed as usual on the *vec2* -input vector. The index :math:`N` of that value within *vec2* is also stored. -Then, instead of performing a min/max on the *vec1* input vector, the -stored index is used to select the :math:`N`\ th element of the *vec1* vector. +specified, where each index ranges from 1 to the number of input +values. The replace keyword can only be used if the *mode* is *min* +or *max*\ . It works as follows. A min/max is computed as usual on +the *vec2* input vector. The index :math:`N` of that value within +*vec2* is also stored. Then, instead of performing a min/max on the +*vec1* input vector, the stored index is used to select the :math:`N`\ +th element of the *vec1* vector. Thus, for example, if you wish to use this compute to find the bond with maximum stretch, you can do it as follows: @@ -190,6 +199,14 @@ information in this context, the *replace* keywords will extract the atom IDs for the two atoms in the bond of maximum stretch. These atom IDs and the bond stretch will be printed with thermodynamic output. +The *inputs* keyword allows selection of whether all the inputs are +per-atom or local quantities. As noted above, all the inputs must be +the same kind (per-atom or local). Per-atom is the default setting. +If a compute or fix is specified as an input, it must produce per-atom +or local data to match this setting. If it produces both, e.g. for +the :doc:`compute voronoi/atom ` command, then +this keyword selects between them. + ---------- If a single input is specified this compute produces a global scalar @@ -197,34 +214,35 @@ value. If multiple inputs are specified, this compute produces a global vector of values, the length of which is equal to the number of inputs specified. -As discussed below, for the *sum*, *sumabs*, and *sumsq* modes, the value(s) -produced by this compute are all "extensive", meaning their value -scales linearly with the number of atoms involved. If normalized -values are desired, this compute can be accessed by the +As discussed below, for the *sum*, *sumabs*, and *sumsq* modes, the +value(s) produced by this compute are all "extensive", meaning their +value scales linearly with the number of atoms involved. If +normalized values are desired, this compute can be accessed by the :doc:`thermo_style custom ` command with -:doc:`thermo_modify norm yes ` set as an option. -Or it can be accessed by a -:doc:`variable ` that divides by the appropriate atom count. +:doc:`thermo_modify norm yes ` set as an option. Or it +can be accessed by a :doc:`variable ` that divides by the +appropriate atom count. ---------- Output info """"""""""" -This compute calculates a global scalar if a single input value is specified -or a global vector of length :math:`N`, where :math:`N` is the number of -inputs, and which can be accessed by indices 1 to :math:`N`. These values can -be used by any command that uses global scalar or vector values from a -compute as input. See the :doc:`Howto output ` doc page -for an overview of LAMMPS output options. +This compute calculates a global scalar if a single input value is +specified or a global vector of length :math:`N`, where :math:`N` is +the number of inputs, and which can be accessed by indices 1 to +:math:`N`. These values can be used by any command that uses global +scalar or vector values from a compute as input. See the :doc:`Howto +output ` doc page for an overview of LAMMPS output +options. All the scalar or vector values calculated by this compute are "intensive", except when the *sum*, *sumabs*, or *sumsq* modes are used on per-atom or local vectors, in which case the calculated values are "extensive". -The scalar or vector values will be in whatever :doc:`units ` the -quantities being reduced are in. +The scalar or vector values will be in whatever :doc:`units ` +the quantities being reduced are in. Restrictions """""""""""" @@ -238,4 +256,4 @@ Related commands Default """"""" -none +The default value for the *inputs* keyword is peratom. diff --git a/doc/src/fix_rigid.rst b/doc/src/fix_rigid.rst index 89759da817f..a50e2156812 100644 --- a/doc/src/fix_rigid.rst +++ b/doc/src/fix_rigid.rst @@ -843,7 +843,7 @@ stress/atom ` commands. The former can be accessed by :doc:`thermodynamic output `. The default setting for this fix is :doc:`fix_modify virial yes `. -All of the *rigid* styles (not the *rigid/small* styles) compute a +All of the *rigid* styles (but not the *rigid/small* styles) compute a global array of values which can be accessed by various :doc:`output commands `. Similar information about the bodies defined by the *rigid/small* styles can be accessed via the @@ -887,7 +887,8 @@ Restrictions """""""""""" These fixes are all part of the RIGID package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. Assigning a temperature via the :doc:`velocity create ` command to a system with :doc:`rigid bodies ` may not have diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 6b27498eb70..8565ddb1c9b 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -31,12 +31,16 @@ using namespace LAMMPS_NS; +enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduceRegion + #define BIG 1.0e20 //---------------------------------------------------------------- + void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) { // r is the already reduced value, n is the new value + double n = std::fabs(*(double *) in), r = *(double *) inout; double m; @@ -47,9 +51,11 @@ void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) } *(double *) inout = m; } + void abs_min(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) { // r is the already reduced value, n is the new value + double n = std::fabs(*(double *) in), r = *(double *) inout; double m; @@ -68,6 +74,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : owner(nullptr), idregion(nullptr), region(nullptr), varatom(nullptr) { int iarg = 0; + if (strcmp(style, "reduce") == 0) { if (narg < 5) utils::missing_cmd_args(FLERR, "compute reduce", error); iarg = 3; @@ -128,42 +135,52 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : // parse values + input_mode = UNDECIDED; + values.clear(); nvalues = 0; for (int iarg = 0; iarg < nargnew; ++iarg) { value_t val; val.id = ""; - val.flavor = 0; val.val.c = nullptr; if (strcmp(arg[iarg], "x") == 0) { + input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 0; } else if (strcmp(arg[iarg], "y") == 0) { + input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 1; } else if (strcmp(arg[iarg], "z") == 0) { + input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 2; } else if (strcmp(arg[iarg], "vx") == 0) { + input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 0; } else if (strcmp(arg[iarg], "vy") == 0) { + input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 1; } else if (strcmp(arg[iarg], "vz") == 0) { + input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 2; } else if (strcmp(arg[iarg], "fx") == 0) { + input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 0; } else if (strcmp(arg[iarg], "fy") == 0) { + input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 1; } else if (strcmp(arg[iarg], "fz") == 0) { + input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 2; @@ -207,6 +224,14 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Compute {} replace column already used for another replacement"); replace[col1] = col2; iarg += 2; + } else if (strcmp(arg[iarg], "inputs") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, mycmd + " inputs", error); + if (strcmp(arg[iarg+1], "peratom") == 0) input_mode = PERATOM; + else if (strcmp(arg[iarg+1], "local") == 0) { + if (input_mode == PERATOM) + error->all(FLERR,"Compute {} inputs must be all peratom or all local"); + input_mode = LOCAL; + } } else error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]); } @@ -231,66 +256,64 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : // setup and error check for (auto &val : values) { - if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) - val.flavor = PERATOM; - - else if (val.which == ArgInfo::COMPUTE) { + if (val.which == ArgInfo::COMPUTE) { val.val.c = modify->get_compute_by_id(val.id); if (!val.val.c) error->all(FLERR, "Compute ID {} for compute {} does not exist", val.id, style); - if (val.val.c->peratom_flag) { - val.flavor = PERATOM; + + if (input_mode == PERATOM) { + if (!val.val.c->peratom_flag) + error->all(FLERR, "Compute {} compute {} does not calculate per-atom values", style, val.id); if (val.argindex == 0 && val.val.c->size_peratom_cols != 0) - error->all(FLERR, "Compute {} compute {} does not calculate a per-atom vector", style, - val.id); + error->all(FLERR, "Compute {} compute {} does not calculate a per-atom vector", style, val.id); if (val.argindex && val.val.c->size_peratom_cols == 0) - error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style, - val.id); + error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style, val.id); if (val.argindex && val.argindex > val.val.c->size_peratom_cols) error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id); - } else if (val.val.c->local_flag) { - val.flavor = LOCAL; + + } else if (input_mode == LOCAL) { + if (!val.val.c->peratom_flag) + error->all(FLERR, "Compute {} compute {} does not calculate local values", style, val.id); if (val.argindex == 0 && val.val.c->size_local_cols != 0) - error->all(FLERR, "Compute {} compute {} does not calculate a local vector", style, - val.id); + error->all(FLERR, "Compute {} compute {} does not calculate a local vector", style, val.id); if (val.argindex && val.val.c->size_local_cols == 0) - error->all(FLERR, "Compute {} compute {} does not calculate a local array", style, - val.id); + error->all(FLERR, "Compute {} compute {} does not calculate a local array", style, val.id); if (val.argindex && val.argindex > val.val.c->size_local_cols) error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id); - } else - error->all(FLERR, "Compute {} compute {} calculates global values", style, val.id); + } } else if (val.which == ArgInfo::FIX) { val.val.f = modify->get_fix_by_id(val.id); if (!val.val.f) error->all(FLERR, "Fix ID {} for compute {} does not exist", val.id, style); - if (val.val.f->peratom_flag) { - val.flavor = PERATOM; + + if (input_mode == PERATOM) { + if (!val.val.f->peratom_flag) + error->all(FLERR, "Compute {} fix {} does not calculate per-atom values", style, val.id); if (val.argindex == 0 && (val.val.f->size_peratom_cols != 0)) - error->all(FLERR, "Compute {} fix {} does not calculate a per-atom vector", style, - val.id); + error->all(FLERR, "Compute {} fix {} does not calculate a per-atom vector", style, val.id); if (val.argindex && (val.val.f->size_peratom_cols == 0)) error->all(FLERR, "Compute {} fix {} does not calculate a per-atom array", style, val.id); if (val.argindex && (val.argindex > val.val.f->size_peratom_cols)) error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id); - } else if (val.val.f->local_flag) { - val.flavor = LOCAL; + + } else if (input_mode == LOCAL) { + if (!val.val.f->local_flag) + error->all(FLERR, "Compute {} fix {} does not calculate local values", style, val.id); if (val.argindex == 0 && (val.val.f->size_local_cols != 0)) error->all(FLERR, "Compute {} fix {} does not calculate a local vector", style, val.id); if (val.argindex && (val.val.f->size_local_cols == 0)) error->all(FLERR, "Compute {} fix {} does not calculate a local array", style, val.id); if (val.argindex && (val.argindex > val.val.f->size_local_cols)) error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id); - } else - error->all(FLERR, "Compute {} fix {} calculates global values", style, val.id); + } } else if (val.which == ArgInfo::VARIABLE) { + if (input_mode == LOCAL) error->all(FLERR,"Compute {} inputs must be all local"); val.val.v = input->variable->find(val.id.c_str()); if (val.val.v < 0) error->all(FLERR, "Variable name {} for compute {} does not exist", val.id, style); if (input->variable->atomstyle(val.val.v) == 0) error->all(FLERR, "Compute {} variable {} is not atom-style variable", style, val.id); - val.flavor = PERATOM; } } @@ -512,7 +535,7 @@ double ComputeReduce::compute_one(int m, int flag) } else if (val.which == ArgInfo::COMPUTE) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { val.val.c->compute_peratom(); val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; @@ -537,7 +560,7 @@ double ComputeReduce::compute_one(int m, int flag) one = carray_atom[flag][aidxm1]; } - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { val.val.c->compute_local(); val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; @@ -567,7 +590,7 @@ double ComputeReduce::compute_one(int m, int flag) if (update->ntimestep % val.val.f->peratom_freq) error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style); - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { if (aidx == 0) { double *fix_vector = val.val.f->vector_atom; if (flag < 0) { @@ -585,7 +608,7 @@ double ComputeReduce::compute_one(int m, int flag) one = fix_array[flag][aidxm1]; } - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { if (aidx == 0) { double *fix_vector = val.val.f->vector_local; int n = val.val.f->size_local_rows; @@ -632,18 +655,18 @@ bigint ComputeReduce::count(int m) if ((val.which == ArgInfo::X) || (val.which == ArgInfo::V) || (val.which == ArgInfo::F)) return group->count(igroup); else if (val.which == ArgInfo::COMPUTE) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { return group->count(igroup); - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { bigint ncount = val.val.c->size_local_rows; bigint ncountall; MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world); return ncountall; } } else if (val.which == ArgInfo::FIX) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { return group->count(igroup); - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { bigint ncount = val.val.f->size_local_rows; bigint ncountall; MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world); diff --git a/src/compute_reduce.h b/src/compute_reduce.h index f8f73cb17a6..f8b652e00c6 100644 --- a/src/compute_reduce.h +++ b/src/compute_reduce.h @@ -37,12 +37,11 @@ class ComputeReduce : public Compute { double memory_usage() override; protected: - int mode, nvalues; + int mode, nvalues, input_mode; struct value_t { int which; int argindex; std::string id; - int flavor; union { class Compute *c; class Fix *f; diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp index efce00ff663..2f5a3de6751 100644 --- a/src/compute_reduce_region.cpp +++ b/src/compute_reduce_region.cpp @@ -26,6 +26,8 @@ using namespace LAMMPS_NS; +enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduce + static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ @@ -97,7 +99,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) // invoke compute if not previously invoked } else if (val.which == ArgInfo::COMPUTE) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { val.val.c->compute_peratom(); val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; @@ -122,7 +124,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) one = compute_array[flag][aidxm1]; } - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { val.val.c->compute_local(); val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; @@ -151,7 +153,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) if (update->ntimestep % val.val.f->peratom_freq) error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style); - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { if (aidx == 0) { double *fix_vector = val.val.f->vector_atom; if (flag < 0) { @@ -171,7 +173,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) one = fix_array[flag][aidxm1]; } - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { if (aidx == 0) { double *fix_vector = val.val.f->vector_local; if (flag < 0) @@ -219,18 +221,18 @@ bigint ComputeReduceRegion::count(int m) if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) return group->count(igroup, region); else if (val.which == ArgInfo::COMPUTE) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { return group->count(igroup, region); - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { bigint ncount = val.val.c->size_local_rows; bigint ncountall; MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world); return ncountall; } } else if (val.which == ArgInfo::FIX) { - if (val.flavor == PERATOM) { + if (input_mode == PERATOM) { return group->count(igroup, region); - } else if (val.flavor == LOCAL) { + } else if (input_mode == LOCAL) { bigint ncount = val.val.f->size_local_rows; bigint ncountall; MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world); From f2901827e6cc74bdacfe8cfb06d1ba62322ed007 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 17 Aug 2023 17:25:27 -0600 Subject: [PATCH 04/22] updates to variable doc page to clarify compute/fix options --- doc/src/compute.rst | 4 +- doc/src/thermo_style.rst | 53 +++++---- doc/src/variable.rst | 246 +++++++++++++++++++++------------------ 3 files changed, 162 insertions(+), 141 deletions(-) diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 317efa4f413..226dc6373bc 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -43,8 +43,8 @@ underscores. ---------- -Computes calculate one or more of four styles of quantities: global, -per-atom, local, or per-atom. A global quantity is one or more +Computes calculate and store any of four styles of quantities: global, +per-atom, local, or per-grid. A global quantity is one or more system-wide values, e.g. the temperature of the system. A per-atom quantity is one or more values per atom, e.g. the kinetic energy of each atom. Per-atom values are set to 0.0 for atoms not in the diff --git a/doc/src/thermo_style.rst b/doc/src/thermo_style.rst index 63ad59e553f..c3c607a4796 100644 --- a/doc/src/thermo_style.rst +++ b/doc/src/thermo_style.rst @@ -385,19 +385,20 @@ creates a global vector with 6 values. The *c_ID* and *c_ID[I]* and *c_ID[I][J]* keywords allow global values calculated by a compute to be output. As discussed on the :doc:`compute ` doc page, computes can calculate global, -per-atom, or local values. Only global values can be referenced by -this command. However, per-atom compute values for an individual atom -can be referenced in a :doc:`variable ` and the variable -referenced by thermo_style custom, as discussed below. See the -discussion above for how the I in *c_ID[I]* can be specified with a -wildcard asterisk to effectively specify multiple values from a global -compute vector. +per-atom, local, and per-grid values. Only global values can be +referenced by this command. However, per-atom compute values for an +individual atom can be referenced in a :doc:`equal-style variable +` and the variable referenced by thermo_style custom, as +discussed below. See the discussion above for how the I in *c_ID[I]* +can be specified with a wildcard asterisk to effectively specify +multiple values from a global compute vector. The ID in the keyword should be replaced by the actual ID of a compute that has been defined elsewhere in the input script. See the -:doc:`compute ` command for details. If the compute calculates -a global scalar, vector, or array, then the keyword formats with 0, 1, -or 2 brackets will reference a scalar value from the compute. +:doc:`compute ` command for details. If the compute +calculates a global scalar, vector, or array, then the keyword formats +with 0, 1, or 2 brackets will reference a scalar value from the +compute. Note that some computes calculate "intensive" global quantities like temperature; others calculate "extensive" global quantities like @@ -410,13 +411,14 @@ norm ` option being used. The *f_ID* and *f_ID[I]* and *f_ID[I][J]* keywords allow global values calculated by a fix to be output. As discussed on the :doc:`fix -` doc page, fixes can calculate global, per-atom, or local -values. Only global values can be referenced by this command. -However, per-atom fix values can be referenced for an individual atom -in a :doc:`variable ` and the variable referenced by -thermo_style custom, as discussed below. See the discussion above for -how the I in *f_ID[I]* can be specified with a wildcard asterisk to -effectively specify multiple values from a global fix vector. +` doc page, fixes can calculate global, per-atom, local, and +per-grid values. Only global values can be referenced by this +command. However, per-atom fix values can be referenced for an +individual atom in a :doc:`equal-style variable ` and the +variable referenced by thermo_style custom, as discussed below. See +the discussion above for how the I in *f_ID[I]* can be specified with +a wildcard asterisk to effectively specify multiple values from a +global fix vector. The ID in the keyword should be replaced by the actual ID of a fix that has been defined elsewhere in the input script. See the @@ -438,14 +440,15 @@ output. The name in the keyword should be replaced by the variable name that has been defined elsewhere in the input script. Only equal-style and vector-style variables can be referenced; the latter requires a bracketed term to specify the Ith element of the vector -calculated by the variable. However, an atom-style variable can be -referenced for an individual atom by an equal-style variable and that -variable referenced. See the :doc:`variable ` command for -details. Variables of style *equal* and *vector* and *atom* define a -formula which can reference per-atom properties or thermodynamic -keywords, or they can invoke other computes, fixes, or variables when -evaluated, so this is a very general means of creating thermodynamic -output. +calculated by the variable. However, an equal-style variable can use +an atom-style variable in its formula indexed by the ID of an +individual atom. This is a way to output a speciic atom's per-atom +coordinates or other per-atom properties in thermo output. See the +:doc:`variable ` command for details. Note that variables +of style *equal* and *vector* and *atom* define a formula which can +reference per-atom properties or thermodynamic keywords, or they can +invoke other computes, fixes, or variables when evaluated, so this is +a very general means of creating thermodynamic output. Note that equal-style and vector-style variables are assumed to produce "intensive" global quantities, which are thus printed as-is, diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 28c0d297995..38e423b6328 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -550,12 +550,11 @@ variables. Most of the formula elements produce a scalar value. Some produce a global or per-atom vector of values. Global vectors can be produced by computes or fixes or by other vector-style variables. Per-atom -vectors are produced by atom vectors, compute references that -represent a per-atom vector, fix references that represent a per-atom -vector, and variables that are atom-style variables. Math functions -that operate on scalar values produce a scalar value; math function -that operate on global or per-atom vectors do so element-by-element -and produce a global or per-atom vector. +vectors are produced by atom vectors, computes or fixes which output a +per-atom vector or array, and variables that are atom-style variables. +Math functions that operate on scalar values produce a scalar value; +math function that operate on global or per-atom vectors do so +element-by-element and produce a global or per-atom vector. A formula for equal-style variables cannot use any formula element that produces a global or per-atom vector. A formula for a @@ -564,12 +563,13 @@ scalar value or a global vector value, but cannot use a formula element that produces a per-atom vector. A formula for an atom-style variable can use formula elements that produce either a scalar value or a per-atom vector, but not one that produces a global vector. + Atom-style variables are evaluated by other commands that define a -:doc:`group ` on which they operate, e.g. a :doc:`dump ` or -:doc:`compute ` or :doc:`fix ` command. When they invoke -the atom-style variable, only atoms in the group are included in the -formula evaluation. The variable evaluates to 0.0 for atoms not in -the group. +:doc:`group ` on which they operate, e.g. a :doc:`dump ` +or :doc:`compute ` or :doc:`fix ` command. When they +invoke the atom-style variable, only atoms in the group are included +in the formula evaluation. The variable evaluates to 0.0 for atoms +not in the group. ---------- @@ -1138,69 +1138,74 @@ only defined if an :doc:`atom_style ` is being used that defines molecule IDs. Note that many other atom attributes can be used as inputs to a -variable by using the :doc:`compute property/atom ` command and then specifying -a quantity from that compute. +variable by using the :doc:`compute property/atom +` command and then specifying a quantity from +that compute. ---------- Compute References ------------------ -Compute references access quantities calculated by a -:doc:`compute `. The ID in the reference should be replaced by -the ID of a compute defined elsewhere in the input script. As -discussed in the page for the :doc:`compute ` command, -computes can produce global, per-atom, or local values. Only global -and per-atom values can be used in a variable. Computes can also -produce a scalar, vector, or array. - -An equal-style variable can only use scalar values, which means a -global scalar, or an element of a global or per-atom vector or array. -A vector-style variable can use scalar values or a global vector of -values, or a column of a global array of values. Atom-style variables -can use global scalar values. They can also use per-atom vector -values, or a column of a per-atom array. See the doc pages for -individual computes to see what kind of values they produce. - -Examples of different kinds of compute references are as follows. -There is typically no ambiguity (see exception below) as to what a -reference means, since computes only produce either global or per-atom -quantities, never both. - -+-------------+-------------------------------------------------------------------------------------------------------+ -| c_ID | global scalar, or per-atom vector | -+-------------+-------------------------------------------------------------------------------------------------------+ -| c_ID[I] | Ith element of global vector, or atom I's value in per-atom vector, or Ith column from per-atom array | -+-------------+-------------------------------------------------------------------------------------------------------+ -| c_ID[I][J] | I,J element of global array, or atom I's Jth value in per-atom array | -+-------------+-------------------------------------------------------------------------------------------------------+ - -For I and J indices, integers can be specified or a variable name, -specified as v_name, where name is the name of the variable. The -rules for this syntax are the same as for the "Atom Values and -Vectors" discussion above. - -One source of ambiguity for compute references is when a vector-style -variable refers to a compute that produces both a global scalar and a -global vector. Consider a compute with ID "foo" that does this, -referenced as follows by variable "a", where "myVec" is another -vector-style variable: - -.. code-block:: LAMMPS - - variable a vector c_foo*v_myVec - -The reference "c_foo" could refer to either the global scalar or -global vector produced by compute "foo". In this case, "c_foo" will -always refer to the global scalar, and "C_foo" can be used to -reference the global vector. Similarly if the compute produces both a -global vector and global array, then "c_foo[I]" will always refer to -an element of the global vector, and "C_foo[I]" can be used to -reference the Ith column of the global array. - -Note that if a variable containing a compute is evaluated directly in -an input script (not during a run), then the values accessed by the -compute must be current. See the discussion below about "Variable +Compute references access quantities calculated by a :doc:`compute +`. The ID in the reference should be replaced by the ID of a +compute defined elsewhere in the input script. + +As discussed on the page for the :doc:`compute ` command, +computes can produce global, per-atom, local, and per-grid values. +Only global and per-atom values can be used in a variable. Computes +can also produce scalars (global only), vectors, and arrays. See the +doc pages for individual computes to see what different kinds of data +they produce. + +An equal-style variable can only use scalar values, either from global +or per-atom data. In the case of per-atom data, this would be a value +for a specific atom. + +A vector-style variable can use scalar values (same as for equal-style +variables), or global vectors of values. The latter can also be a +column of a global array. + +Atom-style variables can use scalar values (same as for equal-style +varaibles), or per-atom vectors of values. The latter can also be a +column of a per-atom array. + +The various allowed compute references in the variable formulas for +equal-, vector-, and atom-style variables are listed in the following +table: + ++--------+------------+--------------------------------------------+ +| equal | c_ID | global scalar | +| equal | c_ID[I] | element of global vector | +| equal | c_ID[I][J] | element of global array | +| equal | C_ID[I] | element of per-atom vector, I = ID of atom | +| equal | C_ID{i}[J] | element of per-atom array, I = ID of atom | ++--------+------------+--------------------------------------------| +| vector | c_ID | global vector | +| vector | c_ID[I] | column of global array | +---------+------------+--------------------------------------------+ +| atom | c_ID | per-atom vector | +| atom | c_ID[I] | column of per-atom array | ++--------+------------+--------------------------------------------+ + +Note that if an equal-style variable formula wishes to access per-atom +data from a compute, it must use capital "C" as the ID prefix and not +lower-case "c". + +Also note that if a vector- or atom-style variable formula needs to +access a scalar value from a compute (i.e. the 5 kinds of values in +the first 5 lines of the table), it can not do so directly. Instead, +it can use a reference to an equal-style variable which stores the +scalar value from the compute. + +The I and J indices in these compute references can be integers or can +be a variable name, specified as v_name, where name is the name of the +variable. The rules for this syntax are the same as for indices in +the "Atom Values and Vectors" discussion above. + +If a variable containing a compute is evaluated directly in an input +script (not during a run), then the values accessed by the compute +should be current. See the discussion below about "Variable Accuracy". ---------- @@ -1208,51 +1213,60 @@ Accuracy". Fix References -------------- -Fix references access quantities calculated by a :doc:`fix `. +Fix references access quantities calculated by a :doc:`fix `. The ID in the reference should be replaced by the ID of a fix defined -elsewhere in the input script. As discussed in the page for the -:doc:`fix ` command, fixes can produce global, per-atom, or local -values. Only global and per-atom values can be used in a variable. -Fixes can also produce a scalar, vector, or array. An equal-style -variable can only use scalar values, which means a global scalar, or -an element of a global or per-atom vector or array. Atom-style -variables can use the same scalar values. They can also use per-atom -vector values. A vector value can be a per-atom vector itself, or a -column of an per-atom array. See the doc pages for individual fixes -to see what kind of values they produce. - -The different kinds of fix references are exactly the same as the -compute references listed in the above table, where "c\_" is replaced -by "f\_". Again, there is typically no ambiguity (see exception below) -as to what a reference means, since fixes only produce either global -or per-atom quantities, never both. - -+-------------+-------------------------------------------------------------------------------------------------------+ -| f_ID | global scalar, or per-atom vector | -+-------------+-------------------------------------------------------------------------------------------------------+ -| f_ID[I] | Ith element of global vector, or atom I's value in per-atom vector, or Ith column from per-atom array | -+-------------+-------------------------------------------------------------------------------------------------------+ -| f_ID[I][J] | I,J element of global array, or atom I's Jth value in per-atom array | -+-------------+-------------------------------------------------------------------------------------------------------+ - -For I and J indices, integers can be specified or a variable name, -specified as v_name, where name is the name of the variable. The -rules for this syntax are the same as for the "Atom Values and -Vectors" discussion above. - -One source of ambiguity for fix references is the same ambiguity -discussed for compute references above. Namely when a vector-style -variable refers to a fix that produces both a global scalar and a -global vector. The solution is the same as for compute references. -For a fix with ID "foo", "f_foo" will always refer to the global -scalar, and "F_foo" can be used to reference the global vector. And -similarly for distinguishing between a fix's global vector versus -global array with "f_foo[I]" versus "F_foo[I]". - -Note that if a variable containing a fix is evaluated directly in an -input script (not during a run), then the values accessed by the fix -should be current. See the discussion below about "Variable -Accuracy". +elsewhere in the input script. + +As discussed on the page for the :doc:`fix ` command, fixes can +produce global, per-atom, local, and per-grid values. Only global and +per-atom values can be used in a variable. Fixes can also produce +scalars (global only), vectors, and arrays. See the doc pages for +individual fixes to see what different kinds of data they produce. + +An equal-style variable can only use scalar values, either from global +or per-atom data. In the case of per-atom data, this would be a value +for a specific atom. + +A vector-style variable can use scalar values (same as for equal-style +variables), or global vectors of values. The latter can also be a +column of a global array. + +Atom-style variables can use scalar values (same as for equal-style +varaibles), or per-atom vectors of values. The latter can also be a +column of a per-atom array. + +The various allowed fix references in the variable formulas for +equal-, vector-, and atom-style variables are listed in the following +table: + ++--------+------------+--------------------------------------------+ +| equal | f_ID | global scalar | +| equal | f_ID[I] | element of global vector | +| equal | f_ID[I][J] | element of global array | +| equal | F_ID[I] | element of per-atom vector, I = ID of atom | +| equal | F_ID{i}[J] | element of per-atom array, I = ID of atom | ++--------+------------+--------------------------------------------| +| vector | f_ID | global vector | +| vector | f_ID[I] | column of global array | +---------+------------+--------------------------------------------+ +| atom | f_ID | per-atom vector | +| atom | f_ID[I] | column of per-atom array | ++--------+------------+--------------------------------------------+ + +Note that if an equal-style variable formula wishes to access per-atom +data from a fix, it must use capital "F" as the ID prefix and not +lower-case "f". + +Also note that if a vector- or atom-style variable formula needs to +access a scalar value from a fix (i.e. the 5 kinds of values in the +first 5 lines of the table), it can not do so directly. Instead, it +can use a reference to an equal-style variable which stores the scalar +value from the fix. + +The I and J indices in these fix references can be integers or can be +a variable name, specified as v_name, where name is the name of the +variable. The rules for this syntax are the same as for indices in +the "Atom Values and Vectors" discussion above. Note that some fixes only generate quantities on certain timesteps. If a variable attempts to access the fix on non-allowed timesteps, an @@ -1260,6 +1274,10 @@ error is generated. For example, the :doc:`fix ave/time ` command may only generate averaged quantities every 100 steps. See the doc pages for individual fix commands for details. +If a variable containing a fix is evaluated directly in an input +script (not during a run), then the values accessed by the fix should +be current. See the discussion below about "Variable Accuracy". + ---------- Variable References @@ -1312,8 +1330,8 @@ produce only a global scalar or global vector or per-atom vector. For the I index, an integer can be specified or a variable name, specified as v_name, where name is the name of the variable. The -rules for this syntax are the same as for the "Atom Values and -Vectors" discussion above. +rules for this syntax are the same as for indices in the "Atom Values +and Vectors" discussion above. ---------- From 95e9e6549f6a2658e7f361c4e40616478f274856 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 18 Aug 2023 09:28:58 -0600 Subject: [PATCH 05/22] simply variable.cpp --- src/variable.cpp | 351 ++++++++++++++++++++++++----------------------- 1 file changed, 182 insertions(+), 169 deletions(-) diff --git a/src/variable.cpp b/src/variable.cpp index cf2e5c3b6fe..5013f3ce550 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1469,8 +1469,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (domain->box_exist == 0) print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar); - // uppercase used to force access of - // global vector vs global scalar, and global array vs global vector + // uppercase used to access of peratom data by equal-style var int lowercase = 1; if (word[0] == 'C') lowercase = 0; @@ -1479,7 +1478,6 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!compute) print_var_error(FLERR,fmt::format("Invalid compute ID '{}' in variable formula", word+2),ivar); - // parse zero or one or two trailing brackets // point i beyond last bracket // nbracket = # of bracket pairs @@ -1501,107 +1499,203 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) } } - // c_ID = scalar from global scalar, must be lowercase + // equal-style variable is being evaluated - if (nbracket == 0 && compute->scalar_flag && lowercase) { + if (style[ivar] == EQUAL) { + + // c_ID = scalar from global scalar - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { - compute->compute_scalar(); - compute->invoked_flag |= Compute::INVOKED_SCALAR; - } + if (lowercase && nbracket == 0) { - value1 = compute->scalar; - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; + if (!compute->scalar_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); - // c_ID[i] = scalar from global vector, must be lowercase + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { + compute->compute_scalar(); + compute->invoked_flag |= Compute::INVOKED_SCALAR; + } - } else if (nbracket == 1 && compute->vector_flag && lowercase) { + value1 = compute->scalar; + argstack[nargstack++] = value1; - if (index1 > compute->size_vector && - compute->size_vector_variable == 0) - print_var_error(FLERR,"Variable formula compute vector is accessed out-of-range",ivar,0); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; - } + // c_ID[i] = scalar from global vector + + } else if (lowercase && nbracket == 1) { + + if (!compute->vector_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (index1 > compute->size_vector && + compute->size_vector_variable == 0) + print_var_error(FLERR,"Variable formula compute vector is accessed out-of-range",ivar,0); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= Compute::INVOKED_VECTOR; + } if (compute->size_vector_variable && index1 > compute->size_vector) value1 = 0.0; else value1 = compute->vector[index1-1]; - if (tree) { + argstack[nargstack++] = value1; + + // c_ID[i][j] = scalar from global array + + } else if (lowercase && nbracket == 2) { + + if (!compute->array_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (index1 > compute->size_array_rows && + compute->size_array_rows_variable == 0) + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); + if (index2 > compute->size_array_cols) + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { + compute->compute_array(); + compute->invoked_flag |= Compute::INVOKED_ARRAY; + } + + if (compute->size_array_rows_variable && + index1 > compute->size_array_rows) value1 = 0.0; + else value1 = compute->array[index1-1][index2-1]; + argstack[nargstack++] = value1; + + // C_ID[i] = scalar element of per-atom vector + + } else if (!lowercase && nbracket == 1) { + + if (!compute->peratom_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (compute->size_peratom_cols) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= Compute::INVOKED_PERATOM; + } + + peratom2global(1,nullptr,compute->vector_atom,1,index1,tree, + treestack,ntreestack,argstack,nargstack); + + // C_ID[i][j] = scalar element of per-atom array + + } else if (!lowercase && nbracket == 2) { + + if (!compute->peratom_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (!compute->size_peratom_cols) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (index2 > compute->size_peratom_cols) + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= Compute::INVOKED_PERATOM; + } + + if (compute->array_atom) + peratom2global(1,nullptr,&compute->array_atom[0][index2-1], + compute->size_peratom_cols,index1, + tree,treestack,ntreestack,argstack,nargstack); + else + peratom2global(1,nullptr,nullptr,compute->size_peratom_cols,index1, + tree,treestack,ntreestack,argstack,nargstack); + + // no other possibilities for equal-style variable, so error + + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + + // vector-style variable is being evaluated + + } else if (style[ivar] == VECTOR) { + + // c_ID = vector from global vector + + if (lowercase && nbracket == 0) { + + if (!compute->vector_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + + // c_ID[i] = vector from global array + + } else if (lowercase && nbracket == 1) { + + if (!compute->array_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + + // no other possibilities for vector-style variable, so error + + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + + // atom-style variable is being evaluated + + } else if (style[ivar] == ATOM) { + + // c_ID = vector from per-atom vector + + if (lowercase && nbracket == 0) { + + if (!compute->peratom_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (compute->size_peratom_cols) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= Compute::INVOKED_PERATOM; + } + auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; + newtree->type = ATOMARRAY; + newtree->array = compute->vector_atom; + newtree->nstride = 1; treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; + + // c_ID[i] = vector from per-atom array + + } else if (lowercase && nbracket == 1) { + + if (!compute->peratom_flag) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (!compute->size_peratom_cols) + print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (index1 > compute->size_peratom_cols) + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= Compute::INVOKED_PERATOM; + } - // c_ID[i][j] = scalar from global array, must be lowercase - - } else if (nbracket == 2 && compute->array_flag && lowercase) { - - if (index1 > compute->size_array_rows && - compute->size_array_rows_variable == 0) - print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); - if (index2 > compute->size_array_cols) - print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { - compute->compute_array(); - compute->invoked_flag |= Compute::INVOKED_ARRAY; - } - - if (compute->size_array_rows_variable && - index1 > compute->size_array_rows) value1 = 0.0; - else value1 = compute->array[index1-1][index2-1]; - if (tree) { auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; + newtree->type = ATOMARRAY; + newtree->array = nullptr; + if (compute->array_atom) + newtree->array = &compute->array_atom[0][index1-1]; + newtree->nstride = compute->size_peratom_cols; treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; - - // c_ID = vector from global vector, lowercase or uppercase - - } else if (nbracket == 0 && compute->vector_flag) { - - if (tree == nullptr) - print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar); - if (treetype == ATOM) - print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar); - if (compute->size_vector == 0) - print_var_error(FLERR,"Variable formula compute vector is zero length",ivar); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; - } - - auto newtree = new Tree(); - newtree->type = VECTORARRAY; - newtree->array = compute->vector; - newtree->nvector = compute->size_vector; - newtree->nstride = 1; - treestack[ntreestack++] = newtree; - - // c_ID[i] = vector from global array, lowercase or uppercase - } else if (nbracket == 1 && compute->array_flag) { + // no other possibilities for atom-style variable, so error + + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + } + if (tree == nullptr) print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar); if (treetype == ATOM) @@ -1623,97 +1717,16 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->nstride = compute->size_array_cols; treestack[ntreestack++] = newtree; - // c_ID[i] = scalar from per-atom vector - - } else if (nbracket == 1 && compute->peratom_flag && - compute->size_peratom_cols == 0) { - - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; - } - - peratom2global(1,nullptr,compute->vector_atom,1,index1,tree, - treestack,ntreestack,argstack,nargstack); - // c_ID[i][j] = scalar from per-atom array - } else if (nbracket == 2 && compute->peratom_flag && - compute->size_peratom_cols > 0) { - if (index2 > compute->size_peratom_cols) - print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; - } - if (compute->array_atom) - peratom2global(1,nullptr,&compute->array_atom[0][index2-1],compute->size_peratom_cols,index1, - tree,treestack,ntreestack,argstack,nargstack); - else - peratom2global(1,nullptr,nullptr,compute->size_peratom_cols,index1, - tree,treestack,ntreestack,argstack,nargstack); - // c_ID = vector from per-atom vector - } else if (nbracket == 0 && compute->peratom_flag && - compute->size_peratom_cols == 0) { - if (tree == nullptr) - print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar); - if (treetype == VECTOR) - print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; - } - - auto newtree = new Tree(); - newtree->type = ATOMARRAY; - newtree->array = compute->vector_atom; - newtree->nstride = 1; - treestack[ntreestack++] = newtree; - - // c_ID[i] = vector from per-atom array - - } else if (nbracket == 1 && compute->peratom_flag && - compute->size_peratom_cols > 0) { - - if (tree == nullptr) - print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar); - if (treetype == VECTOR) - print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar); - if (index1 > compute->size_peratom_cols) - print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; - } - - auto newtree = new Tree(); - newtree->type = ATOMARRAY; - if (compute->array_atom) - newtree->array = &compute->array_atom[0][index1-1]; - newtree->nstride = compute->size_peratom_cols; - treestack[ntreestack++] = newtree; - } else if (nbracket == 1 && compute->local_flag) { - print_var_error(FLERR,"Cannot access local data via indexing",ivar); - } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + // ---------------- // fix // ---------------- From 91d826a5d660e2e00f5a016dfb8db189cbde1985 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 18 Aug 2023 09:34:46 -0600 Subject: [PATCH 06/22] changed compute section of variable formulas --- src/variable.cpp | 62 +++++++++++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/src/variable.cpp b/src/variable.cpp index 5013f3ce550..a41c12d1114 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1626,6 +1626,22 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!compute->vector_flag) print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (compute->size_vector == 0) + print_var_error(FLERR,"Variable formula compute vector is zero length",ivar); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= Compute::INVOKED_VECTOR; + } + + auto newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = compute->vector; + newtree->nvector = compute->size_vector; + newtree->nstride = 1; + treestack[ntreestack++] = newtree; // c_ID[i] = vector from global array @@ -1633,6 +1649,24 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!compute->array_flag) print_var_error(FLERR,"Mismatched compute in variable formula",ivar); + if (compute->size_array_rows == 0) + print_var_error(FLERR,"Variable formula compute array is zero length",ivar); + if (index1 > compute->size_array_cols) + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); + if (!compute->is_initialized()) + print_var_error(FLERR,"Variable formula compute cannot be invoked before " + "initialization by a run",ivar); + if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { + compute->compute_array(); + compute->invoked_flag |= Compute::INVOKED_ARRAY; + } + + auto newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = &compute->array[0][index1-1]; + newtree->nvector = compute->size_array_rows; + newtree->nstride = compute->size_array_cols; + treestack[ntreestack++] = newtree; // no other possibilities for vector-style variable, so error @@ -1695,38 +1729,12 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); } - - if (tree == nullptr) - print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar); - if (treetype == ATOM) - print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar); - if (compute->size_array_rows == 0) - print_var_error(FLERR,"Variable formula compute array is zero length",ivar); - if (!compute->is_initialized()) - print_var_error(FLERR,"Variable formula compute cannot be invoked before " - "initialization by a run",ivar); - if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { - compute->compute_array(); - compute->invoked_flag |= Compute::INVOKED_ARRAY; - } - - auto newtree = new Tree(); - newtree->type = VECTORARRAY; - newtree->array = &compute->array[0][index1-1]; - newtree->nvector = compute->size_array_rows; - newtree->nstride = compute->size_array_cols; - treestack[ntreestack++] = newtree; - - - - - - + // ---------------- // fix // ---------------- From 6e1529ddff9068529f70e8aa29f16c0e49df36a6 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 18 Aug 2023 13:18:50 -0600 Subject: [PATCH 07/22] finish changes to variables --- doc/src/variable.rst | 95 +++++---- src/variable.cpp | 497 +++++++++++++++++++++++-------------------- 2 files changed, 312 insertions(+), 280 deletions(-) diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 38e423b6328..4541de5fa22 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -1174,19 +1174,19 @@ The various allowed compute references in the variable formulas for equal-, vector-, and atom-style variables are listed in the following table: -+--------+------------+--------------------------------------------+ -| equal | c_ID | global scalar | -| equal | c_ID[I] | element of global vector | -| equal | c_ID[I][J] | element of global array | -| equal | C_ID[I] | element of per-atom vector, I = ID of atom | -| equal | C_ID{i}[J] | element of per-atom array, I = ID of atom | -+--------+------------+--------------------------------------------| -| vector | c_ID | global vector | -| vector | c_ID[I] | column of global array | ----------+------------+--------------------------------------------+ -| atom | c_ID | per-atom vector | -| atom | c_ID[I] | column of per-atom array | -+--------+------------+--------------------------------------------+ ++--------+------------+------------------------------------------+ +| equal | c_ID | global scalar | +| equal | c_ID[I] | element of global vector | +| equal | c_ID[I][J] | element of global array | +| equal | C_ID[I] | element of per-atom vector (I = atom ID) | +| equal | C_ID{i}[J] | element of per-atom array (I = atom ID) | ++--------+------------+------------------------------------------+ +| vector | c_ID | global vector | +| vector | c_ID[I] | column of global array | +---------+------------+------------------------------------------+ +| atom | c_ID | per-atom vector | +| atom | c_ID[I] | column of per-atom array | ++--------+------------+------------------------------------------+ Note that if an equal-style variable formula wishes to access per-atom data from a compute, it must use capital "C" as the ID prefix and not @@ -1235,23 +1235,22 @@ Atom-style variables can use scalar values (same as for equal-style varaibles), or per-atom vectors of values. The latter can also be a column of a per-atom array. -The various allowed fix references in the variable formulas for -equal-, vector-, and atom-style variables are listed in the following -table: - -+--------+------------+--------------------------------------------+ -| equal | f_ID | global scalar | -| equal | f_ID[I] | element of global vector | -| equal | f_ID[I][J] | element of global array | -| equal | F_ID[I] | element of per-atom vector, I = ID of atom | -| equal | F_ID{i}[J] | element of per-atom array, I = ID of atom | -+--------+------------+--------------------------------------------| -| vector | f_ID | global vector | -| vector | f_ID[I] | column of global array | ----------+------------+--------------------------------------------+ -| atom | f_ID | per-atom vector | -| atom | f_ID[I] | column of per-atom array | -+--------+------------+--------------------------------------------+ +The allowed fix references in variable formulas for equal-, vector-, +and atom-style variables are listed in the following table: + ++--------+------------+------------------------------------------+ +| equal | f_ID | global scalar | +| equal | f_ID[I] | element of global vector | +| equal | f_ID[I][J] | element of global array | +| equal | F_ID[I] | element of per-atom vector (I = atom ID) | +| equal | F_ID{i}[J] | element of per-atom array (I = atom ID) | ++--------+------------+------------------------------------------+ +| vector | f_ID | global vector | +| vector | f_ID[I] | column of global array | +---------+------------+------------------------------------------+ +| atom | f_ID | per-atom vector | +| atom | f_ID[I] | column of per-atom array | ++--------+------------+------------------------------------------+ Note that if an equal-style variable formula wishes to access per-atom data from a fix, it must use capital "F" as the ID prefix and not @@ -1312,21 +1311,27 @@ including other atom-style or atomfile-style variables. If it uses a vector-style variable, a subscript must be used to access a single value from the vector-style variable. -Examples of different kinds of variable references are as follows. -There is no ambiguity as to what a reference means, since variables -produce only a global scalar or global vector or per-atom vector. - -+------------+----------------------------------------------------------------------+ -| v_name | global scalar from equal-style variable | -+------------+----------------------------------------------------------------------+ -| v_name | global vector from vector-style variable | -+------------+----------------------------------------------------------------------+ -| v_name | per-atom vector from atom-style or atomfile-style variable | -+------------+----------------------------------------------------------------------+ -| v_name[I] | Ith element of a global vector from vector-style variable | -+------------+----------------------------------------------------------------------+ -| v_name[I] | value of atom with ID = I from atom-style or atomfile-style variable | -+------------+----------------------------------------------------------------------+ +The allowed variable references in variable formulas for equal-, +vector-, and atom-style variables are listed in the following table. +Note that there is no ambiguity as to what a reference means, since +referenced variables produce only a global scalar or global vector or +per-atom vector. + ++--------+-----------+-----------------------------------------------------------------------------------+ +| equal | v_name | global scalar from an equal-style variable | +| equal | v_name[I] | element of global vector from a vector-style variable | +| equal | v_name[I] | element of per-atom vector (I = atom ID) from an atom- or atomfile-style variable | ++--------+-----------+-----------------------------------------------------------------------------------+ +| vector | v_name | global scalar from an equal-style variable | +| vector | v_name | global vector from a vector-style variable | +| vector | v_name[I] | element of global vector from a vector-style variable | +| vector | v_name[I] | element of per-atom vector (I = atom ID) from an atom- or atomfile-style variable | ++--------+-----------+-----------------------------------------------------------------------------------+ +| atom | v_name | global scalar from an equal-style variable | +| atom | v_name | per-atom vector from an atom-style or atomfile-style variable | +| atom | v_name[I] | element of global vector from a vector-style variable | +| atom | v_name[I] | element of per-atom vector (I = atom ID) from an atom- or atomfile-style variable | ++--------+-----------+-----------------------------------------------------------------------------------+ For the I index, an integer can be specified or a variable name, specified as v_name, where name is the name of the variable. The diff --git a/src/variable.cpp b/src/variable.cpp index a41c12d1114..ce8f16cd68b 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1509,7 +1509,6 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!compute->scalar_flag) print_var_error(FLERR,"Mismatched compute in variable formula",ivar); - if (!compute->is_initialized()) print_var_error(FLERR,"Variable formula compute cannot be invoked before " "initialization by a run",ivar); @@ -1567,7 +1566,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) else value1 = compute->array[index1-1][index2-1]; argstack[nargstack++] = value1; - // C_ID[i] = scalar element of per-atom vector + // C_ID[i] = scalar element of per-atom vector, note uppercase "C" } else if (!lowercase && nbracket == 1) { @@ -1586,7 +1585,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) peratom2global(1,nullptr,compute->vector_atom,1,index1,tree, treestack,ntreestack,argstack,nargstack); - // C_ID[i][j] = scalar element of per-atom array + // C_ID[i][j] = scalar element of per-atom array, note uppercase "C" } else if (!lowercase && nbracket == 2) { @@ -1728,12 +1727,6 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); } - - - - - - // ---------------- // fix @@ -1753,7 +1746,6 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!fix) print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable formula",word+2),ivar); - // parse zero or one or two trailing brackets // point i beyond last bracket // nbracket = # of bracket pairs @@ -1775,181 +1767,200 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) } } - // f_ID = scalar from global scalar, must be lowercase - - if (nbracket == 0 && fix->scalar_flag && lowercase) { - - if (update->whichflag > 0 && update->ntimestep % fix->global_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - - value1 = fix->compute_scalar(); - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; - - // f_ID[i] = scalar from global vector, must be lowercase - - } else if (nbracket == 1 && fix->vector_flag && lowercase) { + // equal-style variable is being evaluated - if (index1 > fix->size_vector && - fix->size_vector_variable == 0) - print_var_error(FLERR,"Variable formula fix vector is accessed out-of-range",ivar,0); - if (update->whichflag > 0 && update->ntimestep % fix->global_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + if (style[ivar] == EQUAL) { + + // f_ID = scalar from global scalar - value1 = fix->compute_vector(index1-1); - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; + if (lowercase && nbracket == 0) { - // f_ID[i][j] = scalar from global array, must be lowercase + if (!fix->scalar_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - } else if (nbracket == 2 && fix->array_flag && lowercase) { + value1 = fix->compute_scalar(); + argstack[nargstack++] = value1; - if (index1 > fix->size_array_rows && - fix->size_array_rows_variable == 0) - print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); - if (index2 > fix->size_array_cols) - print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); - if (update->whichflag > 0 && update->ntimestep % fix->global_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + // f_ID[i] = scalar from global vector - value1 = fix->compute_array(index1-1,index2-1); - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; + } else if (lowercase && nbracket == 1) { - // f_ID = vector from global vector, lowercase or uppercase + if (!fix->vector_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (index1 > fix->size_vector && + fix->size_vector_variable == 0) + print_var_error(FLERR,"Variable formula fix vector is accessed out-of-range",ivar,0); + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - } else if (nbracket == 0 && fix->vector_flag) { + value1 = fix->compute_vector(index1-1); + argstack[nargstack++] = value1; - if (update->whichflag > 0 && update->ntimestep % fix->global_freq) - print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); - if (tree == nullptr) - print_var_error(FLERR,"Fix global vector in equal-style variable formula",ivar); - if (treetype == ATOM) - print_var_error(FLERR,"Fix global vector in atom-style variable formula",ivar); - if (fix->size_vector == 0) - print_var_error(FLERR,"Variable formula fix vector is zero length",ivar); + // f_ID[i][j] = scalar from global array - int nvec = fix->size_vector; - double *vec; - memory->create(vec,nvec,"variable:values"); - for (int m = 0; m < nvec; m++) - vec[m] = fix->compute_vector(m); + } else if (lowercase && nbracket == 2) { - auto newtree = new Tree(); - newtree->type = VECTORARRAY; - newtree->array = vec; - newtree->nvector = nvec; - newtree->nstride = 1; - newtree->selfalloc = 1; - treestack[ntreestack++] = newtree; + if (!fix->array_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (index1 > fix->size_array_rows && + fix->size_array_rows_variable == 0) + print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); + if (index2 > fix->size_array_cols) + print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + + value1 = fix->compute_array(index1-1,index2-1); + argstack[nargstack++] = value1; - // f_ID[i] = vector from global array, lowercase or uppercase + // F_ID[i] = scalar element of per-atom vector, note uppercase "F" - } else if (nbracket == 1 && fix->array_flag) { + } else if (!lowercase && nbracket == 1) { - if (update->whichflag > 0 && update->ntimestep % fix->global_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - if (tree == nullptr) - print_var_error(FLERR,"Fix global vector in equal-style variable formula",ivar); - if (treetype == ATOM) - print_var_error(FLERR,"Fix global vector in atom-style variable formula",ivar); - if (fix->size_array_rows == 0) - print_var_error(FLERR,"Variable formula fix array is zero length",ivar); + if (!fix->peratom_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (fix->size_peratom_cols) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (update->whichflag > 0 && + update->ntimestep % fix->peratom_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + + peratom2global(1,nullptr,fix->vector_atom,1,index1,tree, + treestack,ntreestack,argstack,nargstack); - int nvec = fix->size_array_rows; - double *vec; - memory->create(vec,nvec,"variable:values"); - for (int m = 0; m < nvec; m++) - vec[m] = fix->compute_array(m,index1-1); + // F_ID[i][j] = scalar element of per-atom array, note uppercase "F" - auto newtree = new Tree(); - newtree->type = VECTORARRAY; - newtree->array = vec; - newtree->nvector = nvec; - newtree->nstride = 1; - newtree->selfalloc = 1; - treestack[ntreestack++] = newtree; + } else if (!lowercase && nbracket == 2) { + + if (!fix->peratom_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (!fix->size_peratom_cols) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (index2 > fix->size_peratom_cols) + print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); + if (update->whichflag > 0 && update->ntimestep % fix->peratom_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + + if (fix->array_atom) + peratom2global(1,nullptr,&fix->array_atom[0][index2-1], + fix->size_peratom_cols,index1, + tree,treestack,ntreestack,argstack,nargstack); + else + peratom2global(1,nullptr,nullptr,fix->size_peratom_cols,index1, + tree,treestack,ntreestack,argstack,nargstack); - // f_ID[i] = scalar from per-atom vector + // no other possibilities for equal-style variable, so error + + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); - } else if (nbracket == 1 && fix->peratom_flag && - fix->size_peratom_cols == 0) { + // vector-style variable is being evaluated - if (update->whichflag > 0 && - update->ntimestep % fix->peratom_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + } else if (style[ivar] == VECTOR) { + + // f_ID = vector from global vector - peratom2global(1,nullptr,fix->vector_atom,1,index1, - tree,treestack,ntreestack,argstack,nargstack); + if (lowercase && nbracket == 0) { - // f_ID[i][j] = scalar from per-atom array + if (!fix->vector_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (fix->size_vector == 0) + print_var_error(FLERR,"Variable formula fix vector is zero length",ivar); + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); + + int nvec = fix->size_vector; + double *vec; + memory->create(vec,nvec,"variable:values"); + for (int m = 0; m < nvec; m++) + vec[m] = fix->compute_vector(m); + + auto newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + newtree->selfalloc = 1; + treestack[ntreestack++] = newtree; + + // f_ID[i] = vector from global array - } else if (nbracket == 2 && fix->peratom_flag && - fix->size_peratom_cols > 0) { + } else if (lowercase && nbracket == 1) { - if (index2 > fix->size_peratom_cols) - print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); - if (update->whichflag > 0 && - update->ntimestep % fix->peratom_freq) - print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + if (!fix->array_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (fix->size_array_rows == 0) + print_var_error(FLERR,"Variable formula fix array is zero length",ivar); + if (index1 > fix->size_array_cols) + print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); + + int nvec = fix->size_array_rows; + double *vec; + memory->create(vec,nvec,"variable:values"); + for (int m = 0; m < nvec; m++) + vec[m] = fix->compute_array(m,index1-1); - if (fix->array_atom) - peratom2global(1,nullptr,&fix->array_atom[0][index2-1],fix->size_peratom_cols,index1, - tree,treestack,ntreestack,argstack,nargstack); - else - peratom2global(1,nullptr,nullptr,fix->size_peratom_cols,index1, - tree,treestack,ntreestack,argstack,nargstack); + auto newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + newtree->selfalloc = 1; + treestack[ntreestack++] = newtree; + + // no other possibilities for vector-style variable, so error + + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); - // f_ID = vector from per-atom vector + // atom-style variable is being evaluated - } else if (nbracket == 0 && fix->peratom_flag && - fix->size_peratom_cols == 0) { + } else if (style[ivar] == ATOM) { + + // f_ID = vector from per-atom vector - if (tree == nullptr) - print_var_error(FLERR,"Per-atom fix in equal-style variable formula",ivar); - if (update->whichflag > 0 && - update->ntimestep % fix->peratom_freq) - print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); + if (lowercase && nbracket == 0) { - auto newtree = new Tree(); - newtree->type = ATOMARRAY; - newtree->array = fix->vector_atom; - newtree->nstride = 1; - treestack[ntreestack++] = newtree; + if (!fix->peratom_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (fix->size_peratom_cols) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (update->whichflag > 0 && update->ntimestep % fix->peratom_freq) + print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); - // f_ID[i] = vector from per-atom array + auto newtree = new Tree(); + newtree->type = ATOMARRAY; + newtree->array = fix->vector_atom; + newtree->nstride = 1; + treestack[ntreestack++] = newtree; + + // f_ID[i] = vector from per-atom array - } else if (nbracket == 1 && fix->peratom_flag && - fix->size_peratom_cols > 0) { + } else if (lowercase && nbracket == 1) { - if (tree == nullptr) - print_var_error(FLERR,"Per-atom fix in equal-style variable formula",ivar); - if (index1 > fix->size_peratom_cols) - print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); - if (update->whichflag > 0 && - update->ntimestep % fix->peratom_freq) - print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); + if (!fix->peratom_flag) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (!fix->size_peratom_cols) + print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + if (index1 > fix->size_peratom_cols) + print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); + if (update->whichflag > 0 && update->ntimestep % fix->peratom_freq) + print_var_error(FLERR,"Fix in variable not computed at compatible time",ivar); - auto newtree = new Tree(); - newtree->type = ATOMARRAY; - if (fix->array_atom) - newtree->array = &fix->array_atom[0][index1-1]; - newtree->nstride = fix->size_peratom_cols; - treestack[ntreestack++] = newtree; + auto newtree = new Tree(); + newtree->type = ATOMARRAY; + newtree->array = nullptr; + if (fix->array_atom) + newtree->array = &fix->array_atom[0][index1-1]; + newtree->nstride = fix->size_peratom_cols; + treestack[ntreestack++] = newtree; - } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + // no other possibilities for atom-style variable, so error + + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); + } // ---------------- // variable @@ -1979,124 +1990,140 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) i = ptr-str+1; } - // v_name = scalar from internal-style variable - // access value directly + // vname with no bracket - if (nbracket == 0 && style[ivar] == INTERNAL) { + if (nbracket == 0) { - value1 = dvalue[ivar]; - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value1; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value1; + // scalar from internal-style variable + // access value directly + + if (style[ivar] = INTERNAL) { - // v_name = scalar from non atom/atomfile & non vector-style variable - // access value via retrieve() + value1 = dvalue[ivar]; + if (tree) { + auto newtree = new Tree(); + newtree->type = VALUE; + newtree->value = value1; + treestack[ntreestack++] = newtree; + } else argstack[nargstack++] = value1; + + // scalar from any style variable except VECTOR, ATOM, ATOMFILE + // access value via retrieve() - } else if (nbracket == 0 && style[ivar] != ATOM && - style[ivar] != ATOMFILE && style[ivar] != VECTOR) { + } else if (style[ivar] != ATOM && style[ivar] != ATOMFILE && style[ivar] != VECTOR) { - char *var = retrieve(word+2); - if (var == nullptr) - print_var_error(FLERR,"Invalid variable evaluation in variable formula",ivar); - if (utils::is_double(var)) { + char *var = retrieve(word+2); + if (var == nullptr) + print_var_error(FLERR,"Invalid variable evaluation in variable formula",ivar); + if (!utils::is_double(var)) + print_var_error(FLERR,"Non-numeric variable value in variable formula",ivar); if (tree) { auto newtree = new Tree(); newtree->type = VALUE; newtree->value = atof(var); treestack[ntreestack++] = newtree; } else argstack[nargstack++] = atof(var); - } else print_var_error(FLERR,"Non-numeric variable value in variable formula",ivar); - // v_name = per-atom vector from atom-style variable - // evaluate the atom-style variable as newtree + // vector from vector-style variable + // evaluate the vector-style variable, put result in newtree - } else if (nbracket == 0 && style[ivar] == ATOM) { + } else if (style[ivar] == VECTOR) { - if (tree == nullptr) - print_var_error(FLERR,"Atom-style variable in equal-style variable formula",ivar); - if (treetype == VECTOR) - print_var_error(FLERR,"Atom-style variable in vector-style variable formula",ivar); + if (tree == nullptr) + print_var_error(FLERR,"Vector-style variable in equal-style variable formula",ivar); + if (treetype == ATOM) + print_var_error(FLERR,"Vector-style variable in atom-style variable formula",ivar); - Tree *newtree = nullptr; - evaluate(data[ivar][0],&newtree,ivar); - treestack[ntreestack++] = newtree; + double *vec; + int nvec = compute_vector(ivar,&vec); - // v_name = per-atom vector from atomfile-style variable + auto newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + treestack[ntreestack++] = newtree; - } else if (nbracket == 0 && style[ivar] == ATOMFILE) { + // vector from atom-style variable + // evaluate the atom-style variable as newtree - if (tree == nullptr) - print_var_error(FLERR,"Atomfile-style variable in equal-style variable formula",ivar); - if (treetype == VECTOR) - print_var_error(FLERR,"Atomfile-style variable in vector-style variable formula",ivar); + } else if (style[ivar] == ATOM) { - auto newtree = new Tree(); - newtree->type = ATOMARRAY; - newtree->array = reader[ivar]->fixstore->vstore; - newtree->nstride = 1; - treestack[ntreestack++] = newtree; + if (tree == nullptr) + print_var_error(FLERR,"Atom-style variable in equal-style variable formula",ivar); + if (treetype == VECTOR) + print_var_error(FLERR,"Atom-style variable in vector-style variable formula",ivar); + + Tree *newtree = nullptr; + evaluate(data[ivar][0],&newtree,ivar); + treestack[ntreestack++] = newtree; + + // vector from atomfile-style variable + // point to the values in FixStore instance + + } else if (style[ivar] == ATOMFILE) { - // v_name = vector from vector-style variable - // evaluate the vector-style variable, put result in newtree + if (tree == nullptr) + print_var_error(FLERR,"Atomfile-style variable in equal-style variable formula",ivar); + if (treetype == VECTOR) + print_var_error(FLERR,"Atomfile-style variable in vector-style variable formula",ivar); + + auto newtree = new Tree(); + newtree->type = ATOMARRAY; + newtree->array = reader[ivar]->fixstore->vstore; + newtree->nstride = 1; + treestack[ntreestack++] = newtree; - } else if (nbracket == 0 && style[ivar] == VECTOR) { + // no other possibilities for variable with no bracket - if (tree == nullptr) - print_var_error(FLERR,"Vector-style variable in equal-style variable formula",ivar); - if (treetype == ATOM) - print_var_error(FLERR,"Vector-style variable in atom-style variable formula",ivar); + } else print_var_error(FLERR,"Mismatched variable in variable formula",ivar); - double *vec; - int nvec = compute_vector(ivar,&vec); + // vname[i] with one bracket - auto newtree = new Tree(); - newtree->type = VECTORARRAY; - newtree->array = vec; - newtree->nvector = nvec; - newtree->nstride = 1; - treestack[ntreestack++] = newtree; + } else if (nbracket == 1) { - // v_name[N] = scalar from atom-style variable - // compute the per-atom variable in result - // use peratom2global to extract single value from result + // scalar from vector-style variable + // compute the vector-style variable, extract single value - } else if (nbracket && style[ivar] == ATOM) { + if (style[ivar] == VECTOR) { - double *result; - memory->create(result,atom->nlocal,"variable:result"); - compute_atom(ivar,0,result,1,0); - peratom2global(1,nullptr,result,1,index,tree,treestack,ntreestack,argstack,nargstack); - memory->destroy(result); + double *vec; + int nvec = compute_vector(ivar,&vec); + if (index <= 0 || index > nvec) + print_var_error(FLERR,"Invalid index into vector-style variable",ivar); + int m = index; // convert from tagint to int - // v_name[N] = scalar from atomfile-style variable + if (tree) { + auto newtree = new Tree(); + newtree->type = VALUE; + newtree->value = vec[m-1]; + treestack[ntreestack++] = newtree; + } else argstack[nargstack++] = vec[m-1]; - } else if (nbracket && style[ivar] == ATOMFILE) { + // scalar from atom-style variable + // compute the per-atom variable in result + // use peratom2global to extract single value from result - peratom2global(1,nullptr,reader[ivar]->fixstore->vstore,1,index, - tree,treestack,ntreestack,argstack,nargstack); + } else if (style[ivar] == ATOM) { - // v_name[N] = scalar from vector-style variable - // compute the vector-style variable, extract single value + double *result; + memory->create(result,atom->nlocal,"variable:result"); + compute_atom(ivar,0,result,1,0); + peratom2global(1,nullptr,result,1,index,tree,treestack,ntreestack,argstack,nargstack); + memory->destroy(result); - } else if (nbracket && style[ivar] == VECTOR) { + // scalar from atomfile-style variable + // use peratom2global to extract single value from FixStore instance - double *vec; - int nvec = compute_vector(ivar,&vec); - if (index <= 0 || index > nvec) - print_var_error(FLERR,"Invalid index into vector-style variable",ivar); - int m = index; // convert from tagint to int + } else if (style[ivar] == ATOMFILE) { - if (tree) { - auto newtree = new Tree(); - newtree->type = VALUE; - newtree->value = vec[m-1]; - treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = vec[m-1]; + peratom2global(1,nullptr,reader[ivar]->fixstore->vstore,1,index, + tree,treestack,ntreestack,argstack,nargstack); - } else print_var_error(FLERR,"Mismatched variable in variable formula",ivar); + // no other possibilities for variable with one bracket + + } else print_var_error(FLERR,"Mismatched variable in variable formula",ivar); + } // ---------------- // math/group/special/labelmap function or atom value/vector or From c6233547a531cee52049780592e8c3589960633d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 21 Aug 2023 09:39:00 -0600 Subject: [PATCH 08/22] update compute and fix doc pages for new generality --- doc/src/compute.rst | 124 ++++++++++++++++++++++++-------------------- doc/src/fix.rst | 97 +++++++++++++++++++--------------- 2 files changed, 122 insertions(+), 99 deletions(-) diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 226dc6373bc..2780cac3686 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -27,58 +27,62 @@ Examples Description """"""""""" -Define a computation that will be performed on a group of atoms. -Quantities calculated by a compute are instantaneous values, meaning -they are calculated from information about atoms on the current -timestep or iteration, though a compute may internally store some -information about a previous state of the system. Defining a compute -does not perform a computation. Instead computes are invoked by other -LAMMPS commands as needed (e.g., to calculate a temperature needed for -a thermostat fix or to generate thermodynamic or dump file output). -See the :doc:`Howto output ` page for a summary of -various LAMMPS output options, many of which involve computes. +Define a diagnostic computation that will be performed on a group of +atoms. Quantities calculated by a compute are instantaneous values, +meaning they are calculated from information about atoms on the +current timestep or iteration, though internally a compute may store +some information about a previous state of the system. Defining a +compute does not perform the computation. Instead computes are +invoked by other LAMMPS commands as needed (e.g., to calculate a +temperature needed for a thermostat fix or to generate thermodynamic +or dump file output). See the :doc:`Howto output ` page +for a summary of various LAMMPS output options, many of which involve +computes. The ID of a compute can only contain alphanumeric characters and underscores. ---------- -Computes calculate and store any of four styles of quantities: global, -per-atom, local, or per-grid. A global quantity is one or more -system-wide values, e.g. the temperature of the system. A per-atom -quantity is one or more values per atom, e.g. the kinetic energy of -each atom. Per-atom values are set to 0.0 for atoms not in the -specified compute group. Local quantities are calculated by each -processor based on the atoms it owns, but there may be zero or more -per atom, e.g. a list of bond distances. Per-grid quantities are -calculated on a regular 2d or 3d grid which overlays a 2d or 3d -simulation domain. The grid points and the data they store are -distributed across processors; each processor owns the grid points -which fall within its subdomain. - -Computes that produce per-atom quantities have the word "atom" at the -end of their style, e.g. *ke/atom*\ . Computes that produce local -quantities have the word "local" at the end of their style, -e.g. *bond/local*\ . Computes that produce per-grid quantities have -the word "grid" at the end of their style, e.g. *property/grid*\ . -Styles with neither "atom" or "local" or "grid" at the end of their -style name produce global quantities. - -Note that a single compute typically produces either global or -per-atom or local or per-grid values. It does not compute both global -and per-atom values. It can produce local values or per-grid values -in tandem with global or per-atom quantities. The compute doc page -will explain the details. - -Global, per-atom, local, and per-grid quantities come in three kinds: -a single scalar value, a vector of values, or a 2d array of values. -The doc page for each compute describes the style and kind of values -it produces, e.g. a per-atom vector. Some computes produce more than -one kind of a single style, e.g. a global scalar and a global vector. - -When a compute quantity is accessed, as in many of the output commands -discussed below, it can be referenced via the following bracket -notation, where ID is the ID of the compute: +Computes calculate and store any of four *styles* of quantities: +global, per-atom, local, or per-grid. + +A global quantity is one or more system-wide values, e.g. the +temperature of the system. A per-atom quantity is one or more values +per atom, e.g. the kinetic energy of each atom. Per-atom values are +set to 0.0 for atoms not in the specified compute group. Local +quantities are calculated by each processor based on the atoms it +owns, but there may be zero or more per atom, e.g. a list of bond +distances. Per-grid quantities are calculated on a regular 2d or 3d +grid which overlays a 2d or 3d simulation domain. The grid points and +the data they store are distributed across processors; each processor +owns the grid points which fall within its subdomain. + +As a general rule of thumb, computes that produce per-atom quantities +have the word "atom" at the end of their style, e.g. *ke/atom*\ . +Computes that produce local quantities have the word "local" at the +end of their style, e.g. *bond/local*\ . Computes that produce +per-grid quantities have the word "grid" at the end of their style, +e.g. *property/grid*\ . And styles with neither "atom" or "local" or +"grid" at the end of their style name produce global quantities. + +Global, per-atom, local, and per-grid quantities can also be of three +*kinds*: a single scalar value (global only), a vector of values, or a +2d array of values. For per-atom, local, and per-grid quantities, a +"vector" means a single value for each atom, each local entity +(e.g. bond), or grid cell. Likewise an "array", means multiple values +for each atom, each local entity, or each grid cell. + +Note that a single compute can produce any combination of global, +per-atom, local, or per-grid values. Likewise it can prouduce any +combination of scalar, vector, or array output for each style. The +exception is that for per-atom, local, and per-grid output, either a +vector or array can be produced, but not both. The doc page for each +compute explains the values it produces. + +When a compute output is accessed by another input script command it +is referenced via the following bracket notation, where ID is the ID +of the compute: +-------------+--------------------------------------------+ | c_ID | entire scalar, vector, or array | @@ -89,17 +93,23 @@ notation, where ID is the ID of the compute: +-------------+--------------------------------------------+ In other words, using one bracket reduces the dimension of the -quantity once (vector :math:`\to` scalar, array :math:`\to` vector). Using two -brackets reduces the dimension twice (array :math:`\to` scalar). Thus a -command that uses scalar compute values as input can also process elements of a -vector or array. - -Note that commands and :doc:`variables ` which use compute -quantities typically do not allow for all kinds (e.g., a command may -require a vector of values, not a scalar). This means there is no -ambiguity about referring to a compute quantity as c_ID even if it -produces, for example, both a scalar and vector. The doc pages for -various commands explain the details. +quantity once (vector :math:`\to` scalar, array :math:`\to` vector). +Using two brackets reduces the dimension twice (array :math:`\to` +scalar). Thus, for example, a command that uses global scalar compute +values as input can also process elements of a vector or array. +Depending on the command, this can either be done directly using the +syntax in the table, or by first defining a :doc:`variable ` +of the appropriate style to store the quantity, then using the +variable as an input to the command. + +Note that commands and :doc:`variables ` which take compute +outputs as input typically do not allow for all styles and kinds of +data (e.g., a command may require global but not per-atom values, or +it may require a vector of values, not a scalar). This means there is +typically no ambiguity about referring to a compute output as c_ID +even if it produces, for example, both a scalar and vector. The doc +pages for various commands explain the details, including how any +ambiguities are resolved. ---------- diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 09fc05d5007..a879a45e05b 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -77,35 +77,44 @@ for individual fixes for info on which ones can be restarted. ---------- -Some fixes calculate one or more of four styles of quantities: global, -per-atom, local, or per-grid, which can be used by other commands or -output as described below. A global quantity is one or more -system-wide values, e.g. the energy of a wall interacting with -particles. A per-atom quantity is one or more values per atom, -e.g. the displacement vector for each atom since time 0. Per-atom -values are set to 0.0 for atoms not in the specified fix group. Local -quantities are calculated by each processor based on the atoms it -owns, but there may be zero or more per atoms. Per-grid quantities -are calculated on a regular 2d or 3d grid which overlays a 2d or 3d -simulation domain. The grid points and the data they store are -distributed across processors; each processor owns the grid points -which fall within its subdomain. - -Note that a single fix typically produces either global or per-atom or -local or per-grid values (or none at all). It does not produce both -global and per-atom. It can produce local or per-grid values in -tandem with global or per-atom values. The fix doc page will explain -the details. - -Global, per-atom, local, and per-grid quantities come in three kinds: -a single scalar value, a vector of values, or a 2d array of values. -The doc page for each fix describes the style and kind of values it -produces, e.g. a per-atom vector. Some fixes produce more than one -kind of a single style, e.g. a global scalar and a global vector. - -When a fix quantity is accessed, as in many of the output commands -discussed below, it can be referenced via the following bracket -notation, where ID is the ID of the fix: +Some fixes calculate and store any of four *styles* of quantities: +global, per-atom, local, or per-grid. + +A global quantity is one or more system-wide values, e.g. the energy +of a wall interacting with particles. A per-atom quantity is one or +more values per atom, e.g. the original coordinates of each atom at +time 0. Per-atom values are set to 0.0 for atoms not in the specified +fix group. Local quantities are calculated by each processor based on +the atoms it owns, but there may be zero or more per atom, e.g. values +for each bond. Per-grid quantities are calculated on a regular 2d or +3d grid which overlays a 2d or 3d simulation domain. The grid points +and the data they store are distributed across processors; each +processor owns the grid points which fall within its subdomain. + +As a general rule of thumb, fixes that produce per-atom quantities +have the word "atom" at the end of their style, e.g. *ave/atom*\ . +Fixes that produce local quantities have the word "local" at the end +of their style, e.g. *store/local*\ . Fixes that produce per-grid +quantities have the word "grid" at the end of their style, +e.g. *ave/grid*\ . + +Global, per-atom, local, and per-grid quantities can also be of three +*kinds*: a single scalar value (global only), a vector of values, or a +2d array of values. For per-atom, local, and per-grid quantities, a +"vector" means a single value for each atom, each local entity +(e.g. bond), or grid cell. Likewise an "array", means multiple values +for each atom, each local entity, or each grid cell. + +Note that a single fix can produce any combination of global, +per-atom, local, or per-grid values. Likewise it can prouduce any +combination of scalar, vector, or array output for each style. The +exception is that for per-atom, local, and per-grid output, either a +vector or array can be produced, but not both. The doc page for each +fix explains the values it produces, if any. + +When a fix output is accessed by another input script command it is +referenced via the following bracket notation, where ID is the ID of +the fix: +-------------+--------------------------------------------+ | f_ID | entire scalar, vector, or array | @@ -116,19 +125,23 @@ notation, where ID is the ID of the fix: +-------------+--------------------------------------------+ In other words, using one bracket reduces the dimension of the -quantity once (vector :math:`\to` scalar, array :math:`\to` vector). Using two -brackets reduces the dimension twice (array :math:`\to` scalar). Thus, a -command that uses scalar fix values as input can also process elements of a -vector or array. - -Note that commands and :doc:`variables ` that use fix -quantities typically do not allow for all kinds (e.g., a command may -require a vector of values, not a scalar), and even if they do, the context -in which they are called can be used to resolve which output is being -requested. This means there is no -ambiguity about referring to a fix quantity as f_ID even if it -produces, for example, both a scalar and vector. The doc pages for -various commands explain the details. +quantity once (vector :math:`\to` scalar, array :math:`\to` vector). +Using two brackets reduces the dimension twice (array :math:`\to` +scalar). Thus, for example, a command that uses global scalar fix +values as input can also process elements of a vector or array. +Depending on the command, this can either be done directly using the +syntax in the table, or by first defining a :doc:`variable ` +of the appropriate style to store the quantity, then using the +variable as an input to the command. + +Note that commands and :doc:`variables ` which take fix +outputs as input typically do not allow for all styles and kinds of +data (e.g., a command may require global but not per-atom values, or +it may require a vector of values, not a scalar). This means there is +typically no ambiguity about referring to a fix output as c_ID even if +it produces, for example, both a scalar and vector. The doc pages for +various commands explain the details, including how any ambiguities +are resolved. ---------- From ab2b83f65430892d1f45a928248f01507b7ddfed Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 21 Aug 2023 10:54:42 -0600 Subject: [PATCH 09/22] clarify doc for fix ave/histo command --- doc/src/compute_reduce.rst | 26 +++--- doc/src/fix_ave_histo.rst | 42 ++++++---- src/compute_reduce_region.cpp | 149 +++++++++++----------------------- src/fix_ave_histo.cpp | 2 +- 4 files changed, 88 insertions(+), 131 deletions(-) diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 31591d4419e..4692e161b49 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -63,9 +63,11 @@ Description """"""""""" Define a calculation that "reduces" one or more vector inputs into -scalar values, one per listed input. The inputs can be per-atom or -local quantities and must all be the same kind (per-atom or local); -see discussion of the optional *inputs* keyword below. +scalar values, one per listed input. For the compute reduce command, +the inputs can be either per-atom or local quantities and must all be +of the same kind (per-atom or local); see discussion of the optional +*inputs* keyword below. The compute reduce/region command can only be +used with per-atom inputs. Atom attributes are per-atom quantities, :doc:`computes ` and :doc:`fixes ` can generate either per-atom or local quantities, @@ -92,13 +94,13 @@ values. Each listed input is operated on independently. For per-atom inputs, the group specified with this command means only atoms within the -group contribute to the result. For per-atom inputs, if the compute -reduce/region command is used, the atoms must also currently be within -the region. Note that an input that produces per-atom quantities may -define its own group which affects the quantities it returns. For -example, if a compute is used as an input which generates a per-atom -vector, it will generate values of 0.0 for atoms that are not in the -group specified for that compute. +group contribute to the result. Likewise for per-atom inputs, if the +compute reduce/region command is used, the atoms must also currently +be within the region. Note that an input that produces per-atom +quantities may define its own group which affects the quantities it +returns. For example, if a compute is used as an input which +generates a per-atom vector, it will generate values of 0.0 for atoms +that are not in the group specified for that compute. Each listed input can be an atom attribute (position, velocity, force component) or can be the result of a :doc:`compute ` or @@ -246,7 +248,9 @@ the quantities being reduced are in. Restrictions """""""""""" - none + +As noted above, the compute reduce/region command can only be used +with per-atom inputs. Related commands """""""""""""""" diff --git a/doc/src/fix_ave_histo.rst b/doc/src/fix_ave_histo.rst index 8bb66f06153..31e5476f9ef 100644 --- a/doc/src/fix_ave_histo.rst +++ b/doc/src/fix_ave_histo.rst @@ -79,9 +79,10 @@ Description Use one or more values as inputs every few timesteps to create a single histogram. The histogram can then be averaged over longer -timescales. The resulting histogram can be used by other :doc:`output commands `, and can also be written to a file. The -fix ave/histo/weight command has identical syntax to fix ave/histo, -except that exactly two values must be specified. See details below. +timescales. The resulting histogram can be used by other :doc:`output +commands `, and can also be written to a file. The fix +ave/histo/weight command has identical syntax to fix ave/histo, except +that exactly two values must be specified. See details below. The group specified with this command is ignored for global and local input values. For per-atom input values, only atoms in the group @@ -96,14 +97,18 @@ different ways; see the discussion of the *beyond* keyword below. Each input value can be an atom attribute (position, velocity, force component) or can be the result of a :doc:`compute ` or -:doc:`fix ` or the evaluation of an equal-style or vector-style or -atom-style :doc:`variable `. The set of input values can be -either all global, all per-atom, or all local quantities. Inputs of -different kinds (e.g. global and per-atom) cannot be mixed. Atom -attributes are per-atom vector values. See the page for -individual "compute" and "fix" commands to see what kinds of -quantities they generate. See the optional *kind* keyword below for -how to force the fix ave/histo command to disambiguate if necessary. +:doc:`fix ` or the evaluation of an equal-style or vector-style +or atom-style :doc:`variable `. The set of input values can +be either all global, all per-atom, or all local quantities. Inputs +of different kinds (e.g. global and per-atom) cannot be mixed. Atom +attributes are per-atom vector values. See the page for individual +"compute" and "fix" commands to see what kinds of quantities they +generate. + +Note that a compute or fix can produce multiple kinds of data (global, +per-atom, local). If LAMMPS cannot unambiguosly determine which kind +of data to use, the optional *kind* keyword discussed below can force +the desired disambiguation. Note that the output of this command is a single histogram for all input values combined together, not one histogram per input value. @@ -258,13 +263,14 @@ keyword is set to *vector*, then all input values must be global or per-atom or local vectors, or columns of global or per-atom or local arrays. -The *kind* keyword only needs to be set if a compute or fix produces -more than one kind of output (global, per-atom, local). If this is -not the case, then LAMMPS will determine what kind of input is -provided and whether all the input arguments are consistent. If a -compute or fix produces more than one kind of output, the *kind* -keyword should be used to specify which output will be used. The -remaining input arguments must still be consistent. +The *kind* keyword only needs to be used if any of the specfied input +computes or fixes produce more than one kind of output (global, +per-atom, local). If not, LAMMPS will determine the kind of data all +the inputs produce and verify it is all the same kind. If not, an +error will be triggered. If a compute or fix produces more than one +kind of output, the *kind* keyword should be used to specify which +output will be used. The other input arguments must still be +consistent. The *beyond* keyword determines how input values that fall outside the *lo* to *hi* bounds are treated. Values such that *lo* :math:`\le` value diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp index 2f5a3de6751..d0a32b8adff 100644 --- a/src/compute_reduce_region.cpp +++ b/src/compute_reduce_region.cpp @@ -35,13 +35,15 @@ static constexpr double BIG = 1.0e20; ComputeReduceRegion::ComputeReduceRegion(LAMMPS *lmp, int narg, char **arg) : ComputeReduce(lmp, narg, arg) { + if (input_mode == LOCAL) + error->all(FLERR,"Compute reduce/region cannot use local data as input"); } /* ---------------------------------------------------------------------- calculate reduced value for one input M and return it if flag = -1: sum/min/max/ave all values in vector - for per-atom quantities, limit to atoms in group and region + limit to atoms in group and region if mode = MIN or MAX, also set index to which vector value wins if flag >= 0: simply return vector[flag] ------------------------------------------------------------------------- */ @@ -59,6 +61,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) // initialization in case it has not yet been run, e.g. when // the compute was invoked right after it has been created + if ((val.which == ArgInfo::COMPUTE) || (val.which == ArgInfo::FIX)) { if (val.val.c == nullptr) init(); } @@ -99,52 +102,29 @@ double ComputeReduceRegion::compute_one(int m, int flag) // invoke compute if not previously invoked } else if (val.which == ArgInfo::COMPUTE) { - if (input_mode == PERATOM) { - if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { - val.val.c->compute_peratom(); - val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; - } - - if (aidx == 0) { - double *compute_vector = val.val.c->vector_atom; - if (flag < 0) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) - combine(one, compute_vector[i], i); - } else - one = compute_vector[flag]; - } else { - double **compute_array = val.val.c->array_atom; - int aidxm1 = aidx - 1; - if (flag < 0) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) - combine(one, compute_array[i][aidxm1], i); - } else - one = compute_array[flag][aidxm1]; - } - - } else if (input_mode == LOCAL) { - if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { - val.val.c->compute_local(); - val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; - } - - if (aidx == 0) { - double *compute_vector = val.val.c->vector_local; - if (flag < 0) - for (int i = 0; i < val.val.c->size_local_rows; i++) combine(one, compute_vector[i], i); - else - one = compute_vector[flag]; - } else { - double **compute_array = val.val.c->array_local; - int aidxm1 = aidx - 1; - if (flag < 0) - for (int i = 0; i < val.val.c->size_local_rows; i++) + + if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { + val.val.c->compute_peratom(); + val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; + } + + if (aidx == 0) { + double *compute_vector = val.val.c->vector_atom; + if (flag < 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) + combine(one, compute_vector[i], i); + } else + one = compute_vector[flag]; + } else { + double **compute_array = val.val.c->array_atom; + int aidxm1 = aidx - 1; + if (flag < 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) combine(one, compute_array[i][aidxm1], i); - else - one = compute_array[flag][aidxm1]; - } + } else + one = compute_array[flag][aidxm1]; } // check if fix frequency is a match @@ -153,45 +133,26 @@ double ComputeReduceRegion::compute_one(int m, int flag) if (update->ntimestep % val.val.f->peratom_freq) error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style); - if (input_mode == PERATOM) { - if (aidx == 0) { - double *fix_vector = val.val.f->vector_atom; - if (flag < 0) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) - combine(one, fix_vector[i], i); - } else - one = fix_vector[flag]; - } else { - double **fix_array = val.val.f->array_atom; - int aidxm1 = aidx - 1; - if (flag < 0) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) - combine(one, fix_array[i][aidxm1], i); - } else - one = fix_array[flag][aidxm1]; - } - - } else if (input_mode == LOCAL) { - if (aidx == 0) { - double *fix_vector = val.val.f->vector_local; - if (flag < 0) - for (int i = 0; i < val.val.f->size_local_rows; i++) combine(one, fix_vector[i], i); - else - one = fix_vector[flag]; - } else { - double **fix_array = val.val.f->array_local; - int aidxm1 = aidx - 1; - if (flag < 0) - for (int i = 0; i < val.val.f->size_local_rows; i++) + if (aidx == 0) { + double *fix_vector = val.val.f->vector_atom; + if (flag < 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) + combine(one, fix_vector[i], i); + } else + one = fix_vector[flag]; + } else { + double **fix_array = val.val.f->array_atom; + int aidxm1 = aidx - 1; + if (flag < 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) combine(one, fix_array[i][aidxm1], i); - else - one = fix_array[flag][aidxm1]; - } + } else + one = fix_array[flag][aidxm1]; } - // evaluate atom-style variable + // evaluate atom-style variable } else if (val.which == ArgInfo::VARIABLE) { if (atom->nmax > maxatom) { @@ -220,25 +181,11 @@ bigint ComputeReduceRegion::count(int m) if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) return group->count(igroup, region); - else if (val.which == ArgInfo::COMPUTE) { - if (input_mode == PERATOM) { - return group->count(igroup, region); - } else if (input_mode == LOCAL) { - bigint ncount = val.val.c->size_local_rows; - bigint ncountall; - MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world); - return ncountall; - } - } else if (val.which == ArgInfo::FIX) { - if (input_mode == PERATOM) { - return group->count(igroup, region); - } else if (input_mode == LOCAL) { - bigint ncount = val.val.f->size_local_rows; - bigint ncountall; - MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world); - return ncountall; - } - } else if (val.which == ArgInfo::VARIABLE) + else if (val.which == ArgInfo::COMPUTE) + return group->count(igroup, region); + else if (val.which == ArgInfo::FIX) + return group->count(igroup, region); + else if (val.which == ArgInfo::VARIABLE) return group->count(igroup, region); bigint dummy = 0; diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp index 0a2975bb2e9..4503ad56f4e 100644 --- a/src/fix_ave_histo.cpp +++ b/src/fix_ave_histo.cpp @@ -164,7 +164,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : } // check input args for kind consistency - // all inputs must all be global, per-atom, or local + // inputs must all be all either global, per-atom, or local if (nevery <= 0) error->all(FLERR,"Illegal {} nevery value: {}", mycmd, nevery); From aad232ffc64240a24f5421c46df0a5d84aceff0a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 22 Aug 2023 11:46:57 -0400 Subject: [PATCH 10/22] fix typo --- src/variable.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/variable.cpp b/src/variable.cpp index ce8f16cd68b..f3c987f00cc 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1997,7 +1997,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // scalar from internal-style variable // access value directly - if (style[ivar] = INTERNAL) { + if (style[ivar] == INTERNAL) { value1 = dvalue[ivar]; if (tree) { From 7d9c068da09d8135c2e672f2de21b1ed20f4510d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 22 Aug 2023 11:50:54 -0400 Subject: [PATCH 11/22] whitespace --- src/VORONOI/compute_voronoi_atom.cpp | 12 +++---- src/compute_reduce.cpp | 8 ++--- src/compute_reduce_region.cpp | 2 +- src/variable.cpp | 54 ++++++++++++++-------------- 4 files changed, 38 insertions(+), 38 deletions(-) diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp index eb4f53986f8..b4f1aa30551 100644 --- a/src/VORONOI/compute_voronoi_atom.cpp +++ b/src/VORONOI/compute_voronoi_atom.cpp @@ -391,13 +391,13 @@ void ComputeVoronoi::checkOccupation() int i, j, k; double rx, ry, rz; - + int nlocal = atom->nlocal; int nall = atom->nghost + nlocal; double **x = atom->x; // prepare destination buffer for variable evaluation - + if (atom->nmax > lmax) { memory->destroy(lnext); lmax = atom->nmax; @@ -432,7 +432,7 @@ void ComputeVoronoi::checkOccupation() } // MPI sum occupation - + #ifdef NOTINPLACE memcpy(sendocc, occvec, oldnatoms*sizeof(*occvec)); MPI_Allreduce(sendocc, occvec, oldnatoms, MPI_INT, MPI_SUM, world); @@ -441,7 +441,7 @@ void ComputeVoronoi::checkOccupation() #endif // determine the total number of atoms in this atom's currently occupied cell - + int c; for (i=0; itag[i]; if (mytag > oldmaxtag) @@ -479,7 +479,7 @@ void ComputeVoronoi::checkOccupation() void ComputeVoronoi::loopCells() { // invoke voro++ and fetch results for owned atoms in group - + voronoicell_neighbor c; int i; if (faces_flag) nfaces = 0; diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 8565ddb1c9b..24fdc4a991b 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -40,7 +40,7 @@ enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduceRegion void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) { // r is the already reduced value, n is the new value - + double n = std::fabs(*(double *) in), r = *(double *) inout; double m; @@ -55,7 +55,7 @@ void abs_max(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) void abs_min(void *in, void *inout, int * /*len*/, MPI_Datatype * /*type*/) { // r is the already reduced value, n is the new value - + double n = std::fabs(*(double *) in), r = *(double *) inout; double m; @@ -270,7 +270,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style, val.id); if (val.argindex && val.argindex > val.val.c->size_peratom_cols) error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id); - + } else if (input_mode == LOCAL) { if (!val.val.c->peratom_flag) error->all(FLERR, "Compute {} compute {} does not calculate local values", style, val.id); @@ -295,7 +295,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Compute {} fix {} does not calculate a per-atom array", style, val.id); if (val.argindex && (val.argindex > val.val.f->size_peratom_cols)) error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id); - + } else if (input_mode == LOCAL) { if (!val.val.f->local_flag) error->all(FLERR, "Compute {} fix {} does not calculate local values", style, val.id); diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp index d0a32b8adff..15280af544c 100644 --- a/src/compute_reduce_region.cpp +++ b/src/compute_reduce_region.cpp @@ -61,7 +61,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) // initialization in case it has not yet been run, e.g. when // the compute was invoked right after it has been created - + if ((val.which == ArgInfo::COMPUTE) || (val.which == ArgInfo::FIX)) { if (val.val.c == nullptr) init(); } diff --git a/src/variable.cpp b/src/variable.cpp index f3c987f00cc..264dcf62586 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1502,7 +1502,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // equal-style variable is being evaluated if (style[ivar] == EQUAL) { - + // c_ID = scalar from global scalar if (lowercase && nbracket == 0) { @@ -1588,7 +1588,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // C_ID[i][j] = scalar element of per-atom array, note uppercase "C" } else if (!lowercase && nbracket == 2) { - + if (!compute->peratom_flag) print_var_error(FLERR,"Mismatched compute in variable formula",ivar); if (!compute->size_peratom_cols) @@ -1612,13 +1612,13 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) tree,treestack,ntreestack,argstack,nargstack); // no other possibilities for equal-style variable, so error - + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); // vector-style variable is being evaluated } else if (style[ivar] == VECTOR) { - + // c_ID = vector from global vector if (lowercase && nbracket == 0) { @@ -1641,7 +1641,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->nvector = compute->size_vector; newtree->nstride = 1; treestack[ntreestack++] = newtree; - + // c_ID[i] = vector from global array } else if (lowercase && nbracket == 1) { @@ -1666,15 +1666,15 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->nvector = compute->size_array_rows; newtree->nstride = compute->size_array_cols; treestack[ntreestack++] = newtree; - + // no other possibilities for vector-style variable, so error - + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); // atom-style variable is being evaluated } else if (style[ivar] == ATOM) { - + // c_ID = vector from per-atom vector if (lowercase && nbracket == 0) { @@ -1696,7 +1696,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->array = compute->vector_atom; newtree->nstride = 1; treestack[ntreestack++] = newtree; - + // c_ID[i] = vector from per-atom array } else if (lowercase && nbracket == 1) { @@ -1724,10 +1724,10 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) treestack[ntreestack++] = newtree; // no other possibilities for atom-style variable, so error - + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); } - + // ---------------- // fix // ---------------- @@ -1770,7 +1770,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // equal-style variable is being evaluated if (style[ivar] == EQUAL) { - + // f_ID = scalar from global scalar if (lowercase && nbracket == 0) { @@ -1826,14 +1826,14 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (update->whichflag > 0 && update->ntimestep % fix->peratom_freq) print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - + peratom2global(1,nullptr,fix->vector_atom,1,index1,tree, treestack,ntreestack,argstack,nargstack); // F_ID[i][j] = scalar element of per-atom array, note uppercase "F" } else if (!lowercase && nbracket == 2) { - + if (!fix->peratom_flag) print_var_error(FLERR,"Mismatched fix in variable formula",ivar); if (!fix->size_peratom_cols) @@ -1852,13 +1852,13 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) tree,treestack,ntreestack,argstack,nargstack); // no other possibilities for equal-style variable, so error - + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); // vector-style variable is being evaluated } else if (style[ivar] == VECTOR) { - + // f_ID = vector from global vector if (lowercase && nbracket == 0) { @@ -1875,7 +1875,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) memory->create(vec,nvec,"variable:values"); for (int m = 0; m < nvec; m++) vec[m] = fix->compute_vector(m); - + auto newtree = new Tree(); newtree->type = VECTORARRAY; newtree->array = vec; @@ -1883,7 +1883,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->nstride = 1; newtree->selfalloc = 1; treestack[ntreestack++] = newtree; - + // f_ID[i] = vector from global array } else if (lowercase && nbracket == 1) { @@ -1910,15 +1910,15 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->nstride = 1; newtree->selfalloc = 1; treestack[ntreestack++] = newtree; - + // no other possibilities for vector-style variable, so error - + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); // atom-style variable is being evaluated } else if (style[ivar] == ATOM) { - + // f_ID = vector from per-atom vector if (lowercase && nbracket == 0) { @@ -1935,7 +1935,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) newtree->array = fix->vector_atom; newtree->nstride = 1; treestack[ntreestack++] = newtree; - + // f_ID[i] = vector from per-atom array } else if (lowercase && nbracket == 1) { @@ -1958,7 +1958,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) treestack[ntreestack++] = newtree; // no other possibilities for atom-style variable, so error - + } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar); } @@ -2053,21 +2053,21 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) print_var_error(FLERR,"Atom-style variable in equal-style variable formula",ivar); if (treetype == VECTOR) print_var_error(FLERR,"Atom-style variable in vector-style variable formula",ivar); - + Tree *newtree = nullptr; evaluate(data[ivar][0],&newtree,ivar); treestack[ntreestack++] = newtree; // vector from atomfile-style variable // point to the values in FixStore instance - + } else if (style[ivar] == ATOMFILE) { if (tree == nullptr) print_var_error(FLERR,"Atomfile-style variable in equal-style variable formula",ivar); if (treetype == VECTOR) print_var_error(FLERR,"Atomfile-style variable in vector-style variable formula",ivar); - + auto newtree = new Tree(); newtree->type = ATOMARRAY; newtree->array = reader[ivar]->fixstore->vstore; @@ -2121,7 +2121,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) tree,treestack,ntreestack,argstack,nargstack); // no other possibilities for variable with one bracket - + } else print_var_error(FLERR,"Mismatched variable in variable formula",ivar); } From dd6b847a5c6f88b24904c74f1a84ae4354400cb5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 22 Aug 2023 16:29:14 -0400 Subject: [PATCH 12/22] mention that "peratom" is no longer required and was removed --- doc/src/compute_voronoi_atom.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/src/compute_voronoi_atom.rst b/doc/src/compute_voronoi_atom.rst index 3e67bb6cbfa..37e53863411 100644 --- a/doc/src/compute_voronoi_atom.rst +++ b/doc/src/compute_voronoi_atom.rst @@ -190,6 +190,10 @@ Voro++ software in the src/VORONOI/README file. Output info """"""""""" +.. deprecated:: TBD + + The *peratom* keyword was removed as it is no longer required. + This compute calculates a per-atom array with two columns. In regular dynamic tessellation mode the first column is the Voronoi volume, the second is the neighbor count, as described above (read above for the From ad33a018f48f228f4d4b9f5815c85a93bd25663f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 22 Aug 2023 15:52:47 -0600 Subject: [PATCH 13/22] update variable syntax in several example input scripts --- examples/snap/in.snap.compute | 2 +- examples/snap/in.snap.compute.quadratic | 2 +- examples/voronoi/in.voronoi | 8 ++++---- examples/voronoi/in.voronoi.data | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/snap/in.snap.compute b/examples/snap/in.snap.compute index b0c7314882e..8d2ffe8b969 100644 --- a/examples/snap/in.snap.compute +++ b/examples/snap/in.snap.compute @@ -70,7 +70,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*] # fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector compute vbsum all reduce sum c_vb[*] # fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector -variable db_2_25 equal c_db[2][25] +variable db_2_25 equal C_db[2][25] # set up compute snap generating global array diff --git a/examples/snap/in.snap.compute.quadratic b/examples/snap/in.snap.compute.quadratic index e03d4af3bfe..20d5ed3039f 100644 --- a/examples/snap/in.snap.compute.quadratic +++ b/examples/snap/in.snap.compute.quadratic @@ -70,7 +70,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*] # fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector compute vbsum all reduce sum c_vb[*] # fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector -variable db_2_100 equal c_db[2][100] +variable db_2_100 equal C_db[2][100] # set up compute snap generating global array diff --git a/examples/voronoi/in.voronoi b/examples/voronoi/in.voronoi index 5254969fbd9..79b6c6efec6 100644 --- a/examples/voronoi/in.voronoi +++ b/examples/voronoi/in.voronoi @@ -146,10 +146,10 @@ variable i2 equal 257 compute v1 all voronoi/atom occupation compute r0 all reduce sum c_v1[1] compute r1 all reduce sum c_v1[2] -variable d5a equal c_v1[${i1}][1] -variable d5b equal c_v1[${i2}][1] -variable d5c equal c_v1[${i1}][2] -variable d5d equal c_v1[${i2}][2] +variable d5a equal C_v1[${i1}][1] +variable d5b equal C_v1[${i2}][1] +variable d5c equal C_v1[${i1}][2] +variable d5d equal C_v1[${i2}][2] thermo_style custom c_r0 c_r1 v_d5a v_d5b v_d5c v_d5d run 0 diff --git a/examples/voronoi/in.voronoi.data b/examples/voronoi/in.voronoi.data index 853c2c2bd10..da00b44e093 100644 --- a/examples/voronoi/in.voronoi.data +++ b/examples/voronoi/in.voronoi.data @@ -67,7 +67,7 @@ undump dlocal # local and global quantities, but # not per-atom quantities -compute v2 all voronoi/atom neighbors yes edge_histo 6 peratom no +compute v2 all voronoi/atom neighbors yes edge_histo 6 # write voronoi local quantities to a file @@ -75,7 +75,7 @@ dump d2 all local 1 dump.neighbors2 index c_v2[1] c_v2[2] c_v2[3] # sum up a voronoi local quantity -compute sumarea all reduce sum c_v2[3] +compute sumarea all reduce sum c_v2[3] inputs local # output voronoi global quantities From 71ca6ee47ca24e0db97e0137b2ad293d33e95647 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 22 Aug 2023 16:02:28 -0600 Subject: [PATCH 14/22] fix one more example input script --- examples/voronoi/in.voronoi.data | 7 +------ src/compute_reduce.cpp | 3 +-- src/compute_reduce.h | 2 +- src/compute_reduce_region.cpp | 2 -- 4 files changed, 3 insertions(+), 11 deletions(-) diff --git a/examples/voronoi/in.voronoi.data b/examples/voronoi/in.voronoi.data index da00b44e093..e5d925c4982 100644 --- a/examples/voronoi/in.voronoi.data +++ b/examples/voronoi/in.voronoi.data @@ -63,9 +63,7 @@ undump dlocal # TEST 2: # -# This compute voronoi generates -# local and global quantities, but -# not per-atom quantities +# This compute voronoi generates peratom and local and global quantities compute v2 all voronoi/atom neighbors yes edge_histo 6 @@ -83,6 +81,3 @@ thermo_style custom c_sumarea c_v2[3] c_v2[4] c_v2[5] c_v2[6] c_v2[7] thermo 1 run 0 - - - diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 24fdc4a991b..59834455176 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -31,8 +31,6 @@ using namespace LAMMPS_NS; -enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduceRegion - #define BIG 1.0e20 //---------------------------------------------------------------- @@ -232,6 +230,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Compute {} inputs must be all peratom or all local"); input_mode = LOCAL; } + iarg += 2; } else error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]); } diff --git a/src/compute_reduce.h b/src/compute_reduce.h index f8b652e00c6..64322bc6ac2 100644 --- a/src/compute_reduce.h +++ b/src/compute_reduce.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS { class ComputeReduce : public Compute { public: enum { SUM, SUMSQ, SUMABS, MINN, MAXX, AVE, AVESQ, AVEABS, MINABS, MAXABS }; - enum { PERATOM, LOCAL }; + enum { UNDECIDED, PERATOM, LOCAL }; ComputeReduce(class LAMMPS *, int, char **); ~ComputeReduce() override; diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp index 15280af544c..bd850e902c8 100644 --- a/src/compute_reduce_region.cpp +++ b/src/compute_reduce_region.cpp @@ -26,8 +26,6 @@ using namespace LAMMPS_NS; -enum{UNDECIDED,PERATOM,LOCAL}; // same as in ComputeReduce - static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ From 17dd04b4dea143fe271188b235edbe13e473956e Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 22 Aug 2023 16:22:57 -0600 Subject: [PATCH 15/22] tweak variable doc page --- doc/src/variable.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 4541de5fa22..f1a316da1f9 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -1179,7 +1179,7 @@ table: | equal | c_ID[I] | element of global vector | | equal | c_ID[I][J] | element of global array | | equal | C_ID[I] | element of per-atom vector (I = atom ID) | -| equal | C_ID{i}[J] | element of per-atom array (I = atom ID) | +| equal | C_ID[I][J] | element of per-atom array (I = atom ID) | +--------+------------+------------------------------------------+ | vector | c_ID | global vector | | vector | c_ID[I] | column of global array | @@ -1243,7 +1243,7 @@ and atom-style variables are listed in the following table: | equal | f_ID[I] | element of global vector | | equal | f_ID[I][J] | element of global array | | equal | F_ID[I] | element of per-atom vector (I = atom ID) | -| equal | F_ID{i}[J] | element of per-atom array (I = atom ID) | +| equal | F_ID[I][J] | element of per-atom array (I = atom ID) | +--------+------------+------------------------------------------+ | vector | f_ID | global vector | | vector | f_ID[I] | column of global array | From 3e22eb83555a484388a433b20be14203dcbb4269 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 22 Aug 2023 16:40:25 -0600 Subject: [PATCH 16/22] adjust version date --- src/version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/version.h b/src/version.h index 572a2740536..35780aa785e 100644 --- a/src/version.h +++ b/src/version.h @@ -1,2 +1,2 @@ -#define LAMMPS_VERSION "2 Aug 2023" +#define LAMMPS_VERSION "3 Aug 2023" #define LAMMPS_UPDATE "Development" From e6b98f5942856c11ec34b3206ee59cb175a8dc83 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 23 Aug 2023 09:47:36 -0600 Subject: [PATCH 17/22] fix logic issue in compute reduce --- src/compute_reduce.cpp | 23 ++++++----------------- src/compute_reduce.h | 2 +- 2 files changed, 7 insertions(+), 18 deletions(-) diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 59834455176..5385554f33f 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -133,8 +133,6 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : // parse values - input_mode = UNDECIDED; - values.clear(); nvalues = 0; for (int iarg = 0; iarg < nargnew; ++iarg) { @@ -144,41 +142,32 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : val.val.c = nullptr; if (strcmp(arg[iarg], "x") == 0) { - input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 0; } else if (strcmp(arg[iarg], "y") == 0) { - input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 1; } else if (strcmp(arg[iarg], "z") == 0) { - input_mode = PERATOM; val.which = ArgInfo::X; val.argindex = 2; } else if (strcmp(arg[iarg], "vx") == 0) { - input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 0; } else if (strcmp(arg[iarg], "vy") == 0) { - input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 1; } else if (strcmp(arg[iarg], "vz") == 0) { - input_mode = PERATOM; val.which = ArgInfo::V; val.argindex = 2; } else if (strcmp(arg[iarg], "fx") == 0) { - input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 0; } else if (strcmp(arg[iarg], "fy") == 0) { - input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 1; } else if (strcmp(arg[iarg], "fz") == 0) { - input_mode = PERATOM; val.which = ArgInfo::F; val.argindex = 2; @@ -203,6 +192,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : nvalues = values.size(); replace = new int[nvalues]; for (int i = 0; i < nvalues; ++i) replace[i] = -1; + input_mode = PERATOM; std::string mycmd = "compute "; mycmd += style; @@ -225,11 +215,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg], "inputs") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, mycmd + " inputs", error); if (strcmp(arg[iarg+1], "peratom") == 0) input_mode = PERATOM; - else if (strcmp(arg[iarg+1], "local") == 0) { - if (input_mode == PERATOM) - error->all(FLERR,"Compute {} inputs must be all peratom or all local"); - input_mode = LOCAL; - } + else if (strcmp(arg[iarg+1], "local") == 0) input_mode = LOCAL; iarg += 2; } else error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]); @@ -255,7 +241,10 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : // setup and error check for (auto &val : values) { - if (val.which == ArgInfo::COMPUTE) { + if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) { + if (input_mode == LOCAL) error->all(FLERR,"Compute {} inputs must be all local"); + + } else if (val.which == ArgInfo::COMPUTE) { val.val.c = modify->get_compute_by_id(val.id); if (!val.val.c) error->all(FLERR, "Compute ID {} for compute {} does not exist", val.id, style); diff --git a/src/compute_reduce.h b/src/compute_reduce.h index 64322bc6ac2..f8b652e00c6 100644 --- a/src/compute_reduce.h +++ b/src/compute_reduce.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS { class ComputeReduce : public Compute { public: enum { SUM, SUMSQ, SUMABS, MINN, MAXX, AVE, AVESQ, AVEABS, MINABS, MAXABS }; - enum { UNDECIDED, PERATOM, LOCAL }; + enum { PERATOM, LOCAL }; ComputeReduce(class LAMMPS *, int, char **); ~ComputeReduce() override; From b4e7d5f0b9a2a8ce329dc5ca9fb383b9b1390861 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 23 Aug 2023 20:11:32 -0400 Subject: [PATCH 18/22] fix whitespace (again) --- src/compute_reduce.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 5385554f33f..3feabf2ec3d 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -243,7 +243,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : for (auto &val : values) { if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) { if (input_mode == LOCAL) error->all(FLERR,"Compute {} inputs must be all local"); - + } else if (val.which == ArgInfo::COMPUTE) { val.val.c = modify->get_compute_by_id(val.id); if (!val.val.c) From 6ccccb5d13bbb9f53033ba64c8beb646eba5ae3f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 24 Aug 2023 09:27:17 -0400 Subject: [PATCH 19/22] add versionadded tag to new inputs keyword docs --- doc/src/compute_reduce.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 4692e161b49..6820d2ee041 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -201,6 +201,8 @@ information in this context, the *replace* keywords will extract the atom IDs for the two atoms in the bond of maximum stretch. These atom IDs and the bond stretch will be printed with thermodynamic output. +.. versionadded:: TBD + The *inputs* keyword allows selection of whether all the inputs are per-atom or local quantities. As noted above, all the inputs must be the same kind (per-atom or local). Per-atom is the default setting. From ad1400ac71ba8b6f305275e97817539c04e054c6 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 16 Oct 2023 12:33:21 -0600 Subject: [PATCH 20/22] Fix broken example --- examples/mliap/in.mliap.quadratic.compute | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute index 929dbf38249..cc9ad331b57 100644 --- a/examples/mliap/in.mliap.quadratic.compute +++ b/examples/mliap/in.mliap.quadratic.compute @@ -65,7 +65,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*] # fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector compute vbsum all reduce sum c_vb[*] # fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector -variable db_2_100 equal c_db[2][100] +variable db_2_100 equal C_db[2][100] # test output: 1: total potential energy # 2: xy component of stress tensor From 8c7493d02a9ce2f2f9171ac9a26fd155703f7a5f Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 16 Oct 2023 15:11:37 -0600 Subject: [PATCH 21/22] Fix more broken examples --- examples/mliap/in.mliap.snap.compute | 2 +- examples/snap/in.grid.snap | 2 +- examples/snap/in.grid.tri | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute index 4cfccedbdf7..c49365f55f1 100644 --- a/examples/mliap/in.mliap.snap.compute +++ b/examples/mliap/in.mliap.snap.compute @@ -65,7 +65,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*] # fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector compute vbsum all reduce sum c_vb[*] # fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector -variable db_2_25 equal c_db[2][25] +variable db_2_25 equal C_db[2][25] thermo 100 diff --git a/examples/snap/in.grid.snap b/examples/snap/in.grid.snap index 08c95a004f1..d37b6ffde43 100644 --- a/examples/snap/in.grid.snap +++ b/examples/snap/in.grid.snap @@ -67,7 +67,7 @@ compute mygridlocal all sna/grid/local grid ${ngrid} ${ngrid} ${ngrid} & # define output -variable B5atom equal c_b[2][5] +variable B5atom equal C_b[2][5] variable B5grid equal c_mygrid[8][8] variable rmse_global equal "sqrt( & diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index 5283957eb82..b34c9dba30c 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -87,7 +87,7 @@ compute mygridlocal all sna/grid/local grid ${ngridx} ${ngridy} ${ngridz} & # define output -variable B5atom equal c_b[7][5] +variable B5atom equal C_b[7][5] variable B5grid equal c_mygrid[13][8] # do not compare x,y,z because assignment of ids From dc67f2527061b036da95c8fd151eb25da57abde3 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 16 Oct 2023 15:17:46 -0600 Subject: [PATCH 22/22] Another tweak --- examples/snap/in.grid.snap | 10 +++++----- examples/snap/in.grid.tri | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/snap/in.grid.snap b/examples/snap/in.grid.snap index d37b6ffde43..da48957d977 100644 --- a/examples/snap/in.grid.snap +++ b/examples/snap/in.grid.snap @@ -74,11 +74,11 @@ variable rmse_global equal "sqrt( & (c_mygrid[8][1] - x[2])^2 + & (c_mygrid[8][2] - y[2])^2 + & (c_mygrid[8][3] - z[2])^2 + & - (c_mygrid[8][4] - c_b[2][1])^2 + & - (c_mygrid[8][5] - c_b[2][2])^2 + & - (c_mygrid[8][6] - c_b[2][3])^2 + & - (c_mygrid[8][7] - c_b[2][4])^2 + & - (c_mygrid[8][8] - c_b[2][5])^2 & + (c_mygrid[8][4] - C_b[2][1])^2 + & + (c_mygrid[8][5] - C_b[2][2])^2 + & + (c_mygrid[8][6] - C_b[2][3])^2 + & + (c_mygrid[8][7] - C_b[2][4])^2 + & + (c_mygrid[8][8] - C_b[2][5])^2 & )" thermo_style custom step v_B5atom v_B5grid v_rmse_global diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index b34c9dba30c..95a14f3bb48 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -94,11 +94,11 @@ variable B5grid equal c_mygrid[13][8] # to atoms is not unnique for different processor grids variable rmse_global equal "sqrt( & - (c_mygrid[13][4] - c_b[7][1])^2 + & - (c_mygrid[13][5] - c_b[7][2])^2 + & - (c_mygrid[13][6] - c_b[7][3])^2 + & - (c_mygrid[13][7] - c_b[7][4])^2 + & - (c_mygrid[13][8] - c_b[7][5])^2 & + (c_mygrid[13][4] - C_b[7][1])^2 + & + (c_mygrid[13][5] - C_b[7][2])^2 + & + (c_mygrid[13][6] - C_b[7][3])^2 + & + (c_mygrid[13][7] - C_b[7][4])^2 + & + (c_mygrid[13][8] - C_b[7][5])^2 & )" thermo_style custom step v_B5atom v_B5grid v_rmse_global