ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

finished first draft of doc pages

Browse Source
This commit is contained in:
Steve Plimpton
2022-08-16 15:42:48 -06:00
parent 206ab141c9
commit 5e935519bf
10 changed files with 702 additions and 449 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
doc/src/compute_modify.rst
View File

@ -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 <dimension>` 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 <compute_temp_asphere>` command for an example.
temperature use this option. The default is 2 or 3 for :doc:`2d or 3d
systems <dimension>` 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 <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 <fix_pour>`, :doc:`fix deposit <fix_deposit>`, and :doc:`fix gcmc <fix_gcmc>` commands) or expect atoms or molecules to be lost
(e.g. due to exiting the simulation box or via :doc:`fix evaporate <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
<fix_pour>`, :doc:`fix deposit <fix_deposit>`, and :doc:`fix gcmc
<fix_gcmc>` commands) or expect atoms or molecules to be lost
(e.g. due to exiting the simulation box or via :doc:`fix evaporate
<fix_evaporate>`), then this option should be used to insure the
temperature is correctly normalized.
.. note::

16
doc/src/compute_property_grid.rst
View File

@ -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 <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 <Howto_output>`
such as :doc:`dump grid <dump>`.
*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 <dump>`,
*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
<units>`.

8
doc/src/compute_temp.rst
View File

@ -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 <thermo_modify>`, :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix npt <fix_nh>`, etc.
computes a temperature, e.g. :doc:`thermo_modify <thermo_modify>`,
:doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix npt <fix_nh>`,
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 <Howto_output>` page for an overview of LAMMPS output
options.
vector values from a compute as input. See the :doc:`Howto output
<Howto_output>` page for an overview of LAMMPS output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".

8
doc/src/dump.rst
View File

@ -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 <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 <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

9
doc/src/fix_ave_chunk.rst
View File

@ -137,10 +137,11 @@ quantities. :doc:`Variables <variable>` 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

436
doc/src/fix_ave_grid.rst Normal file
View File

@ -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 <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 <fix_ave_chunk>` command when it uses chunks
from the :doc:`compute chunk/atom <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 <compute>` or :doc:`fix <fix>` or the
evaluation of an atom-style :doc:`variable <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
<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 <compute>` that produce per-atom quantities are those
which have the word *atom* in their style name. See the doc pages for
individual :doc:`fixes <fix>` to determine which ones produce per-atom
quantities. :doc:`Variables <variable>` 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 <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 <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 <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 <compute_temp_partial>`, or
:doc:`fix shake <fix_shake>` or :doc:`fix rigid <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 <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 <Modify>`. 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 <Modify>`. 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 <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 <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 <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 <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 <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 <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
<restart>`. None of the :doc:`fix_modify <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 <Howto_output>`. 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 <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 <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`fix ave/atom <fix_ave_atom>`, :doc:`fix ave/chunk <fix_ave_chunk>`
Default
"""""""
The option defaults are norm = all, ave = one, and bias = none.

7
doc/src/fix_ttm.rst
View File

@ -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
<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 <run>` command. The fixes are not invoked during
:doc:`energy minimization <minimize>`.

555
src/compute_property_grid.cpp
View File

@ -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<LOW,UNSCALED,0>;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_y;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<LOW,UNSCALED,1>;
} 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<LOW,UNSCALED,2>;
} else if (strcmp(arg[iarg], "xs") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_xs;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<LOW,SCALED,0>;
} else if (strcmp(arg[iarg], "ys") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_ys;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<LOW,SCALED,1>;
} 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<LOW,SCALED,2>;
} else if (strcmp(arg[iarg], "xc") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_xc;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<CTR,UNSCALED,0>;
} else if (strcmp(arg[iarg], "yc") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_yc;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<CTR,UNSCALED,1>;
} 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<CTR,UNSCALED,2>;
} else if (strcmp(arg[iarg], "xsc") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_xsc;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<CTR,SCALED,0>;
} else if (strcmp(arg[iarg], "ysc") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_ysc;
pack_choice[jarg] = &ComputePropertyGrid::pack_coords<CTR,SCALED,1>;
} 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<CTR,SCALED,2>;
} 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 <int IDIM> 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 <int POS, int MODE, int IDIM>
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);
// 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++)
vec2d[iy][ix] = boxlo + ix*dx;
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++)
array2d[iy][ix][n] = boxlo + ix*dx;
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;
}
} 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;
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 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;
}
}
}
/* ---------------------------------------------------------------------- */
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;
lamda[2] = 0.0;
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;
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];
}
}
}
}
// 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++)
vec3d[iz][iy][ix] = iy*dy;
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++)
array3d[iz][iy][ix][n] = iy*dy;
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
void ComputePropertyGrid::pack_zs(int n)
{
} else {
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++)
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;
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];
}
}
/* ---------------------------------------------------------------------- */
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;
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_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;
}
}

25
src/compute_property_grid.h
View File

@ -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 <int IDIM> void pack_indices(int);
template <int POS, int MODE, int IDIM> void pack_coords(int);
};
} // namespace LAMMPS_NS

10
src/fix_ave_grid.cpp
View File

@ -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) {