diff --git a/doc/src/compute_modify.rst b/doc/src/compute_modify.rst index a68f3c14e1..2e00157d5e 100644 --- a/doc/src/compute_modify.rst +++ b/doc/src/compute_modify.rst @@ -40,21 +40,26 @@ all compute styles support all parameters. The *extra/dof* or *extra* keyword refers to how many degrees-of-freedom are subtracted (typically from 3N) as a normalizing factor in a temperature computation. Only computes that compute a -temperature use this option. The default is 2 or 3 for :doc:`2d or 3d systems ` which is a correction factor for an ensemble -of velocities with zero total linear momentum. For compute -temp/partial, if one or more velocity components are excluded, the -value used for *extra* is scaled accordingly. You can use a negative -number for the *extra* parameter if you need to add -degrees-of-freedom. See the :doc:`compute temp/asphere ` command for an example. +temperature use this option. The default is 2 or 3 for :doc:`2d or 3d +systems ` which is a correction factor for an ensemble of +velocities with zero total linear momentum. For compute temp/partial, +if one or more velocity components are excluded, the value used for +*extra* is scaled accordingly. You can use a negative number for the +*extra* parameter if you need to add degrees-of-freedom. See the +:doc:`compute temp/asphere ` command for an +example. The *dynamic/dof* or *dynamic* keyword determines whether the number of atoms N in the compute group and their associated degrees of freedom are re-computed each time a temperature is computed. Only compute styles that calculate a temperature use this option. By default, N and their DOF are assumed to be constant. If you are -adding atoms or molecules to the system (see the :doc:`fix pour `, :doc:`fix deposit `, and :doc:`fix gcmc ` commands) or expect atoms or molecules to be lost -(e.g. due to exiting the simulation box or via :doc:`fix evaporate `), then this option should be used to -insure the temperature is correctly normalized. +adding atoms or molecules to the system (see the :doc:`fix pour +`, :doc:`fix deposit `, and :doc:`fix gcmc +` commands) or expect atoms or molecules to be lost +(e.g. due to exiting the simulation box or via :doc:`fix evaporate +`), then this option should be used to insure the +temperature is correctly normalized. .. note:: diff --git a/doc/src/compute_property_grid.rst b/doc/src/compute_property_grid.rst index 86c0ffbc46..7471637b48 100644 --- a/doc/src/compute_property_grid.rst +++ b/doc/src/compute_property_grid.rst @@ -8,10 +8,11 @@ Syntax .. parsed-literal:: - compute ID group-ID property/grid input1 input2 ... + compute ID group-ID property/grid Nx Ny Nz input1 input2 ... * ID, group-ID are documented in :doc:`compute ` command * property/grid = style name of this compute command +* Nx, Ny, Nz = grid size in each dimension * input1,etc = one or more attributes .. parsed-literal:: @@ -46,6 +47,12 @@ This compute stores the specified attributes of grids as per-grid data so they can be accessed by other :doc:`output commands ` such as :doc:`dump grid `. +*Nx*, *Ny*, and *Nz* define the size of the grid. For a 2d simulation +*Nz* must be 1. When this compute is used by :doc:`dump grid `, +*to output per-grid values from other computes of fixes, the grid size +*specified for this command must be consistent with the grid sizes +*used by the other commands. + The *id* attribute stores the grid ID for each grid point. For a global grid of size Nx by Ny by Nz (in 3d simulations) the grid IDs range from 1 to Nx*Ny*Nz. They are ordered with the X index of the 3d @@ -83,10 +90,9 @@ Output info This compute calculates a per-grid vector or array depending on the number of input values. The length of the vector or number of array -rows for each processor is the number of grid points it owns. - -For access by other commands, the name of the grid produced -by this command is "grid". The name of its data is "data". +rows (distributed across all processors) is Nx * Ny * Nz. For access +by other commands, the name of the grid produced by this command is +"grid". The name of its data is "data". The (x,y,z) and (xc,yc,zc) coordinates are in distance :doc:`units `. diff --git a/doc/src/compute_temp.rst b/doc/src/compute_temp.rst index 3163de2272..c4984d97cc 100644 --- a/doc/src/compute_temp.rst +++ b/doc/src/compute_temp.rst @@ -29,7 +29,9 @@ Description Define a computation that calculates the temperature of a group of atoms. A compute of this style can be used by any command that -computes a temperature, e.g. :doc:`thermo_modify `, :doc:`fix temp/rescale `, :doc:`fix npt `, etc. +computes a temperature, e.g. :doc:`thermo_modify `, +:doc:`fix temp/rescale `, :doc:`fix npt `, +etc. The temperature is calculated by the formula KE = dim/2 N k T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m v\^2), @@ -79,8 +81,8 @@ Output info This compute calculates a global scalar (the temperature) and a global vector of length 6 (KE tensor), which can be accessed by indices 1-6. 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 ` page for an overview of LAMMPS output -options. +vector values from a compute as input. See the :doc:`Howto output +` page for an overview of LAMMPS output options. The scalar value calculated by this compute is "intensive". The vector values are "extensive". diff --git a/doc/src/dump.rst b/doc/src/dump.rst index a141b66e79..78f019adb9 100644 --- a/doc/src/dump.rst +++ b/doc/src/dump.rst @@ -162,7 +162,7 @@ Syntax possible attributes = c_ID:gname:dname, c_ID:gname:dname[I], f_ID:gname:dname, f_ID:gname:dname[I] gname = name of grid defined by compute or fix - ename = name of data field defined by compute or fix + dname = name of data field defined by compute or fix c_ID = per-grid vector calculated by a compute with ID c_ID[I] = Ith column of per-grid array calculated by a compute with ID, I can include wildcard (see below) f_ID = per-grid vector calculated by a fix with ID @@ -818,7 +818,7 @@ the distance and energy of each bond: This section explains the per-grid attributes that can be specified as part of the *grid* style. -The attributes that begin with *c_ID** and *f_ID* both take +The attributes that begin with *c_ID* and *f_ID* both take colon-separated fields *gname* and *dname*. These refer to a grid name and data field name which is defined by the compute or fix. Note that a compute or fix can define one or more grids (of different @@ -829,9 +829,7 @@ The *c_ID:gname:dname* and *c_ID:gname:dname[I]* attributes allow per-grid vectors or arrays calculated by a :doc:`compute ` to be output. The ID in the attribute should be replaced by the actual ID of the compute that has been defined previously in the input -script. See the :doc:`compute ` command for details. There -are computes for calculating local information such as indices, types, -and energies for bonds and angles. +script. If *c_ID:gname:dname* is used as a attribute, then the per-grid vector calculated by the compute is printed. If *c_ID:gname:dname[I]* is diff --git a/doc/src/fix_ave_chunk.rst b/doc/src/fix_ave_chunk.rst index cda706a217..5acdb89cd8 100644 --- a/doc/src/fix_ave_chunk.rst +++ b/doc/src/fix_ave_chunk.rst @@ -137,10 +137,11 @@ quantities. :doc:`Variables ` of style *atom* are the only ones that can be used with this fix since all other styles of variable produce global quantities. -Note that for values from a compute or fix, the bracketed index I can -be specified using a wildcard asterisk with the index to effectively -specify multiple values. This takes the form "\*" or "\*n" or "n\*" or -"m\*n". If N = the size of the vector (for *mode* = scalar) or the +Note that for values from a compute or fix that produces a per-atom +array (multiple values per atom), the bracketed index I can be +specified using a wildcard asterisk with the index to effectively +specify multiple values. This takes the form "\*" or "\*n" or "n\*" +or "m\*n". If N = the size of the vector (for *mode* = scalar) or the number of columns in the array (for *mode* = vector), then an asterisk with no numeric values means all indices from 1 to N. A leading asterisk means all indices from 1 to n (inclusive). A trailing diff --git a/doc/src/fix_ave_grid.rst b/doc/src/fix_ave_grid.rst new file mode 100644 index 0000000000..1c454a609e --- /dev/null +++ b/doc/src/fix_ave_grid.rst @@ -0,0 +1,436 @@ +.. index:: fix ave/grid + +fix ave/grid command +===================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID ave/grid Nevery Nrepeat Nfreq Nx Ny Nz value1 value2 ... keyword args ... + +* ID, group-ID are documented in :doc:`fix ` command +* ave/grid = style name of this fix command +* Nevery = use input values every this many timesteps +* Nrepeat = # of times to use input values for calculating averages +* Nfreq = calculate averages every this many timesteps +* Nx, Ny, Nz = grid size in each dimension +* one or more per-atom or per-grid input values can be listed +* per-atom value = vx, vy, vz, fx, fy, fz, density/mass, density/number, mass, temp, c_ID, c_ID[I], f_ID, f_ID[I], v_name + + .. parsed-literal:: + + vx,vy,vz,fx,fy,fz,mass = atom attribute (velocity, force component, mass) + density/number, density/mass = number or mass density (per volume) + temp = temperature + c_ID = per-atom vector calculated by a compute with ID + c_ID[I] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below) + f_ID = per-atom vector calculated by a fix with ID + f_ID[I] = Ith column of per-atom array calculated by a fix with ID, I can include wildcard (see below) + v_name = per-atom vector calculated by an atom-style variable with name + +* per-grid value = c_ID:gname:dname, c_ID:gname:dname[I], f_ID:gname:dname, f_ID:gname:dname[I] + + .. parsed-literal:: + + gname = name of grid defined by compute or fix + dname = name of data field defined by compute or fix + c_ID = per-grid vector calculated by a compute with ID + c_ID[I] = Ith column of per-grid array calculated by a compute with ID, I can include wildcard (see below) + f_ID = per-grid vector calculated by a fix with ID + f_ID[I] = Ith column of per-grid array calculated by a fix with ID, I can include wildcard (see below) + +* zero or more keyword/arg pairs may be appended +* keyword = *norm* or *ave* or *bias* or *adof* or *cdof* + + .. parsed-literal:: + + *norm* arg = *all* or *sample* or *none* = how output on *Nfreq* steps is normalized + all = output is sum of atoms across all *Nrepeat* samples, divided by atom count + sample = output is sum of *Nrepeat* sample averages, divided by *Nrepeat* + none = output is sum of *Nrepeat* sample sums, divided by *Nrepeat* + *ave* args = *one* or *running* or *window M* + one = output new average value every Nfreq steps + running = output cumulative average of all previous Nfreq steps + window M = output average of M most recent Nfreq steps + *bias* arg = bias-ID + bias-ID = ID of a temperature compute that removes a velocity bias for temperature calculation + *adof* value = dof_per_atom + dof_per_atom = define this many degrees-of-freedom per atom for temperature calculation + *cdof* value = dof_per_grid_cell + dof_per_grid_cell = add this many degrees-of-freedom per grid_cell for temperature calculation + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all ave/grid 10000 1 10000 10 10 10 fx fy fz c_myMSD[*] + fix 1 flow ave/chunk 100 10 1000 20 20 30 f_TTM:grid:data + +Description +""""""""""" + +Overlay the 2d or 3d simulation box with a uniformly spaced 2d or 3d +grid and use it to either (a) time-average per-atom quantities for the +atoms in each grid cell, or to (b) time-average per-grid quantities +produced by other computes or fixes. This fix operates in either +"per-atom mode" (all input values are per-atom) or in "per-grid mode" +(all input values are per-grid). + +The grid created by this command is distributed; each processor owns +the grid points that are within its subdomain. This is in contrast to +the :doc:`fix ave/chunk ` command when it uses chunks +from the :doc:`compute chunk/atom ` command which +are 2d or 3d regular bins. The per-bin outputs in that case are +global; each processor stores a copy of the entire set of bin data. +Thus it is better to use this command when the grid is large and a +simulation is run on many processors. + +For per-atom mode, only atoms in the specified group contribute to the +summing and averaging calculations. For per-grid mode, the specified +group is ignored. + +---------- + +The *Nevery*, *Nrepeat*, and *Nfreq* arguments specify on what +timesteps the input values will be accessed and contribute to the +average. The final averaged quantities are generated on timesteps +that are a multiples of *Nfreq*\ . The average is over *Nrepeat* +quantities, computed in the preceding portion of the simulation every +*Nevery* timesteps. *Nfreq* must be a multiple of *Nevery* and +*Nevery* must be non-zero even if *Nrepeat* is 1. Also, the timesteps +contributing to the average value cannot overlap, i.e. Nrepeat\*Nevery +can not exceed Nfreq. + +For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on +timesteps 90,92,94,96,98,100 will be used to compute the final average +on timestep 100. Similarly for timesteps 190,192,194,196,198,200 on +timestep 200, etc. If Nrepeat=1 and Nfreq = 100, then no time +averaging is done; values are simply generated on timesteps +100,200,etc. + +In per-atom mode, each input value can also be averaged over the atoms +in each grid cell. The way the averaging is done across the *Nrepeat* +timesteps to produce output on the *Nfreq* timesteps, and across +multiple *Nfreq* outputs, is determined by the *norm* and *ave* +keyword settings, as discussed below. + +---------- + +In both per-atom and per-grid mode, input values from a compute or fix +that produces an array of values (multiple values per atom or per grid +point), the bracketed index I can be specified using a wildcard +asterisk with the index to effectively specify multiple values. This +takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the number of +columns in the array (for *mode* = vector), then an asterisk with no +numeric values means all indices from 1 to N. A leading asterisk +means all indices from 1 to n (inclusive). A trailing asterisk means +all indices from n to N (inclusive). A middle asterisk means all +indices from m to n (inclusive). + +Using a wildcard is the same as if the individual columns of the array +had been listed one by one. E.g. if there were a compute fft/grid +command which produced 3 values for each grid point, these two fix +ave/grid commands would be equivalent: + +.. code-block:: LAMMPS + + compute myFFT all fft/grid 10 10 10 ... + fix 1 all ave/grid 100 1 100 10 10 10 c_myFFT:grid:data[*] + fix 2 all ave/grid 100 1 100 10 10 10 c_myFFT:grid:data[*][1] c_myFFT:grid:data[*][2] c_myFFT:grid:data[3] + +---------- + +*Per-atom mode*: + +Each specified per-atom value can be an atom attribute (velocity, +force component), a number or mass density, a mass or temperature, or +the result of a :doc:`compute ` or :doc:`fix ` or the +evaluation of an atom-style :doc:`variable `. In the latter +cases, the compute, fix, or variable must produce a per-atom quantity, +not a global quantity. Note that the :doc:`compute property/atom +` command provides access to any attribute +defined and stored by atoms. + +The per-atom values of each input vector are summed and averaged +independently of the per-atom values in other input vectors. + +:doc:`Computes ` that produce per-atom quantities are those +which have the word *atom* in their style name. See the doc pages for +individual :doc:`fixes ` to determine which ones produce per-atom +quantities. :doc:`Variables ` of style *atom* are the only +ones that can be used with this fix since all other styles of variable +produce global quantities. + +---------- + +The atom attribute values (vx,vy,vz,fx,fy,fz,mass) are +self-explanatory. As noted above, any other atom attributes can be +used as input values to this fix by using the :doc:`compute +property/atom ` command and then specifying an +input value from that compute. + +The *density/number* value means the number density is computed for +each grid cell, i.e. number/volume. The *density/mass* value means +the mass density is computed for each grid/cell, +i.e. total-mass/volume. The output values are in units of 1/volume or +density (mass/volume). See the :doc:`units ` command page for +the definition of density for each choice of units, e.g. gram/cm\^3. + +The *temp* value means the temperature is computed for each grid cell, +by the formula KE = DOF/2 k T, where KE = total kinetic energy of the +atoms in the grid cell (sum of 1/2 m v\^2), DOF = the total number of +degrees of freedom for all atoms in the grid cell, k = Boltzmann +constant, and T = temperature. + +The DOF is calculated as N\*adof + cdof, where N = number of atoms in +the grid cell, adof = degrees of freedom per atom, and cdof = degrees +of freedom per grid cell. By default adof = 2 or 3 = dimensionality +of system, as set via the :doc:`dimension ` command, and +cdof = 0.0. This gives the usual formula for temperature. + +Note that currently this temperature only includes translational +degrees of freedom for each atom. No rotational degrees of freedom +are included for finite-size particles. Also no degrees of freedom +are subtracted for any velocity bias or constraints that are applied, +such as :doc:`compute temp/partial `, or +:doc:`fix shake ` or :doc:`fix rigid `. This is +because those degrees of freedom (e.g. a constrained bond) could apply +to sets of atoms that are both inside and outside a specific grid +cell, and hence the concept is somewhat ill-defined. In some cases, +you can use the *adof* and *cdof* keywords to adjust the calculated +degrees of freedom appropriately, as explained below. + +Also note that a bias can be subtracted from atom velocities before +they are used in the above formula for KE, by using the *bias* +keyword. This allows, for example, a thermal temperature to be +computed after removal of a flow velocity profile. + +Note that the per-grid-cell temperature calculated by this fix and the +:doc:`compute temp/chunk ` command (using bins) +can be different. The compute calculates the temperature for each +chunk for a single snapshot. This fix can do that but can also time +average those values over many snapshots, or it can compute a +temperature as if the atoms in the grid cell on different timesteps +were collected together as one set of atoms to calculate their +temperature. The compute allows the center-of-mass velocity of each +chunk to be subtracted before calculating the temperature; this fix +does not. + +If a value begins with "c\_", a compute ID must follow which has been +previously defined in the input script. If no bracketed integer is +appended, the per-atom vector calculated by the compute is used. If a +bracketed integer is appended, the Ith column of the per-atom 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 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. If no bracketed integer is +appended, the per-atom vector calculated by the fix is used. If a +bracketed integer is appended, the Ith column of the per-atom array +calculated by the fix is used. Note that some fixes only produce +their values on certain timesteps, which must be compatible with +*Nevery*, else an error results. Users can also write code for their +own fix styles and :doc:`add them to LAMMPS `. See the +discussion above for how 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. Variables of style +*atom* 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 average within grid cells. + +---------- + +*Per-grid mode*: + +The attributes that begin with *c_ID* and *f_ID* both take +colon-separated fields *gname* and *dname*. These refer to a grid +name and data field name which is defined by the compute or fix. Note +that a compute or fix can define one or more grids (of different +sizes) and one or more data fields for each of those grids. The sizes +of all grids used as values for one instance of this fix must be the +same. + +The *c_ID:gname:dname* and *c_ID:gname:dname[I]* attributes allow +per-grid vectors or arrays calculated by a :doc:`compute ` to +be accessed. The ID in the attribute should be replaced by the actual +ID of the compute that has been defined previously in the input +script. + +If *c_ID:gname:dname* is used as a attribute, then the per-grid vector +calculated by the compute is accessed. If *c_ID:gname:dname[I]* is +used, then I must be in the range from 1-M, which will acccess the Ith +column of the per-grid array with M columns calculated by the compute. +See the discussion above for how I can be specified with a wildcard +asterisk to effectively specify multiple values. + +The *f_ID:gname:dname* and *f_ID:gname:dname[I]* attributes allow +per-grid vectors or arrays calculated by a :doc:`fix ` to be +output. The ID in the attribute should be replaced by the actual ID +of the fix that has been defined previously in the input script. + +If *f_ID:gname:dname* is used as a attribute, then the per-grid vector +calculated by the fix is printed. If *f_ID:gname:dname[I]* is used, +then I must be in the range from 1-M, which will print the Ith column +of the per-grid with M columns calculated by the fix. See the +discussion above for how I can be specified with a wildcard asterisk +to effectively specify multiple values. + +---------- + +Additional optional keywords also affect the operation of this fix and +its outputs. Some are only applicable to per-atom mode. Some are +applicable to both per-atom and per-grid mode. + +The *norm* keyword is only applicable to per-atom mode. It affects +how averaging is done for the per-grid values that are output once +every *Nfreq* timesteps when *Nrepeat* samples contribute to the +output. It has 3 possible settings: *all* or *sample* or *none*. +*All* is the default. + +In per-atom mode, *norm all* means the output grid value is summed +over all atoms in all *Nrepeat* samples, as is the count of atoms in +each grid cell. The averaged output value for a grid cell on the +*Nfreq* timesteps is Total-sum / Total-count. In other words it is an +average over atoms across the entire *Nfreq* timescale. For the +*density/number* and *density/mass* values, the grid cell volume used +in the final normalization will be the volume at the final *Nfreq* +timestep. For the *temp* values, degrees of freedom and kinetic energy +are summed separately across the entire *Nfreq* timescale, and the +output value is calculated by dividing those two sums. + +In per-atom mode, *norm sample* means the output grid value is summed +over atoms for each sample, as is the count, and an "average sample +value" is computed for each sample, i.e. Sample-sum / Sample-count. +The output grid value on the *Nfreq* timesteps is the average of the +*Nrepeat* "average sample values", i.e. the sum of *Nrepeat* "average +sample values" divided by *Nrepeat*\ . In other words it is an +average of an average. For the *density/number* and *density/mass* +values, the grid cell volume used in the per-sample normalization will +be the current grid cell volume at each sampling step. + +In per-atom mode, *norm none* perfomrma a similar computation as *norm +sample*, except the individual "average sample values" are "summed +sample values". A summed sample value is simply the grid value summed +over atoms in the sample, without dividing by the number of atoms in +the sample. The output grid value on the *Nfreq* timesteps is the +average of the *Nrepeat* "summed sample values", i.e. the sum of +*Nrepeat* "summed sample values" divided by *Nrepeat*\ . For the +*density/number* and *density/mass* values, the grid cell volume used +in the per-sample sum normalization will be the current grid cell +volume at each sampling step. + +In per-grid mode, all the *norm* keyword options operate the same. +The output grid value is summed over the grid value in each of the +*Nrepeat* samples and then divided by *Nrepeat*. + +The *ave* keyword is applicated to both per-atom and per-grid mode. +Itdetermines how the per-grid values produced once every *Nfreq* steps +are averaged with values produced on previous steps that were +multiples of *Nfreq*, before they are accessed by another output +command. + +If the *ave* setting is *one*, which is the default, then the grid +values produced on *Nfreq* timesteps are independent of each other; +they are output as-is without further averaging. + +If the *ave* setting is *running*, then the grid values produced on +*Nfreq* timesteps are summed and averaged in a cumulative sense before +being output. Each output grid value is thus the average of the grid +value produced on that timestep with all preceding values for the same +grid value. This running average begins when the fix is defined; it +can only be restarted by deleting the fix via the :doc:`unfix ` +command, or re-defining the fix by re-specifying it. + +If the *ave* setting is *window*, then the grid values produced on +*Nfreq* timesteps are summed and averaged within a moving "window" of +time, so that the last M values for the same grid are used to produce +the output. E.g. if M = 3 and Nfreq = 1000, then the grid value +output on step 10000 will be the average of the grid values on steps +8000,9000,10000. Outputs on early steps will average over less than M +values if they are not available. + +The *bias*, *adof*, and *cdof* keywords are only applicable to +per-atom mode. + +The *bias* keyword specifies the ID of a temperature compute that +removes a "bias" velocity from each atom, specified as *bias-ID*\ . +It is only used when the *temp* value is calculated, to compute the +thermal temperature of each grid cell after the translational kinetic +energy components have been altered in a prescribed way, e.g. to +remove a flow velocity profile. See the doc pages for individual +computes that calculate a temperature to see which ones implement a +bias. + +The *adof* and *cdof* keywords define the values used in the degree of +freedom (DOF) formula described above for temperature calculation for +each grid cell. They are only used when the *temp* value is +calculated. They can be used to calculate a more appropriate +temperature in some cases. Here are 3 examples: + +If grid cells contain some number of water molecules and :doc:`fix +shake ` is used to make each molecule rigid, then you could +calculate a temperature with 6 degrees of freedom (DOF) (3 +translational, 3 rotational) per molecule by setting *adof* to 2.0. + +If :doc:`compute temp/partial ` is used with the +*bias* keyword to only allow the x component of velocity to contribute +to the temperature, then *adof* = 1.0 would be appropriate. + +Using *cdof* = -2 or -3 (for 2d or 3d simulations) will subtract out 2 +or 3 degrees of freedom for each grid cell, similar to how the +:doc:`compute temp ` command subtracts out 3 DOF for the +entire system. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. + +This fix computes a global array of values which can be accessed by +various :doc:`output commands `. The values can only be +accessed on timesteps that are multiples of *Nfreq* since that is when +averaging is performed. The global array has # of rows = the number +of grids *grid* as calculated by the specified :doc:`compute +grid/atom ` command. The # of columns = +M+1+Nvalues, where M = 1 to 4, depending on whether the optional +columns for OrigID and CoordN are used, as explained above. Following +the optional columns, the next column contains the count of atoms in +the grid, and the remaining columns are the Nvalue quantities. When +the array is accessed with a row I that exceeds the current number of +grids, than a 0.0 is returned by the fix instead of an error, since +the number of grids can vary as a simulation runs depending on how +that value is computed by the compute grid/atom command. + +The array values calculated by this fix are treated as "intensive", +since they are typically already normalized by the count of atoms in +each grid. + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. + +Restrictions +"""""""""""" +none + +Related commands +"""""""""""""""" + +:doc:`fix ave/atom `, :doc:`fix ave/chunk ` + +Default +""""""" + +The option defaults are norm = all, ave = one, and bias = none. diff --git a/doc/src/fix_ttm.rst b/doc/src/fix_ttm.rst index d146321395..5ecb2f30f4 100644 --- a/doc/src/fix_ttm.rst +++ b/doc/src/fix_ttm.rst @@ -371,6 +371,13 @@ electronic subsystem energies reported at the end of the timestep. The vector values calculated are "extensive". +Th fix ttm/grid command also calculates a per-grid vector which store +the electron temperature for each grid cell in temperature :doc:`units +`. The length of the vector (distributed across all +processors) is Nx * Ny * Nz. For access by other commands, the name +of the grid produced by fix ttm/grid is "grid". The name of its data +is "data". + No parameter of the fixes can be used with the *start/stop* keywords of the :doc:`run ` command. The fixes are not invoked during :doc:`energy minimization `. diff --git a/src/compute_property_grid.cpp b/src/compute_property_grid.cpp index 1e83c0eb28..8a33eab9e2 100644 --- a/src/compute_property_grid.cpp +++ b/src/compute_property_grid.cpp @@ -24,7 +24,8 @@ using namespace LAMMPS_NS; -enum { ID, X, Y, Z, XS, YS, ZS, XC, YC, ZC, XSC, YSC, ZSC }; +enum { LOW, CTR }; +enum { UNSCALED, SCALED }; #define DELTA 10000 @@ -61,49 +62,49 @@ ComputePropertyGrid::ComputePropertyGrid(LAMMPS *lmp, int narg, char **arg) : pack_choice[jarg] = &ComputePropertyGrid::pack_id; } else if (strcmp(arg[iarg], "ix") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_ix; + pack_choice[jarg] = &ComputePropertyGrid::pack_indices<0>; } else if (strcmp(arg[iarg], "iy") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_iy; + pack_choice[jarg] = &ComputePropertyGrid::pack_indices<1>; } else if (strcmp(arg[iarg], "iz") == 0) { if (dimension == 2) error->all(FLERR,"Compute property/grid for 2d cannot use z coord"); - pack_choice[jarg] = &ComputePropertyGrid::pack_iz; + pack_choice[jarg] = &ComputePropertyGrid::pack_indices<2>; } else if (strcmp(arg[iarg], "x") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_x; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "y") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_y; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "z") == 0) { if (dimension == 2) error->all(FLERR,"Compute property/grid for 2d cannot use z coord"); - pack_choice[jarg] = &ComputePropertyGrid::pack_z; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "xs") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_xs; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "ys") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_ys; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "zs") == 0) { if (dimension == 2) error->all(FLERR,"Compute property/grid for 2d cannot use z coord"); - pack_choice[jarg] = &ComputePropertyGrid::pack_zs; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; - } else if (strcmp(arg[iarg], "xc") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_xc; + } else if (strcmp(arg[iarg], "xc") == 0) { + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "yc") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_yc; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "zc") == 0) { if (dimension == 2) error->all(FLERR,"Compute property/grid for 2d cannot use z coord"); - pack_choice[jarg] = &ComputePropertyGrid::pack_zc; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "xsc") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_xsc; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "ysc") == 0) { - pack_choice[jarg] = &ComputePropertyGrid::pack_ysc; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else if (strcmp(arg[iarg], "zsc") == 0) { if (dimension == 2) error->all(FLERR,"Compute property/grid for 2d cannot use z coord"); - pack_choice[jarg] = &ComputePropertyGrid::pack_zsc; + pack_choice[jarg] = &ComputePropertyGrid::pack_coords; } else error->all(FLERR, "Illegal compute property/grid command"); } @@ -124,6 +125,13 @@ ComputePropertyGrid::~ComputePropertyGrid() /* ---------------------------------------------------------------------- */ +void ComputePropertyGrid::init() +{ + triclinic = domain->triclinic; +} + +/* ---------------------------------------------------------------------- */ + void ComputePropertyGrid::compute_pergrid() { invoked_pergrid = update->ntimestep; @@ -283,7 +291,7 @@ double ComputePropertyGrid::memory_usage() } /* ---------------------------------------------------------------------- - one method for every keyword compute property/grid can output + compute grid point IDs ------------------------------------------------------------------------- */ void ComputePropertyGrid::pack_id(int n) @@ -313,411 +321,220 @@ void ComputePropertyGrid::pack_id(int n) } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + compute grid indices via templating +------------------------------------------------------------------------- */ -void ComputePropertyGrid::pack_ix(int n) +template void ComputePropertyGrid::pack_indices(int n) { if (dimension == 2) { if (nvalues == 0) { for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = ix + 1; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (IDIM == 0) vec2d[iy][ix] = ix + 1; + if (IDIM == 1) vec2d[iy][ix] = iy + 1; + } } else { for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = ix + 1; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (IDIM == 0) array2d[iy][ix][n] = ix + 1; + if (IDIM == 1) array2d[iy][ix][n] = iy + 1; + } } + } else if (dimension == 3) { if (nvalues == 0) { for (int iz = nzlo_in; iz <= nzhi_in; iz++) for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = ix + 1; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (IDIM == 0) vec3d[iz][iy][ix] = ix + 1; + if (IDIM == 1) vec3d[iz][iy][ix] = iy + 1; + if (IDIM == 2) vec3d[iz][iy][ix] = iz + 1; + } } else { for (int iz = nzlo_in; iz <= nzhi_in; iz++) for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = ix + 1; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (IDIM == 0) array3d[iz][iy][ix][n] = ix + 1; + if (IDIM == 1) array3d[iz][iy][ix][n] = iy + 1; + if (IDIM == 2) array3d[iz][iy][ix][n] = iz + 1; + } } } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + compute LOW/CTR, SCALED/UNSCALED, orthogonal/triclinic grid point coords + via templating +------------------------------------------------------------------------- */ -void ComputePropertyGrid::pack_iy(int n) +template +void ComputePropertyGrid::pack_coords(int n) { - if (dimension == 2) { - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = iy + 1; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = iy + 1; - } - } else if (dimension == 3) { - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = iy + 1; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = iy + 1; - } - } -} + double boxlo,delta; + double lamda[3],xone[3]; -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_iz(int n) -{ - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = iz + 1; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = iz + 1; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_x(int n) -{ - double boxlo,dx; + // 2d grid if (dimension == 2) { - grid2d->get_box(0,boxlo,dx); - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = boxlo + ix*dx; + + // for coords which are orthogonal OR scaled + + if (!triclinic || MODE == SCALED) { + + if (MODE == UNSCALED) grid2d->get_box(IDIM,boxlo,delta); + if (MODE == SCALED) { + boxlo = 0.0; + if (IDIM == 0) delta = 1.0/nxgrid; + if (IDIM == 1) delta = 1.0/nygrid; + } + + if (nvalues == 0) { + for (int iy = nylo_in; iy <= nyhi_in; iy++) + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (POS == LOW) { + if (IDIM == 0) vec2d[iy][ix] = boxlo + ix*delta; + if (IDIM == 1) vec2d[iy][ix] = boxlo + iy*delta; + } + if (POS == CTR) { + if (IDIM == 0) vec2d[iy][ix] = boxlo + (ix+0.5)*delta; + if (IDIM == 1) vec2d[iy][ix] = boxlo + (iy+0.5)*delta; + } + } + + } else { + for (int iy = nylo_in; iy <= nyhi_in; iy++) + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (POS == LOW) { + if (IDIM == 0) array2d[iy][ix][n] = boxlo + ix*delta; + if (IDIM == 1) array2d[iy][ix][n] = boxlo + iy*delta; + } + if (POS == CTR) { + if (IDIM == 0) array2d[iy][ix][n] = boxlo + (ix+0.5)*delta; + if (IDIM == 1) array2d[iy][ix][n] = boxlo + (iy+0.5)*delta; + } + } + } + + // only for coords which are triclinic AND unscaled + } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = boxlo + ix*dx; + + double dx = 1.0/nxgrid; + double dy = 1.0/nygrid; + lamda[2] = 0.0; + + if (nvalues == 0) { + for (int iy = nylo_in; iy <= nyhi_in; iy++) { + lamda[1] = iy*dy; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + lamda[0] = ix*dx; + domain->lamda2x(lamda,xone); + if (IDIM == 0) vec2d[iy][ix] = xone[0]; + if (IDIM == 1) vec2d[iy][ix] = xone[1]; + } + } + + } else { + for (int iy = nylo_in; iy <= nyhi_in; iy++) { + lamda[1] = iy*dy; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + lamda[0] = ix*dx; + domain->lamda2x(lamda,xone); + if (IDIM == 0) array2d[iy][ix][n] = xone[0]; + if (IDIM == 1) array2d[iy][ix][n] = xone[1]; + } + } + } } - } else if (dimension == 3) { - grid3d->get_box(0,boxlo,dx); - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + ix*dx; + + // 3d grid + + } else if (dimension == 3) { + + // for coords which are orthogonal OR scaled + + if (!triclinic || MODE == SCALED) { + + if (MODE == UNSCALED) grid3d->get_box(IDIM,boxlo,delta); + if (MODE == SCALED) { + boxlo = 0.0; + if (IDIM == 0) delta = 1.0/nxgrid; + if (IDIM == 1) delta = 1.0/nygrid; + if (IDIM == 2) delta = 1.0/nzgrid; + } + + if (nvalues == 0) { + for (int iz = nzlo_in; iz <= nzhi_in; iz++) + for (int iy = nylo_in; iy <= nyhi_in; iy++) + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (POS == LOW) { + if (IDIM == 0) vec3d[iz][iy][ix] = boxlo + ix*delta; + if (IDIM == 1) vec3d[iz][iy][ix] = boxlo + iy*delta; + if (IDIM == 2) vec3d[iz][iy][ix] = boxlo + iz*delta; + } + if (POS == CTR) { + if (IDIM == 0) vec3d[iz][iy][ix] = boxlo + (ix+0.5)*delta; + if (IDIM == 1) vec3d[iz][iy][ix] = boxlo + (iy+0.5)*delta; + if (IDIM == 2) vec3d[iz][iy][ix] = boxlo + (iz+0.5)*delta; + } + } + + } else { + for (int iz = nzlo_in; iz <= nzhi_in; iz++) + for (int iy = nylo_in; iy <= nyhi_in; iy++) + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + if (POS == LOW) { + if (IDIM == 0) array3d[iz][iy][ix][n] = boxlo + ix*delta; + if (IDIM == 1) array3d[iz][iy][ix][n] = boxlo + iy*delta; + if (IDIM == 2) array3d[iz][iy][ix][n] = boxlo + iz*delta; + } + if (POS == CTR) { + if (IDIM == 0) array3d[iz][iy][ix][n] = boxlo + (ix+0.5)*delta; + if (IDIM == 1) array3d[iz][iy][ix][n] = boxlo + (iy+0.5)*delta; + if (IDIM == 2) array3d[iz][iy][ix][n] = boxlo + (iz+0.5)*delta; + } + } + } + + // only for coords which are triclinic AND unscaled + } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + ix*dx; + + double dx = 1.0/nxgrid; + double dy = 1.0/nygrid; + double dz = 1.0/nzgrid; + + if (nvalues == 0) { + for (int iz = nzlo_in; iz <= nzhi_in; iz++) { + lamda[2] = iz*dz; + for (int iy = nylo_in; iy <= nyhi_in; iy++) { + lamda[1] = iy*dy; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + lamda[0] = ix*dx; + domain->lamda2x(lamda,xone); + if (IDIM == 0) vec3d[iz][iy][ix] = xone[0]; + if (IDIM == 1) vec3d[iz][iy][ix] = xone[1]; + if (IDIM == 2) vec3d[iz][iy][ix] = xone[2]; + } + } + } + + } else { + for (int iz = nzlo_in; iz <= nzhi_in; iz++) { + lamda[2] = iz*dz; + for (int iy = nylo_in; iy <= nyhi_in; iy++) { + lamda[1] = iy*dy; + for (int ix = nxlo_in; ix <= nxhi_in; ix++) { + lamda[0] = ix*dx; + domain->lamda2x(lamda,xone); + if (IDIM == 0) array3d[iz][iy][ix][n] = xone[0]; + if (IDIM == 1) array3d[iz][iy][ix][n] = xone[1]; + if (IDIM == 2) array3d[iz][iy][ix][n] = xone[2]; + } + } + } + } } } } - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_y(int n) -{ - double boxlo,dy; - - if (dimension == 2) { - grid2d->get_box(1,boxlo,dy); - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = boxlo + iy*dy; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = boxlo + iy*dy; - } - } else if (dimension == 3) { - grid3d->get_box(1,boxlo,dy); - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + iy*dy; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + iy*dy; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_z(int n) -{ - double boxlo,dz; - grid3d->get_box(2,boxlo,dz); - - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + iz*dz; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + iz*dz; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_xs(int n) -{ - double dx = 1.0/nxgrid; - - if (dimension == 2) { - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = ix*dx; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = ix*dx; - } - } else if (dimension == 3) { - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = ix*dx; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = ix*dx; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_ys(int n) -{ - double dy = 1.0/nygrid; - - if (dimension == 2) { - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = iy*dy; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = iy*dy; - } - } else if (dimension == 3) { - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = iy*dy; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = iy*dy; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_zs(int n) -{ - double dz = 1.0/nzgrid; - - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = iz*dz; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = iz*dz; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_xc(int n) -{ - double boxlo,dx; - - if (dimension == 2) { - grid2d->get_box(0,boxlo,dx); - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = boxlo + (ix+0.5)*dx; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = boxlo + (ix+0.5)*dx; - } - } else if (dimension == 3) { - grid3d->get_box(0,boxlo,dx); - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + (ix+0.5)*dx; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + (ix+0.5)*dx; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_yc(int n) -{ - double boxlo,dy; - - if (dimension == 2) { - grid2d->get_box(1,boxlo,dy); - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = boxlo + (iy+0.5)*dy; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = boxlo + (iy+0.5)*dy; - } - } else if (dimension == 3) { - grid3d->get_box(1,boxlo,dy); - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + (iy+0.5)*dy; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + (iy+0.5)*dy; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_zc(int n) -{ - double boxlo,dz; - grid3d->get_box(2,boxlo,dz); - - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = boxlo + (iz+0.5)*dz; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = boxlo + (iz+0.5)*dz; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_xsc(int n) -{ - double dx = 1.0/nxgrid; - - if (dimension == 2) { - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = (ix+0.5)*dx; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = (ix+0.5)*dx; - } - } else if (dimension == 3) { - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = (ix+0.5)*dx; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = (ix+0.5)*dx; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_ysc(int n) -{ - double dy = 1.0/nygrid; - - if (dimension == 2) { - if (nvalues == 0) { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec2d[iy][ix] = (iy+0.5)*dy; - } else { - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array2d[iy][ix][n] = (iy+0.5)*dy; - } - } else if (dimension == 3) { - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = (iy+0.5)*dy; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = (iy+0.5)*dy; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyGrid::pack_zsc(int n) -{ - double dz = 1.0/nzgrid; - - if (nvalues == 0) { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - vec3d[iz][iy][ix] = (iz+0.5)*dz; - } else { - for (int iz = nzlo_in; iz <= nzhi_in; iz++) - for (int iy = nylo_in; iy <= nyhi_in; iy++) - for (int ix = nxlo_in; ix <= nxhi_in; ix++) - array3d[iz][iy][ix][n] = (iz+0.5)*dz; - } -} diff --git a/src/compute_property_grid.h b/src/compute_property_grid.h index 4948888ca7..929346ad7a 100644 --- a/src/compute_property_grid.h +++ b/src/compute_property_grid.h @@ -28,7 +28,7 @@ class ComputePropertyGrid : public Compute { public: ComputePropertyGrid(class LAMMPS *, int, char **); ~ComputePropertyGrid() override; - void init() override {} + void init() override; void compute_pergrid() override; void reset_grid() override; @@ -44,6 +44,7 @@ class ComputePropertyGrid : public Compute { int nxgrid,nygrid,nzgrid; int nvalues; int dimension; + int triclinic; class Grid2d *grid2d; class Grid3d *grid3d; @@ -64,26 +65,8 @@ class ComputePropertyGrid : public Compute { FnPtrPack *pack_choice; // ptrs to pack functions void pack_id(int); - - void pack_ix(int); - void pack_iy(int); - void pack_iz(int); - - void pack_x(int); - void pack_y(int); - void pack_z(int); - - void pack_xs(int); - void pack_ys(int); - void pack_zs(int); - - void pack_xc(int); - void pack_yc(int); - void pack_zc(int); - - void pack_xsc(int); - void pack_ysc(int); - void pack_zsc(int); + template void pack_indices(int); + template void pack_coords(int); }; } // namespace LAMMPS_NS diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index f6d2a3707e..2b63d04e3c 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -64,7 +64,7 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : pergrid_freq = utils::inumeric(FLERR,arg[5],false,lmp); time_depend = 1; - // NOTE: allow Dxyz as well + // NOTE: allow Dxyz as well at some point nxgrid = utils::inumeric(FLERR,arg[6],false,lmp); nygrid = utils::inumeric(FLERR,arg[7],false,lmp); @@ -957,14 +957,12 @@ void FixAveGrid::atom2grid() n = value2index[m]; j = argindex[m]; - // X,V,F adds coord,velocity,force to value + // V,F adds velocity,force to value - if (which[m] == ArgInfo::X || which[m] == ArgInfo::V || - which[m] == ArgInfo::F) { + if (which[m] == ArgInfo::V || which[m] == ArgInfo::F) { double **attribute; - if (which[m] == ArgInfo::X) attribute = atom->x; - else if (which[m] == ArgInfo::V) attribute = atom->v; + if (which[m] == ArgInfo::V) attribute = atom->v; else if (which[m] == ArgInfo::F) attribute = atom->f; if (dimension == 2) {