enable per-atom custom arrays in addition to vectors

This commit is contained in:
Steve Plimpton
2020-09-04 11:33:49 -06:00
parent d00807ee9a
commit d1442b0538
23 changed files with 1390 additions and 722 deletions

View File

@ -465,12 +465,18 @@ atom-style variables can reference the position of a particle, its
velocity, the volume of its Voronoi cell, etc. velocity, the volume of its Voronoi cell, etc.
The *store* weight style does not compute a weight factor. Instead it The *store* weight style does not compute a weight factor. Instead it
stores the current accumulated weights in a custom per-atom property stores the current accumulated weights in a custom per-atom vector
specified by *name*\ . This must be a property defined as *d_name* via specified by *name*\ . This must be a vector defined as *d_name* via
the :doc:`fix property/atom <fix_property_atom>` command. Note that the :doc:`fix property/atom <fix_property_atom>` command. This means
these custom per-atom properties can be output in a :doc:`dump <dump>` the values in the vector can be read as part of a data file with the
file, so this is a way to examine, debug, or visualize the :doc:`read_data <read_data>` command or specified with the :doc:`set
per-particle weights computed during the load-balancing operation. <set>` command. These weights can also be output in a :doc:`dump
<dump>` file, so this is a way to examine, debug, or visualize the
per-particle weights used during the load-balancing operation.
Note that the name of the custom per-atom vector is specified just
as *name*, not as *d_name* as it is for other commands that use
different kinds of custom atom vectors or arrays as arguments.
---------- ----------

View File

@ -20,7 +20,8 @@ Syntax
x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz, x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
vx, vy, vz, fx, fy, fz, vx, vy, vz, fx, fy, fz,
q, mux, muy, muz, mu, q, mux, muy, muz, mu,
sp, spx, spy, spz, fmx, fmy, fmz, spx, spy, spz, sp, fmx, fmy, fmz,
nbonds,
radius, diameter, omegax, omegay, omegaz, radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, angmomx, angmomy, angmomz,
shapex,shapey, shapez, shapex,shapey, shapez,
@ -29,42 +30,42 @@ Syntax
corner1x, corner1y, corner1z, corner1x, corner1y, corner1z,
corner2x, corner2y, corner2z, corner2x, corner2y, corner2z,
corner3x, corner3y, corner3z, corner3x, corner3y, corner3z,
nbonds, i_name, d_name, i2_name[I], d2_name[I],
buckling, vfrac, s0, spin, eradius, ervel, erforce,
vfrac, s0, rho, drho, e, de, cv, buckling
spin, eradius, ervel, erforce,
rho, drho, e, de, cv,
i_name, d_name
.. parsed-literal:: .. parsed-literal::
id = atom ID *id* = atom ID
mol = molecule ID *mol* = molecule ID
proc = ID of processor that owns atom *proc* = ID of processor that owns atom
type = atom type *type* = atom type
mass = atom mass *mass* = atom mass
x,y,z = unscaled atom coordinates *x,y,z* = unscaled atom coordinates
xs,ys,zs = scaled atom coordinates *xs,ys,zs* = scaled atom coordinates
xu,yu,zu = unwrapped atom coordinates *xu,yu,zu* = unwrapped atom coordinates
ix,iy,iz = box image that the atom is in *ix,iy,iz* = box image that the atom is in
vx,vy,vz = atom velocities *vx,vy,vz* = atom velocities
fx,fy,fz = forces on atoms *fx,fy,fz* = forces on atoms
q = atom charge *q* = atom charge
mux,muy,muz = orientation of dipole moment of atom *mux,muy,muz* = orientation of dipole moment of atom
mu = magnitude of dipole moment of atom *mu* = magnitude of dipole moment of atom
sp = atomic magnetic spin moment *spx, spy, spz* = direction of the atomic magnetic spin
spx, spy, spz = direction of the atomic magnetic spin *sp* = magintude of atomic magnetic spin moment
fmx, fmy, fmz = magnetic force *fmx, fmy, fmz* = magnetic force
radius,diameter = radius,diameter of spherical particle *nbonds* = number of bonds assigned to an atom
omegax,omegay,omegaz = angular velocity of spherical particle *radius,diameter* = radius,diameter of spherical particle
angmomx,angmomy,angmomz = angular momentum of aspherical particle *omegax,omegay,omegaz* = angular velocity of spherical particle
shapex,shapey,shapez = 3 diameters of aspherical particle *angmomx,angmomy,angmomz* = angular momentum of aspherical particle
quatw,quati,quatj,quatk = quaternion components for aspherical or body particles *shapex,shapey,shapez* = 3 diameters of aspherical particle
tqx,tqy,tqz = torque on finite-size particles *quatw,quati,quatj,quatk* = quaternion components for aspherical or body particles
end12x, end12y, end12z = end points of line segment *tqx,tqy,tqz* = torque on finite-size particles
corner123x, corner123y, corner123z = corner points of triangle *end12x, end12y, end12z* = end points of line segment
nbonds = number of bonds assigned to an atom *corner123x, corner123y, corner123z* = corner points of triangle
buckling = buckling flag used in mesoscopic simulation of nanotubes *i_name* = custom integer vector with name
*d_name* = custom floating point vector with name
*i2_name[I]* = Ith column of custom integer array with name
*d2_name[I]* = Ith column of custom floating-point array with name
.. parsed-literal:: .. parsed-literal::
@ -91,9 +92,8 @@ Syntax
.. parsed-literal:: .. parsed-literal::
:doc:`fix property/atom <fix_property_atom>` per-atom properties: USER-MESONT package per-atom properties:
i_name = custom integer vector with name buckling = buckling flag used in mesoscopic simulation of nanotubes
d_name = custom integer vector with name
Examples Examples
"""""""" """"""""
@ -104,6 +104,7 @@ Examples
compute 2 all property/atom type compute 2 all property/atom type
compute 1 all property/atom ix iy iz compute 1 all property/atom ix iy iz
compute 3 all property/atom sp spx spy spz compute 3 all property/atom sp spx spy spz
compute 1 all property/atom i_myFlag d_Sxyz[1] d_Sxyz[3]
Description Description
""""""""""" """""""""""
@ -116,20 +117,37 @@ ave/atom <fix_ave_atom>`, :doc:`fix ave/histo <fix_ave_histo>`,
:doc:`fix ave/chunk <fix_ave_chunk>`, and :doc:`atom-style variable :doc:`fix ave/chunk <fix_ave_chunk>`, and :doc:`atom-style variable
<variable>` commands. <variable>` commands.
The list of possible attributes is the same as that used by the The list of possible attributes is essentially the same as that used
:doc:`dump custom <dump>` command, which describes their meaning, with by the :doc:`dump custom <dump>` command, which describes their
some additional quantities that are only defined for certain meaning, with some additional quantities that are only defined for
:doc:`atom styles <atom_style>`. Basically, this augmented list gives certain :doc:`atom styles <atom_style>`. The goal of this augmented
an input script access to any per-atom quantity stored by LAMMPS. list gives an input script access to any per-atom quantity stored by
LAMMPS.
The values are stored in a per-atom vector or array as discussed The values are stored in a per-atom vector or array as discussed
below. Zeroes are stored for atoms not in the specified group or for below. Zeroes are stored for atoms not in the specified group or for
quantities that are not defined for a particular particle in the group quantities that are not defined for a particular particle in the group
(e.g. *shapex* if the particle is not an ellipsoid). (e.g. *shapex* if the particle is not an ellipsoid).
Attributes *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom
per-atom integer and floating-point vectors or arrays that have been
added via the :doc:`fix property/atom <fix_property_atom>` command.
When that command is used specific names are given to each attribute
which are the "name" portion of these attributes. For arrays *i2_name*
and *d2_name*, the column of the array must also be included following
the name in brackets: e.g. d2_xyz[2], i2_mySpin[3].
The additional quantities only accessible via this command, and not The additional quantities only accessible via this command, and not
directly via the :doc:`dump custom <dump>` command, are as follows. directly via the :doc:`dump custom <dump>` command, are as follows.
*Nbonds* is available for all molecular atom styles and refers to the
number of explicit bonds assigned to an atom. Note that if the
:doc:`newton bond <newton>` command is set to *on*\ , which is the
default, then every bond in the system is assigned to only one of the
two atoms in the bond. Thus a bond between atoms I,J may be tallied
for either atom I or atom J. If :doc:`newton bond off <newton>` is
set, it will be tallied with both atom I and atom J.
*Shapex*\ , *shapey*\ , and *shapez* are defined for ellipsoidal particles *Shapex*\ , *shapey*\ , and *shapez* are defined for ellipsoidal particles
and define the 3d shape of each particle. and define the 3d shape of each particle.
@ -146,19 +164,8 @@ line segment.
*corner2z*\ , *corner3x*\ , *corner3y*\ , *corner3z*\ , are defined for *corner2z*\ , *corner3x*\ , *corner3y*\ , *corner3z*\ , are defined for
triangular particles and define the corner points of each triangle. triangular particles and define the corner points of each triangle.
*Nbonds* is available for all molecular atom styles and refers to the In addition, the various per-atom quantities listed above for specific
number of explicit bonds assigned to an atom. Note that if the packages are only accessible by this command.
:doc:`newton bond <newton>` command is set to *on*\ , which is the
default, then every bond in the system is assigned to only one of the
two atoms in the bond. Thus a bond between atoms I,J may be tallied
for either atom I or atom J. If :doc:`newton bond off <newton>` is
set, it will be tallied with both atom I and atom J.
The *i_name* and *d_name* attributes refer to custom integer and
floating-point properties that have been added to each atom via the
:doc:`fix property/atom <fix_property_atom>` command. When that
command is used specific names are given to each attribute which are
what is specified as the "name" portion of *i_name* or *d_name*.
Output info Output info
""""""""""" """""""""""

View File

@ -78,7 +78,8 @@ Syntax
q, mux, muy, muz, mu, q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz, radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz, angmomx, angmomy, angmomz, tqx, tqy, tqz,
c_ID, c_ID[N], f_ID, f_ID[N], v_name c_ID, c_ID[I], f_ID, f_ID[I], v_name,
i_name, d_name, i2_name[I], d2_name[I]
.. parsed-literal:: .. parsed-literal::
@ -108,8 +109,10 @@ Syntax
f_ID = per-atom vector calculated by a fix with ID 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) 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 v_name = per-atom vector calculated by an atom-style variable with name
d_name = per-atom floating point vector with name, managed by fix property/atom i_name = custom integer vector with name
i_name = per-atom integer vector with name, managed by fix property/atom d_name = custom floating point vector with name
i2_name[I] = Ith column of custom integer array with name, I can include wildcard (see below)
d2_name[I] = Ith column of custom floating point vector with name, I can include wildcard (see below)
* *local* args = list of local attributes * *local* args = list of local attributes
@ -134,7 +137,7 @@ Examples
dump 2 subgroup atom 50 dump.run.mpiio.bin dump 2 subgroup atom 50 dump.run.mpiio.bin
dump 4a all custom 100 dump.myforce.* id type x y vx fx dump 4a all custom 100 dump.myforce.* id type x y vx fx
dump 4b flow custom 100 dump.%.myforce id type c_myF[3] v_ke dump 4b flow custom 100 dump.%.myforce id type c_myF[3] v_ke
dump 4b flow custom 100 dump.%.myforce id type c_myF[\*] v_ke dump 4b flow custom 100 dump.%.myforce id type c_myF[*] v_ke
dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz
dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress[2] dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress[2]
dump 1 all xtc 1000 file.xtc dump 1 all xtc 1000 file.xtc
@ -465,16 +468,15 @@ styles.
---------- ----------
Note that in the discussion which follows, for styles which can Note that in the discussion which follows, for styles which can
reference values from a compute or fix, like the *custom*\ , *cfg*\ , or reference values from a compute or fix or custom atom property, like
*local* styles, the bracketed index I can be specified using a the *custom*\ , *cfg*\ , or *local* styles, the bracketed index I can
wildcard asterisk with the index to effectively specify multiple be specified using a wildcard asterisk with the index to effectively
values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the specify multiple values. This takes the form "\*" or "\*n" or "n\*"
size of the vector (for *mode* = scalar) or the number of columns in or "m\*n". If N = the number of columns in the array, then an
the array (for *mode* = vector), then an asterisk with no numeric asterisk with no numeric values means all column indices from 1 to N.
values means all indices from 1 to N. A leading asterisk means all A leading asterisk means all indices from 1 to n (inclusive). A
indices from 1 to n (inclusive). A trailing asterisk means all trailing asterisk means all indices from n to N (inclusive). A middle
indices from n to N (inclusive). A middle asterisk means all indices asterisk means all indices from m to n (inclusive).
from m to n (inclusive).
Using a wildcard is the same as if the individual columns of the array Using a wildcard is the same as if the individual columns of the array
had been listed one by one. E.g. these 2 dump commands are had been listed one by one. E.g. these 2 dump commands are
@ -512,8 +514,9 @@ bonds and angles.
Note that computes which calculate global or per-atom quantities, as Note that computes which calculate global or per-atom quantities, as
opposed to local quantities, cannot be output in a dump local command. opposed to local quantities, cannot be output in a dump local command.
Instead, global quantities can be output by the :doc:`thermo_style custom <thermo_style>` command, and per-atom quantities can be Instead, global quantities can be output by the :doc:`thermo_style
output by the dump custom command. custom <thermo_style>` command, and per-atom quantities can be output
by the dump custom command.
If *c_ID* is used as a attribute, then the local vector calculated by If *c_ID* is used as a attribute, then the local vector calculated by
the compute is printed. If *c_ID[I]* is used, then I must be in the the compute is printed. If *c_ID[I]* is used, then I must be in the
@ -660,9 +663,13 @@ invoke other computes, fixes, or variables when they are evaluated, so
this is a very general means of creating quantities to output to a this is a very general means of creating quantities to output to a
dump file. dump file.
The *d_name* and *i_name* attributes allow to output custom per atom The *i_name*, *d_name*, *i2_name*, *d2_name* attributes refer to
floating point or integer properties that are managed by per-atom integer and floating-point vectors or arrays that have been
:doc:`fix property/atom <fix_property_atom>`. added via the :doc:`fix property/atom <fix_property_atom>` command.
When that command is used specific names are given to each attribute
which are the "name" portion of these keywords. For arrays *i2_name*
and *d2_name*, the column of the array must also be included following
the name in brackets: e.g. d2_xyz[2], i2_mySpin[3].
See the :doc:`Modify <Modify>` doc page for information on how to add See the :doc:`Modify <Modify>` doc page for information on how to add
new compute and fix styles to LAMMPS to calculate per-atom quantities new compute and fix styles to LAMMPS to calculate per-atom quantities

View File

@ -11,11 +11,11 @@ Syntax
.. parsed-literal:: .. parsed-literal::
fix ID group-ID property/atom vec1 vec2 ... keyword value ... fix ID group-ID property/atom name1 name2 ... keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command * ID, group-ID are documented in :doc:`fix <fix>` command
* property/atom = style name of this fix command * property/atom = style name of this fix command
* vec1,vec2,... = *mol* or *q* or *rmass* or *i_name* or *d_name* * name1,name2,... = *mol* or *q* or *rmass* or *i_name* or *d_name* or *i2_name* or *d2_name*
.. parsed-literal:: .. parsed-literal::
@ -24,6 +24,10 @@ Syntax
*rmass* = per-atom mass *rmass* = per-atom mass
*i_name* = new integer vector referenced by name *i_name* = new integer vector referenced by name
*d_name* = new floating-point vector referenced by name *d_name* = new floating-point vector referenced by name
*i2_name* = new integer array referenced by name
i2_name arg = N = number of columns in the array
*d2_name* = new floating-point array referenced by name
d2_name arg = N = number of columns in the array
* zero of more keyword/value pairs may be appended * zero of more keyword/value pairs may be appended
* keyword = *ghost* * keyword = *ghost*
@ -39,58 +43,64 @@ Examples
fix 1 all property/atom mol fix 1 all property/atom mol
fix 1 all property/atom i_myflag1 i_myflag2 fix 1 all property/atom i_myflag1 i_myflag2
fix 1 all property/atom d_sx d_sy d_sz fix 1 all property/atom d2_sxyz 3 ghost yes
Description Description
""""""""""" """""""""""
Create one or more additional per-atom vectors to store information Create one or more additional per-atom vectors or arrays to store
about atoms and to use during a simulation. The specified *group-ID* information about atoms and to use during a simulation. The specified
is ignored by this fix. *group-ID* is ignored by this fix.
The atom style used for a simulation defines a set of per-atom The atom style used for a simulation defines a set of per-atom
properties, as explained on the :doc:`atom_style <atom_style>` and properties, as explained on the :doc:`atom_style <atom_style>` and
:doc:`read_data <read_data>` doc pages. The latter command allows these :doc:`read_data <read_data>` doc pages. The latter command defines
properties to be defined for each atom in the system when a data file these properties for each atom in the system when a data file is read.
is read. This fix will augment the set of properties with new custom This fix augments the set of per-atom properties with new custom
ones. This can be useful in several scenarios. ones. This can be useful in several scenarios.
If the atom style does not define molecule IDs, per-atom charge, If the atom style does not define molecule IDs, per-atom charge, or
or per-atom mass, they can be added using the *mol*\ , *q* or *rmass* per-atom mass, they can be added using the *mol*\ , *q* or *rmass*
keywords. This can be useful, e.g, to define "molecules" to use as keywords. This ciykd be useful to define "molecules" to use as rigid
rigid bodies with the :doc:`fix rigid <fix_rigid>` command, or just to bodies with the :doc:`fix rigid <fix_rigid>` command, or to carry
carry around an extra flag with the atoms (stored as a molecule ID) around an extra flag with atoms (stored as a molecule ID) that can be
that can be used to group atoms without having to use the group used by various commands like :doc:`compute chunk/atom
<compute_chunk_atom>` to group atoms without having to use the group
command (which is limited to a total of 32 groups including *all*\ ). command (which is limited to a total of 32 groups including *all*\ ).
Another application would be to use the *rmass* flag in order to have Another application is to use the *rmass* flag in order to have
per-atom masses instead of per-type masses, for example this can be per-atom masses instead of per-type masses. This could be used to
useful to study isotope effects with partial isotope substitution. study isotope effects with partial isotope substitution. :ref:`See
Please :ref:`see below <isotopes>` for an example of simulating a mixture below <isotopes>` for an example of simulating a mixture of light and
of light and heavy water with the TIP4P water potential. heavy water with the TIP4P water potential.
An alternative to using fix *property/atom* in these ways is to An alternative to using fix *property/atom* for these examples is to
use an atom style that does define molecule IDs or charge or per-atom use an atom style that does define molecule IDs or charge or per-atom
mass (indirectly via diameter and density) or to use a hybrid atom mass (indirectly via diameter and density) or to use a hybrid atom
style that combines two or more atom styles style that combines two or more atom styles to provide the union of
to provide the union of all atom properties. However, this has two all their atom properties. However, this has two practical drawbacks:
practical drawbacks: first, it typically necessitates changing the first, it typically necessitates changing the format of the Atoms
format of the data file, which can be tedious for large systems; section in the data file and second, it may define additional
and second, it may define additional properties that are not needed properties that are not needed such as bond lists, which incurs some
such as bond lists, which has some overhead when there are no bonds. overhead when there are no bonds.
In the future, we may add additional per-atom properties similar to In the future, we may add additional existing per-atom properties to
*mol*\ , *q* or *rmass*\ , which "turn-on" specific properties defined fix property/atom, similar to *mol*\ , *q* or *rmass*\ , which
by some atom styles, so they can be used by atom styles that do not "turn-on" specific properties defined by some atom styles, so they can
define them. be easily used by atom styles that do not define them.
More generally, the *i_name* and *d_name* vectors allow one or more More generally, the *i_name* and *d_name* options allow one or more
new custom per-atom properties to be defined. Each name must be new custom per-atom vectors to be defined. Likewise the *i2_name* and
unique and can use alphanumeric or underscore characters. These *d2_name* options allow one or more custom per-atom arrays to be
vectors can store whatever values you decide are useful in your defined. The *i2_name* and *d2_name* options take an argument *N*
simulation. As explained below there are several ways to initialize which specifies the number of columns in the per-atom array, i.e. the
and access and output these values, both via input script commands and number of attributes associated with each atom. *N* >= 1 is required.
in new code that you add to LAMMPS.
Each name must be unique and can use alphanumeric or underscore
characters. These vectors and arrays can store whatever values you
decide are useful in your simulation. As explained below there are
several ways to initialize, access, and output these values, via input
script commands, data files, and in new code you add to LAMMPS.
This is effectively a simple way to add per-atom properties to a model This is effectively a simple way to add per-atom properties to a model
without needing to write code for a new :doc:`atom style <atom_style>` without needing to write code for a new :doc:`atom style <atom_style>`
@ -107,42 +117,38 @@ new properties are also defined for the ghost atoms.
.. note:: .. note::
If you use this command with the *mol*\ , *q* or *rmass* vectors, If you use the *mol*\ , *q* or *rmass* names, you most likely want
then you most likely want to set *ghost* yes, since these properties to set *ghost* yes, since these properties are stored with ghost
are stored with ghost atoms if you use an :doc:`atom_style <atom_style>` atoms if you use an :doc:`atom_style <atom_style>` that defines
that defines them, and many LAMMPS operations that use molecule IDs or them. Many LAMMPS operations that use molecule IDs or charge, such
charge, such as neighbor lists and pair styles, will expect ghost as neighbor lists and pair styles, will expect ghost atoms to have
atoms to have these values. LAMMPS will issue a warning it you define these values. LAMMPS will issue a warning it you define those
those vectors but do not set *ghost* yes. vectors but do not set *ghost* yes.
.. note:: .. note::
The properties for ghost atoms are not updated every timestep, The specified properties for ghost atoms are not updated every
but only once every few steps when neighbor lists are re-built. Thus timestep, but only once every few steps when neighbor lists are
the *ghost* keyword is suitable for static properties, like molecule re-built. Thus the *ghost* keyword is suitable for static
IDs, but not for dynamic properties that change every step. For the properties, like molecule IDs, but not for dynamic properties that
latter, the code you add to LAMMPS to change the properties will also change every step. For the latter, the code you add to LAMMPS to
need to communicate their new values to/from ghost atoms, an operation change the properties will also need to communicate their new
that can be invoked from within a :doc:`pair style <pair_style>` or values to/from ghost atoms, an operation that can be invoked from
:doc:`fix <fix>` or :doc:`compute <compute>` that you write. within a :doc:`pair style <pair_style>` or :doc:`fix <fix>` or
:doc:`compute <compute>` that you write.
.. note::
If this fix is defined **after** the simulation box is created,
a 'run 0' command should be issued to properly initialize the storage
created by this fix.
---------- ----------
This fix is one of a small number that can be defined in an input This fix is one of a small number that can be defined in an input
script before the simulation box is created or atoms are defined. script before the simulation box is created or atoms are defined.
This is so it can be used with the :doc:`read_data <read_data>` command This is so it can be used with the :doc:`read_data <read_data>`
as described below. command as described next.
Per-atom properties that are defined by the :doc:`atom style <atom_style>` are initialized when atoms are created, e.g. by Per-atom properties that are defined by the :doc:`atom style
the :doc:`read_data <read_data>` or :doc:`create_atoms <create_atoms>` <atom_style>` are initialized when atoms are created, e.g. by the
:doc:`read_data <read_data>` or :doc:`create_atoms <create_atoms>`
commands. The per-atom properties defined by this fix are not. So commands. The per-atom properties defined by this fix are not. So
you need to initialize them explicitly. This can be done by the you need to initialize them explicitly. One way to do thisis
:doc:`read_data <read_data>` command, using its *fix* keyword and :doc:`read_data <read_data>` command, using its *fix* keyword and
passing it the fix-ID of this fix. passing it the fix-ID of this fix.
@ -167,15 +173,22 @@ would allow a data file to have a section like this:
... ...
N 763 4.5 N 763 4.5
where N is the number of atoms, and the first field on each line is where N is the number of atoms, the first field on each line is the
the atom-ID, followed by a molecule-ID and a floating point value that atom-ID, the next two are a molecule-ID and a floating point value
will be stored in a new property called "flag". Note that the list of that will be stored in a new property called "flag". If a per-atom
per-atom properties can be in any order. array was specified in the fix property/atom commmand then the *N*
values for that array must be specified consecutively for that
property on each line.
Another way of initializing the new properties is via the Note that the the lines of per-atom properties can be listed in any
:doc:`set <set>` command. For example, if you wanted molecules order. Also note that all the per-atom properties specified by the
defined for every set of 10 atoms, based on their atom-IDs, fix ID (prop in this case) must be included on each line in the
these commands could be used: specified data file section (Molecules in this case).
Another way of initializing the new properties is via the :doc:`set
<set>` command. For example, if you wanted molecules defined for
every set of 10 atoms, based on their atom-IDs, these commands could
be used:
.. code-block:: LAMMPS .. code-block:: LAMMPS
@ -185,53 +198,59 @@ these commands could be used:
The :doc:`atom-style variable <variable>` will create values for atoms The :doc:`atom-style variable <variable>` will create values for atoms
with IDs 31,32,33,...40 that are 4.0,4.1,4.2,...,4.9. When the with IDs 31,32,33,...40 that are 4.0,4.1,4.2,...,4.9. When the
:doc:`set <set>` commands assigns them to the molecule ID for each atom, :doc:`set <set>` commands assigns them to the molecule ID for each
they will be truncated to an integer value, so atoms 31-40 will all be atom, they will be truncated to an integer value, so atoms 31-40 will
assigned a molecule ID of 4. all be assigned a molecule ID of 4.
Note that :doc:`atomfile-style variables <variable>` can also be used in Note that :doc:`atomfile-style variables <variable>` can also be used
place of atom-style variables, which means in this case that the in place of atom-style variables, which means in this case that the
molecule IDs could be read-in from a separate file and assigned by the molecule IDs could be read-in from a separate file and assigned by the
:doc:`set <set>` command. This allows you to initialize new per-atom :doc:`set <set>` command. This allows you to initialize new per-atom
properties in a completely general fashion. properties in a completely general fashion.
---------- ----------
For new atom properties specified as *i_name* or *d_name*, the For new atom properties specified as *i_name*, *d_name*, *i2_name*, or
:doc:`compute property/atom <compute_property_atom>` command can access *d2_name*, the :doc:`dump custom <dump>` and :doc:`compute
their values. This means that the values can be output via the :doc:`dump custom <dump>` command, accessed by fixes like :doc:`fix ave/atom <fix_ave_atom>`, accessed by other computes like :doc:`compute reduce <compute_reduce>`, or used in :doc:`atom-style variables <variable>`. property/atom <compute_property_atom>` commands can access their
values. This means that the values can be used accessed by fixes like
:doc:`fix ave/atom <fix_ave_atom>`, accessed by other computes like
:doc:`compute reduce <compute_reduce>`, or used in :doc:`atom-style
variables <variable>`.
For example, these commands will output two new properties to a custom For example, these commands will output both the instantanous and
dump file: time-averaged values of two new properties to a custom dump file:
.. code-block:: LAMMPS .. code-block:: LAMMPS
fix prop all property/atom i_flag1 d_flag2 fix myprops all property/atom i_flag1 d_flag2
compute 1 all property/atom i_flag1 d_flag2 compute 1 all property/atom i_flag1 d_flag2
dump 1 all custom 100 tmp.dump id x y z c_1[1] c_1[2] fix 1 all ave/atom 10 10 100 c_1[1] c_1[2]
dump 1 all custom 100 tmp.dump id x y z i_flag1 d_flag2 f_1[1] f_1[2]
---------- ----------
If you wish to add new :doc:`pair styles <pair_style>`, If you wish to add new :doc:`pair styles <pair_style>`, :doc:`fixes
:doc:`fixes <fix>`, or :doc:`computes <compute>` that use the per-atom <fix>`, or :doc:`computes <compute>` that use the per-atom properties
properties defined by this fix, see the :doc:`Modify atom <Modify_atom>` defined by this fix, see the :doc:`Modify atom <Modify_atom>` doc page
doc page which has details on how the properties can be accessed from which has details on how the custom properties of this fix can be
added classes. accessed from added classes.
---------- ----------
.. _isotopes: .. _isotopes:
Example for using per-atom masses with TIP4P water to Here is an example of using per-atom masses with TIP4P water to study
study isotope effects. When setting up simulations with the :doc:`TIP4P pair styles <Howto_tip4p>` for water, you have to provide exactly isotope effects. When setting up simulations with the :doc:`TIP4P pair
one atom type each to identify the water oxygen and hydrogen styles <Howto_tip4p>` for water, you have to provide exactly one atom
atoms. Since the atom mass is normally tied to the atom type, this type each to identify the water oxygen and hydrogen atoms. Since the
makes it impossible to study multiple isotopes in the same simulation. atom mass is normally tied to the atom type, this makes it impossible
With *fix property/atom rmass* however, the per-type masses are to study multiple isotopes in the same simulation. With *fix
replaced by per-atom masses. Asumming you have a working input deck property/atom rmass* however, the per-type masses are replaced by
for regular TIP4P water, where water oxygen is atom type 1 and water per-atom masses. Asumming you have a working input deck for regular
hydrogen is atom type 2, the following lines of input script convert TIP4P water, where water oxygen is atom type 1 and water hydrogen is
this to using per-atom masses: atom type 2, the following lines of input script convert this to using
per-atom masses:
.. code-block:: LAMMPS .. code-block:: LAMMPS
@ -239,22 +258,22 @@ this to using per-atom masses:
set type 1 mass 15.9994 set type 1 mass 15.9994
set type 2 mass 1.008 set type 2 mass 1.008
When writing out the system data with the :doc:`write_data <write_data>` When writing out the system data with the :doc:`write_data
command, there will be a new section named with the fix-ID <write_data>` command, there will be a new section named with the
(i.e. *Isotopes* in this case). Alternatively, you can take an fix-ID (i.e. *Isotopes* in this case). Alternatively, you can take an
existing data file and just add this *Isotopes* section with existing data file and just add this *Isotopes* section with one line
one line per atom containing atom-ID and mass. Either way, the per atom containing atom-ID and mass. Either way, the extended data
extended data file can be read back with: file can be read back with:
.. code-block:: LAMMPS .. code-block:: LAMMPS
fix Isotopes all property/atom rmass ghost yes fix Isotopes all property/atom rmass ghost yes
read_data tip4p-isotopes.data fix Isotopes NULL Isotopes read_data tip4p-isotopes.data fix Isotopes NULL Isotopes
Please note that the first *Isotopes* refers to the fix-ID Please note that the first *Isotopes* refers to the fix-ID and the
and the second to the name of the section. The following input second to the name of the section. The following input script code
script code will now change the first 100 water molecules in this will now change the first 100 water molecules in this example to heavy
example to heavy water: water:
.. code-block:: LAMMPS .. code-block:: LAMMPS
@ -271,17 +290,19 @@ example to heavy water:
Restart, fix_modify, output, run start/stop, minimize info Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the per-atom values it stores to :doc:`binary restart files <restart>`, so that the values can be restored when a This fix writes the per-atom values it stores to :doc:`binary restart
simulation is restarted. See the :doc:`read_restart <read_restart>` files <restart>`, so that the values can be restored when a simulation
command for info on how to re-specify a fix in an input script that is restarted. See the :doc:`read_restart <read_restart>` command for
reads a restart file, so that the operation of the fix continues in an info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion. uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this None of the :doc:`fix_modify <fix_modify>` options are relevant to
fix. No global or per-atom quantities are stored by this fix for this fix. No global or per-atom quantities are stored by this fix for
access by various :doc:`output commands <Howto_output>`. No parameter access by various :doc:`output commands <Howto_output>`. No parameter
of this fix can be used with the *start/stop* keywords of the 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>`. :doc:`run <run>` command. This fix is not invoked during :doc:`energy
minimization <minimize>`.
Restrictions Restrictions
"""""""""""" """"""""""""
@ -290,9 +311,10 @@ Restrictions
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`read_data <read_data>`, :doc:`set <set>`, :doc:`compute property/atom <compute_property_atom>` :doc:`read_data <read_data>`, :doc:`set <set>`,
:doc:`compute property/atom <compute_property_atom>`
Default Default
""""""" """""""
The default keyword values are ghost = no. The default keyword value is ghost = no.

View File

@ -23,8 +23,8 @@ Syntax
q, mux, muy, muz, mu, q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz, radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz, angmomx, angmomy, angmomz, tqx, tqy, tqz,
c_ID, c_ID[N], f_ID, f_ID[N], v_name, c_ID, c_ID[I], f_ID, f_ID[I], v_name,
d_name, i_name d_name, i_name, i2_name[I], d2_name[I],
.. parsed-literal:: .. parsed-literal::
@ -46,13 +46,15 @@ Syntax
omegax,omegay,omegaz = angular velocity of spherical particle omegax,omegay,omegaz = angular velocity of spherical particle
angmomx,angmomy,angmomz = angular momentum of aspherical particle angmomx,angmomy,angmomz = angular momentum of aspherical particle
tqx,tqy,tqz = torque on finite-size particles tqx,tqy,tqz = torque on finite-size particles
c_ID = per-atom vector calculated by a compute with ID *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 *c_ID[I]* = Ith column of per-atom array calculated by a compute with ID
f_ID = per-atom vector calculated by a fix with ID *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 *f_ID[I]* = Ith column of per-atom array calculated by a fix with ID
v_name = per-atom vector calculated by an atom-style variable with name *v_name* = per-atom vector calculated by an atom-style variable with name
d_name = per-atom floating point vector name, managed by fix property/atom *i_name* = custom integer vector with name
i_name = per-atom integer vector name, managed by fix property/atom *d_name* = custom floating point vector with name
*i2_name[I]* = Ith column of custom integer array with name
*d2_name[I]* = Ith column of custom floating-point array with name
* zero or more keyword/value pairs may be appended * zero or more keyword/value pairs may be appended
* keyword = *com* * keyword = *com*
@ -92,7 +94,8 @@ steps.
those attributes may require quantities that are not defined in those attributes may require quantities that are not defined in
between runs. between runs.
The list of possible attributes is the same as that used by the :doc:`dump custom <dump>` command, which describes their meaning. The list of possible attributes is the same as that used by the
:doc:`dump custom <dump>` command, which describes their meaning.
If the *com* keyword is set to *yes* then the *xu*\ , *yu*\ , and *zu* If the *com* keyword is set to *yes* then the *xu*\ , *yu*\ , and *zu*
inputs store the position of each atom relative to the center-of-mass inputs store the position of each atom relative to the center-of-mass
@ -105,10 +108,11 @@ group.
Restart, fix_modify, output, run start/stop, minimize info Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the per-atom values it stores to :doc:`binary restart files <restart>`, so that the values can be restored when a This fix writes the per-atom values it stores to :doc:`binary restart
simulation is restarted. See the :doc:`read_restart <read_restart>` files <restart>`, so that the values can be restored when a simulation
command for info on how to re-specify a fix in an input script that is restarted. See the :doc:`read_restart <read_restart>` command for
reads a restart file, so that the operation of the fix continues in an info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion. uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this None of the :doc:`fix_modify <fix_modify>` options are relevant to this
@ -121,7 +125,8 @@ can be accessed by various :doc:`output commands <Howto_output>`. The
per-atom values be accessed on any timestep. per-atom values be accessed on any timestep.
No parameter of this fix can be used with the *start/stop* keywords of 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>`. the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions Restrictions
"""""""""""" """"""""""""

View File

@ -41,7 +41,7 @@ Syntax
keyword = *region* or *var* or *every* keyword = *region* or *var* or *every*
*region* value = region-ID *region* value = region-ID
*var* value = name of variable *var* value = name of variable
*property* value = name of per-atom property *property* value = name of custom integer or floating point vector
*every* value = N = update group every this many timesteps *every* value = N = update group every this many timesteps
*static* = no args *static* = no args
@ -226,18 +226,33 @@ simulation runs. This is in contrast to static groups where atoms are
permanently assigned to the group. The way the assignment occurs is permanently assigned to the group. The way the assignment occurs is
as follows. Only atoms in the group specified as the parent group via as follows. Only atoms in the group specified as the parent group via
the parent-ID are assigned to the dynamic group before the following the parent-ID are assigned to the dynamic group before the following
conditions are applied. If the *region* keyword is used, atoms not in conditions are applied.
the specified region are removed from the dynamic group. If the *var*
keyword is used, the variable name must be an atom-style or If the *region* keyword is used, atoms not in the specified region are
atomfile-style variable. The variable is evaluated and atoms whose removed from the dynamic group.
per-atom values are 0.0, are removed from the dynamic group. If the *property*
keyword is used, the per-atom property name must be a previously defined If the *var* keyword is used, the variable name must be an atom-style
per-atom property. The per-atom property is evaluated and atoms whose or atomfile-style variable. The variable is evaluated and atoms whose
values are 0.0 are removed from the dynamic group. per-atom values are 0.0, are removed from the dynamic group.
If the *property* keyword is used, the name refers to a custom integer
or floating point per-atom vector defined via the :doc:`fix
property/atom <fix_property_atom>` command. This means the values in
the vector can be read as part of a data file with the :doc:`read_data
<read_data>` command or specified with the :doc:`set <set>` command.
Or accessed and changed via the :doc:`library interface to LAMMPS
<Howto_library>`, or by styles you add to LAMMPS (pair, fix, compute,
etc) which access the custom vector and modify its values. Which
means the values can be modified between or during simulations. Atoms
whose values in the custom vector are zero are removed from the
dynamic group. Note that the name of the custom per-atom vector is
specified just as *name*, not as *i_name* or *d_name* as it is for
other commands that use different kinds of custom atom vectors or
arrays as arguments.
The assignment of atoms to a dynamic group is done at the beginning of The assignment of atoms to a dynamic group is done at the beginning of
each run and on every timestep that is a multiple of *N*\ , which is the each run and on every timestep that is a multiple of *N*\ , which is
argument for the *every* keyword (N = 1 is the default). For an the argument for the *every* keyword (N = 1 is the default). For an
energy minimization, via the :doc:`minimize <minimize>` command, an energy minimization, via the :doc:`minimize <minimize>` command, an
assignment is made at the beginning of the minimization, but not assignment is made at the beginning of the minimization, but not
during the iterations of the minimizer. during the iterations of the minimizer.

View File

@ -13,7 +13,17 @@ Syntax
* style = *atom* or *type* or *mol* or *group* or *region* * style = *atom* or *type* or *mol* or *group* or *region*
* ID = atom ID range or type range or mol ID range or group ID or region ID * ID = atom ID range or type range or mol ID range or group ID or region ID
* one or more keyword/value pairs may be appended * one or more keyword/value pairs may be appended
* keyword = *type* or *type/fraction* or *type/ratio* or *type/subset* or *mol* or *x* or *y* or *z* or *charge* or *dipole* or *dipole/random* or *quat* or *spin* or *spin/random* or *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or *theta* or *theta/random* or *angmom* or *omega* or *mass* or *density* or *density/disc* or *volume* or *image* or *bond* or *angle* or *dihedral* or *improper* or *sph/e* or *sph/cv* or *sph/rho* or *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or *edpd/temp* or *edpd/cv* or *cc* or *i_name* or *d_name* * keyword = *type* or *type/fraction* or *type/ratio* or *type/subset*
or *mol* or *x* or *y* or *z* or *charge* or *dipole* or
*dipole/random* or *quat* or *spin* or *spin/random* or
*quat* or *quat/random* or *diameter* or *shape* or
*length* or *tri* or *theta* or *theta/random* or *angmom* or
*omega* or *mass* or *density* or *density/disc* or
*volume* or *image* or *bond* or *angle* or *dihedral* or
*improper* or *sph/e* or *sph/cv* or *sph/rho* or
*smd/contact/radius* or *smd/mass/density* or *dpd/theta* or
*edpd/temp* or *edpd/cv* or *cc* or
*i_name* or *d_name* or *i2_name* or *d2_name*
.. parsed-literal:: .. parsed-literal::
@ -114,8 +124,12 @@ Syntax
*cc* values = index cc *cc* values = index cc
index = index of a chemical species (1 to Nspecies) index = index of a chemical species (1 to Nspecies)
cc = chemical concentration of tDPD particles for a species (mole/volume units) cc = chemical concentration of tDPD particles for a species (mole/volume units)
*i_name* value = value for custom integer vector with name *i_name* value = custom integer vector with name
*d_name* value = value for custom floating-point vector with name *d_name* value = custom floating-point vector with name
*i2_name* value = column of a custom integer array with name
column specified as i2_name[N] where N is 1 to Ncol
*d2_name* value = column of a custom floating-point array with name
column specified as d2_name[N] where N is 1 to Ncol
Examples Examples
"""""""" """"""""
@ -132,6 +146,8 @@ Examples
set atom 100*200 x 0.5 y 1.0 set atom 100*200 x 0.5 y 1.0
set atom 100 vx 0.0 vy 0.0 vz -1.0 set atom 100 vx 0.0 vy 0.0 vz -1.0
set atom 1492 type 3 set atom 1492 type 3
set atom * i_myVal 5
set atom * d2_Sxyz[1] 6.4
Description Description
""""""""""" """""""""""
@ -473,11 +489,13 @@ attribute. An integer for "index" selects a chemical species (1 to
Nspecies) where Nspecies is set by the atom_style command. The value Nspecies) where Nspecies is set by the atom_style command. The value
for the chemical concentration must be >= 0.0. for the chemical concentration must be >= 0.0.
Keywords *i_name* and *d_name* refer to custom integer and Keywords *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom
floating-point properties that have been added to each atom via the per-atom integer and floating-point vectors or arrays that have been
:doc:`fix property/atom <fix_property_atom>` command. When that command added via the :doc:`fix property/atom <fix_property_atom>` command.
is used specific names are given to each attribute which are what is When that command is used specific names are given to each attribute
specified as the "name" portion of *i_name* or *d_name*. which are the "name" portion of these keywords. For arrays *i2_name*
and *d2_name*, the column of the array must also be included following
the name in brackets: e.g. d2_xyz[2], i2_mySpin[3].
Restrictions Restrictions
"""""""""""" """"""""""""

View File

@ -260,6 +260,7 @@ clean:
clean-all: clean-all:
rm -rf Obj_* rm -rf Obj_*
rm style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h
clean-%: clean-%:
@if [ $@ = "clean-serial" ]; \ @if [ $@ = "clean-serial" ]; \

View File

@ -178,10 +178,13 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// custom atom arrays // custom atom arrays
nivector = ndvector = 0; nivector = ndvector = niarray = ndarray = 0;
ivector = NULL; ivector = NULL;
dvector = NULL; dvector = NULL;
iname = dname = NULL; iarray = NULL;
darray = NULL;
icols = dcols = NULL;
ivname = dvname = ianame = daname = NULL;
// initialize atom style and array existence flags // initialize atom style and array existence flags
@ -262,20 +265,32 @@ Atom::~Atom()
// delete custom atom arrays // delete custom atom arrays
for (int i = 0; i < nivector; i++) { for (int i = 0; i < nivector; i++) {
delete [] iname[i]; delete [] ivname[i];
memory->destroy(ivector[i]); memory->destroy(ivector[i]);
} }
if (dvector != NULL) { for (int i = 0; i < ndvector; i++) {
for (int i = 0; i < ndvector; i++) { delete [] dvname[i];
delete [] dname[i]; memory->destroy(dvector[i]);
memory->destroy(dvector[i]); }
} for (int i = 0; i < niarray; i++) {
delete [] ianame[i];
memory->destroy(iarray[i]);
}
for (int i = 0; i < ndarray; i++) {
delete [] daname[i];
memory->destroy(darray[i]);
} }
memory->sfree(iname); memory->sfree(ivname);
memory->sfree(dname); memory->sfree(dvname);
memory->sfree(ianame);
memory->sfree(daname);
memory->sfree(ivector); memory->sfree(ivector);
memory->sfree(dvector); memory->sfree(dvector);
memory->sfree(iarray);
memory->sfree(darray);
memory->sfree(icols);
memory->sfree(dcols);
// delete user-defined molecules // delete user-defined molecules
@ -2282,23 +2297,41 @@ void Atom::update_callback(int ifix)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
find custom per-atom vector with name find custom per-atom vector with name
return index if found, and flag = 0/1 for int/double return index if found, -1 if not found
return -1 if not found lists of names can have NULL entries if previously removed
return flag = 0/1 for int/double
return cols = 0/N for vector/array where N = # of columns
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int Atom::find_custom(const char *name, int &flag) int Atom::find_custom(const char *name, int &flag, int &cols)
{ {
if(name == NULL) return -1; if (name == NULL) return -1;
for (int i = 0; i < nivector; i++) for (int i = 0; i < nivector; i++)
if (iname[i] && strcmp(iname[i],name) == 0) { if (ivname[i] && strcmp(ivname[i],name) == 0) {
flag = 0; flag = 0;
cols = 0;
return i; return i;
} }
for (int i = 0; i < ndvector; i++) for (int i = 0; i < ndvector; i++)
if (dname[i] && strcmp(dname[i],name) == 0) { if (dvname[i] && strcmp(dvname[i],name) == 0) {
flag = 1; flag = 1;
cols = 0;
return i;
}
for (int i = 0; i < niarray; i++)
if (ianame[i] && strcmp(ianame[i],name) == 0) {
flag = 0;
cols = icols[i];
return i;
}
for (int i = 0; i < ndarray; i++)
if (daname[i] && strcmp(daname[i],name) == 0) {
flag = 1;
cols = dcols[i];
return i; return i;
} }
@ -2306,60 +2339,107 @@ int Atom::find_custom(const char *name, int &flag)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add a custom variable with name of type flag = 0/1 for int/double add a custom variable with name
flag = 0/1 for int/double, cols = 0/N for vector/array where N = # of columns
assumes name does not already exist assumes name does not already exist
return index in ivector or dvector of its location vectors of names and data ptrs are always incremented by one
return index in vector/array lists of its location
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int Atom::add_custom(const char *name, int flag) int Atom::add_custom(const char *name, int flag, int cols)
{ {
int index; int index;
if (flag == 0) { if (flag == 0 && cols == 0) {
index = nivector; index = nivector;
nivector++; nivector++;
iname = (char **) memory->srealloc(iname,nivector*sizeof(char *), ivname = (char **) memory->srealloc(ivname,nivector*sizeof(char *),
"atom:iname"); "atom:ivname");
int n = strlen(name) + 1; int n = strlen(name) + 1;
iname[index] = new char[n]; ivname[index] = new char[n];
strcpy(iname[index],name); strcpy(ivname[index],name);
ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *), ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
"atom:ivector"); "atom:ivector");
memory->create(ivector[index],nmax,"atom:ivector"); memory->create(ivector[index],nmax,"atom:ivector");
} else {
} else if (flag == 1 && cols == 0) {
index = ndvector; index = ndvector;
ndvector++; ndvector++;
dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *), dvname = (char **) memory->srealloc(dvname,ndvector*sizeof(char *),
"atom:dname"); "atom:dvname");
int n = strlen(name) + 1; int n = strlen(name) + 1;
dname[index] = new char[n]; dvname[index] = new char[n];
strcpy(dname[index],name); strcpy(dvname[index],name);
dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *), dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),
"atom:dvector"); "atom:dvector");
memory->create(dvector[index],nmax,"atom:dvector"); memory->create(dvector[index],nmax,"atom:dvector");
}
} else if (flag == 0 && cols) {
index = niarray;
niarray++;
ianame = (char **) memory->srealloc(ianame,niarray*sizeof(char *),
"atom:ianame");
int n = strlen(name) + 1;
ianame[index] = new char[n];
strcpy(ianame[index],name);
iarray = (int ***) memory->srealloc(iarray,niarray*sizeof(int **),
"atom:iarray");
memory->create(iarray[index],nmax,cols,"atom:iarray");
icols = (int *) memory->srealloc(icols,niarray*sizeof(int),"atom:icols");
icols[index] = cols;
} else if (flag == 1 && cols) {
index = ndarray;
ndarray++;
daname = (char **) memory->srealloc(daname,ndarray*sizeof(char *),
"atom:daname");
int n = strlen(name) + 1;
daname[index] = new char[n];
strcpy(daname[index],name);
darray = (double ***) memory->srealloc(darray,ndarray*sizeof(double **),
"atom:darray");
memory->create(darray[index],nmax,cols,"atom:darray");
dcols = (int *) memory->srealloc(dcols,ndarray*sizeof(int),"atom:dcols");
dcols[index] = cols;
}
return index; return index;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
remove a custom variable of type flag = 0/1 for int/double at index remove a custom variable stored at index
free memory for vector and name and set ptrs to NULL flag = 0/1 for int/double, cols = 0/1 for vector/array or N for array columns
ivector/dvector and iname/dname lists never shrink free memory for vector/array and name and set ptrs to NULL
vector/array lists and name lists never shrink
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Atom::remove_custom(int flag, int index) void Atom::remove_custom(int index, int flag, int cols)
{ {
if (flag == 0) { if (flag == 0 && cols == 0) {
memory->destroy(ivector[index]); memory->destroy(ivector[index]);
ivector[index] = NULL; ivector[index] = NULL;
delete [] iname[index]; delete [] ivname[index];
iname[index] = NULL; ivname[index] = NULL;
} else {
} else if (flag == 1 && cols == 0) {
memory->destroy(dvector[index]); memory->destroy(dvector[index]);
dvector[index] = NULL; dvector[index] = NULL;
delete [] dname[index]; delete [] dvname[index];
dname[index] = NULL; dvname[index] = NULL;
} else if (flag == 0 && cols) {
memory->destroy(iarray[index]);
iarray[index] = NULL;
delete [] ianame[index];
ianame[index] = NULL;
} else if (flag == 1 && cols) {
memory->destroy(darray[index]);
darray[index] = NULL;
delete [] daname[index];
daname[index] = NULL;
} }
} }
@ -2372,7 +2452,8 @@ void *Atom::extract(char *name)
{ {
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// 4th customization section: customize by adding new variable name // 4th customization section: customize by adding new variable name
// if new variable is from a package, add package comment
if (strcmp(name,"mass") == 0) return (void *) mass; if (strcmp(name,"mass") == 0) return (void *) mass;
if (strcmp(name,"id") == 0) return (void *) tag; if (strcmp(name,"id") == 0) return (void *) tag;
@ -2415,10 +2496,30 @@ void *Atom::extract(char *name)
if (strcmp(name,"cv") == 0) return (void *) cv; if (strcmp(name,"cv") == 0) return (void *) cv;
if (strcmp(name,"vest") == 0) return (void *) vest; if (strcmp(name,"vest") == 0) return (void *) vest;
// USER-MESONT package // custom vectors and arrays
if (strcmp(name,"length") == 0) return (void *) length;
if (strcmp(name,"buckling") == 0) return (void *) buckling; if (strstr(name,"i_") == name || strstr(name,"d_") == name ||
if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt; strstr(name,"i2_") == name || strstr(name,"d2_") == name) {
int which = 0;
if (name[0] == 'd') which = 1;
int array = 0;
if (name[1] == '2') array = 1;
int index,flag,cols;
if (!array) index = find_custom(&name[2],flag,cols);
else index = find_custom(&name[3],flag,cols);
if (index < 0) return NULL;
if (which != flag) return NULL;
if ((!array && cols) || (array && !cols)) return NULL;
if (!which && !array) return (void *) ivector[index];
if (which && !array) return (void *) dvector[index];
if (!which && array) return (void *) iarray[index];
if (which && array) return (void *) darray[index];
}
// USER-SMD package
if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius; if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius;
if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9; if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9;
@ -2429,9 +2530,20 @@ void *Atom::extract(char *name)
return (void *) eff_plastic_strain_rate; return (void *) eff_plastic_strain_rate;
if (strcmp(name, "damage") == 0) return (void *) damage; if (strcmp(name, "damage") == 0) return (void *) damage;
// USER-DPD package
if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta; if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
// USER-MESO package
if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp; if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp;
// USER-MESONT package
if (strcmp(name,"length") == 0) return (void *) length;
if (strcmp(name,"buckling") == 0) return (void *) buckling;
if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt;
// end of customization section // end of customization section
// -------------------------------------------------------------------- // --------------------------------------------------------------------

View File

@ -206,12 +206,13 @@ class Atom : protected Pointers {
PerAtom *peratom; PerAtom *peratom;
int nperatom,maxperatom; int nperatom,maxperatom;
// custom arrays used by fix property/atom // custom vectors and arrays used by fix property/atom
int **ivector; int **ivector,***iarray;
double **dvector; double **dvector,***darray;
char **iname,**dname; int *icols,*dcols;
int nivector,ndvector; char **ivname,**dvname,**ianame,**daname;
int nivector,ndvector,niarray,ndarray;
// molecule templates // molecule templates
// each template can be a set of consecutive molecules // each template can be a set of consecutive molecules
@ -321,9 +322,9 @@ class Atom : protected Pointers {
void delete_callback(const char *, int); void delete_callback(const char *, int);
void update_callback(int); void update_callback(int);
int find_custom(const char *, int &); int find_custom(const char *, int &, int &);
virtual int add_custom(const char *, int); virtual int add_custom(const char *, int, int);
virtual void remove_custom(int, int); virtual void remove_custom(int, int, int);
virtual void sync_modify(ExecutionSpace, unsigned int, unsigned int) {} virtual void sync_modify(ExecutionSpace, unsigned int, unsigned int) {}

View File

@ -24,6 +24,7 @@
#include "update.h" #include "update.h"
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "utils.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
@ -33,7 +34,7 @@ using namespace LAMMPS_NS;
ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg),
index(NULL), pack_choice(NULL) index(NULL), colindex(NULL), pack_choice(NULL)
{ {
if (narg < 4) error->all(FLERR,"Illegal compute property/atom command"); if (narg < 4) error->all(FLERR,"Illegal compute property/atom command");
@ -47,6 +48,7 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
pack_choice = new FnPtrPack[nvalues]; pack_choice = new FnPtrPack[nvalues];
index = new int[nvalues]; index = new int[nvalues];
colindex = new int[nvalues];
int i; int i;
for (int iarg = 3; iarg < narg; iarg++) { for (int iarg = 3; iarg < narg; iarg++) {
@ -141,7 +143,10 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Compute property/atom for " error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated"); "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_mu; pack_choice[i] = &ComputePropertyAtom::pack_mu;
} else if (strcmp(arg[iarg],"spx") == 0) { // pack magnetic variables
// pack magnetic variables
} else if (strcmp(arg[iarg],"spx") == 0) {
if (!atom->sp_flag) if (!atom->sp_flag)
error->all(FLERR,"Compute property/atom for " error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated"); "atom property that isn't allocated");
@ -176,6 +181,17 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Compute property/atom for " error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated"); "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_fmz; pack_choice[i] = &ComputePropertyAtom::pack_fmz;
// bond count
} else if (strcmp(arg[iarg],"nbonds") == 0) {
if (!atom->molecule_flag)
error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_nbonds;
// finite-size particles
} else if (strcmp(arg[iarg],"radius") == 0) { } else if (strcmp(arg[iarg],"radius") == 0) {
if (!atom->radius_flag) if (!atom->radius_flag)
error->all(FLERR,"Compute property/atom for " error->all(FLERR,"Compute property/atom for "
@ -329,6 +345,7 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (!avec_tri) error->all(FLERR,"Compute property/atom for " if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated"); "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner2x; pack_choice[i] = &ComputePropertyAtom::pack_corner2x;
} else if (strcmp(arg[iarg],"corner2y") == 0) { } else if (strcmp(arg[iarg],"corner2y") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri"); avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for " if (!avec_tri) error->all(FLERR,"Compute property/atom for "
@ -355,40 +372,59 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
"atom property that isn't allocated"); "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner3z; pack_choice[i] = &ComputePropertyAtom::pack_corner3z;
} else if (strcmp(arg[iarg],"nbonds") == 0) { // custom per-atom vectors
if (!atom->molecule_flag)
error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_nbonds;
} else if (strstr(arg[iarg],"i_") == arg[iarg]) { } else if (strstr(arg[iarg],"i_") == arg[iarg]) {
int flag; int flag,cols;
index[i] = atom->find_custom(&arg[iarg][2],flag); index[i] = atom->find_custom(&arg[iarg][2],flag,cols);
if (index[i] < 0 || flag != 0) if (index[i] < 0 || flag || cols)
error->all(FLERR,"Compute property/atom integer " error->all(FLERR,"Compute property/atom custom vector does not exist");
"vector does not exist");
pack_choice[i] = &ComputePropertyAtom::pack_iname; pack_choice[i] = &ComputePropertyAtom::pack_iname;
} else if (strstr(arg[iarg],"d_") == arg[iarg]) { } else if (strstr(arg[iarg],"d_") == arg[iarg]) {
int flag; int flag,cols;
index[i] = atom->find_custom(&arg[iarg][2],flag); index[i] = atom->find_custom(&arg[iarg][2],flag,cols);
if (index[i] < 0 || flag != 1) if (index[i] < 0 || !flag || cols)
error->all(FLERR,"Compute property/atom floating point " error->all(FLERR,"Compute property/atom custom vector does not exist");
"vector does not exist");
pack_choice[i] = &ComputePropertyAtom::pack_dname; pack_choice[i] = &ComputePropertyAtom::pack_dname;
}
else if (strcmp(arg[iarg],"buckling") == 0) { // custom per-atom arrays, must include bracketed index
if (!atom->mesont_flag)
error->all(FLERR,"Compute property/atom for " } else if (strstr(arg[iarg],"i2_") == arg[iarg] ||
"atom property that isn't allocated"); strstr(arg[iarg],"d2_") == arg[iarg]) {
pack_choice[i] = &ComputePropertyAtom::pack_buckling; int which = 0;
// check if atom style recognizes keyword if (arg[iarg][0] == 'd') which = 1;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][3]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in set command");
suffix[strlen(suffix)-1] = '\0';
colindex[i] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0';
} else error->all(FLERR,"Compute property/atom custom array is not indexed");
int flag,cols;
index[i] = atom->find_custom(suffix,flag,cols);
delete [] suffix;
if ((!which && (index[i] < 0 || flag || !cols)) ||
(which && (index[i] < 0 || !flag || !cols)))
error->all(FLERR,"Compute property/atom custom array does not exist");
if (!which) pack_choice[i] = &ComputePropertyAtom::pack_i2name;
else pack_choice[i] = &ComputePropertyAtom::pack_d2name;
// anything else must be recognized by atom style
} else { } else {
index[i] = atom->avec->property_atom(arg[iarg]); index[i] = atom->avec->property_atom(arg[iarg]);
if (index[i] < 0) if (index[i] < 0)
error->all(FLERR,"Invalid keyword in compute property/atom command"); error->all(FLERR,"Invalid keyword in compute property/atom command");
pack_choice[i] = &ComputePropertyAtom::pack_property_atom; pack_choice[i] = &ComputePropertyAtom::pack_atom_style;
} }
} }
@ -401,6 +437,7 @@ ComputePropertyAtom::~ComputePropertyAtom()
{ {
delete [] pack_choice; delete [] pack_choice;
delete [] index; delete [] index;
delete [] colindex;
memory->destroy(vector_atom); memory->destroy(vector_atom);
memory->destroy(array_atom); memory->destroy(array_atom);
} }
@ -413,6 +450,9 @@ void ComputePropertyAtom::init()
avec_line = (AtomVecLine *) atom->style_match("line"); avec_line = (AtomVecLine *) atom->style_match("line");
avec_tri = (AtomVecTri *) atom->style_match("tri"); avec_tri = (AtomVecTri *) atom->style_match("tri");
avec_body = (AtomVecBody *) atom->style_match("body"); avec_body = (AtomVecBody *) atom->style_match("body");
// NOTE: could reset custom vector/array indices here, like dump custom does
// in case have been deleted
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -1145,6 +1185,21 @@ void ComputePropertyAtom::pack_fmz(int n)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_nbonds(int n)
{
int *num_bond = atom->num_bond;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = num_bond[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_radius(int n) void ComputePropertyAtom::pack_radius(int n)
{ {
double *radius = atom->radius; double *radius = atom->radius;
@ -1780,36 +1835,6 @@ void ComputePropertyAtom::pack_corner3z(int n)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_buckling(int n)
{
int *buckling = atom->buckling;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = static_cast<double>(buckling[i]);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_nbonds(int n)
{
int *num_bond = atom->num_bond;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = num_bond[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_iname(int n) void ComputePropertyAtom::pack_iname(int n)
{ {
int *ivector = atom->ivector[index[n]]; int *ivector = atom->ivector[index[n]];
@ -1840,7 +1865,39 @@ void ComputePropertyAtom::pack_dname(int n)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_property_atom(int n) void ComputePropertyAtom::pack_i2name(int n)
{
int **iarray = atom->iarray[index[n]];
int icol = colindex[n] - 1;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = iarray[i][icol];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_d2name(int n)
{
double **darray = atom->darray[index[n]];
int icol = colindex[n] - 1;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = darray[i][icol];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_atom_style(int n)
{ {
atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit); atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit);
} }

View File

@ -35,7 +35,7 @@ class ComputePropertyAtom : public Compute {
private: private:
int nvalues; int nvalues;
int nmax; int nmax;
int *index; int *index,*colindex;
double *buf; double *buf;
class AtomVecEllipsoid *avec_ellipsoid; class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line; class AtomVecLine *avec_line;
@ -81,10 +81,6 @@ class ComputePropertyAtom : public Compute {
void pack_muy(int); void pack_muy(int);
void pack_muz(int); void pack_muz(int);
void pack_mu(int); void pack_mu(int);
void pack_radius(int);
void pack_diameter(int);
// pack magnetic variables
void pack_spx(int); void pack_spx(int);
void pack_spy(int); void pack_spy(int);
@ -94,12 +90,17 @@ class ComputePropertyAtom : public Compute {
void pack_fmy(int); void pack_fmy(int);
void pack_fmz(int); void pack_fmz(int);
void pack_nbonds(int);
void pack_radius(int);
void pack_diameter(int);
void pack_omegax(int); void pack_omegax(int);
void pack_omegay(int); void pack_omegay(int);
void pack_omegaz(int); void pack_omegaz(int);
void pack_angmomx(int); void pack_angmomx(int);
void pack_angmomy(int); void pack_angmomy(int);
void pack_angmomz(int); void pack_angmomz(int);
void pack_shapex(int); void pack_shapex(int);
void pack_shapey(int); void pack_shapey(int);
void pack_shapez(int); void pack_shapez(int);
@ -110,12 +111,14 @@ class ComputePropertyAtom : public Compute {
void pack_tqx(int); void pack_tqx(int);
void pack_tqy(int); void pack_tqy(int);
void pack_tqz(int); void pack_tqz(int);
void pack_end1x(int); void pack_end1x(int);
void pack_end1y(int); void pack_end1y(int);
void pack_end1z(int); void pack_end1z(int);
void pack_end2x(int); void pack_end2x(int);
void pack_end2y(int); void pack_end2y(int);
void pack_end2z(int); void pack_end2z(int);
void pack_corner1x(int); void pack_corner1x(int);
void pack_corner1y(int); void pack_corner1y(int);
void pack_corner1z(int); void pack_corner1z(int);
@ -125,14 +128,13 @@ class ComputePropertyAtom : public Compute {
void pack_corner3x(int); void pack_corner3x(int);
void pack_corner3y(int); void pack_corner3y(int);
void pack_corner3z(int); void pack_corner3z(int);
void pack_buckling(int);
void pack_nbonds(int);
void pack_iname(int); void pack_iname(int);
void pack_dname(int); void pack_dname(int);
void pack_i2name(int);
void pack_d2name(int);
void pack_property_atom(int); void pack_atom_style(int);
}; };
} }

View File

@ -30,6 +30,7 @@
#include "update.h" #include "update.h"
#include "variable.h" #include "variable.h"
#include "fmt/format.h" #include "fmt/format.h"
#include "utils.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -44,7 +45,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER, Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
TQX,TQY,TQZ, TQX,TQY,TQZ,
COMPUTE,FIX,VARIABLE,INAME,DNAME}; COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY};
enum{LT,LE,GT,GE,EQ,NEQ,XOR}; enum{LT,LE,GT,GE,EQ,NEQ,XOR};
#define INVOKED_PERATOM 8 #define INVOKED_PERATOM 8
@ -60,7 +61,7 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
earg(NULL), vtype(NULL), vformat(NULL), columns(NULL), choose(NULL), earg(NULL), vtype(NULL), vformat(NULL), columns(NULL), choose(NULL),
dchoose(NULL), clist(NULL), field2index(NULL), argindex(NULL), id_compute(NULL), dchoose(NULL), clist(NULL), field2index(NULL), argindex(NULL), id_compute(NULL),
compute(NULL), id_fix(NULL), fix(NULL), id_variable(NULL), variable(NULL), compute(NULL), id_fix(NULL), fix(NULL), id_variable(NULL), variable(NULL),
vbuf(NULL), id_custom(NULL), flag_custom(NULL), typenames(NULL), vbuf(NULL), id_custom(NULL), custom(NULL), custom_flag(NULL), typenames(NULL),
pack_choice(NULL) pack_choice(NULL)
{ {
if (narg == 5) error->all(FLERR,"No dump custom arguments specified"); if (narg == 5) error->all(FLERR,"No dump custom arguments specified");
@ -119,7 +120,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
ncustom = 0; ncustom = 0;
id_custom = NULL; id_custom = NULL;
flag_custom = NULL; custom = NULL;
custom_flag = NULL;
// process attributes // process attributes
// ioptional = start of additional optional args in expanded args // ioptional = start of additional optional args in expanded args
@ -236,7 +238,8 @@ DumpCustom::~DumpCustom()
for (int i = 0; i < ncustom; i++) delete [] id_custom[i]; for (int i = 0; i < ncustom; i++) delete [] id_custom[i];
memory->sfree(id_custom); memory->sfree(id_custom);
delete [] flag_custom; delete [] custom;
delete [] custom_flag;
memory->destroy(choose); memory->destroy(choose);
memory->destroy(dchoose); memory->destroy(dchoose);
@ -324,7 +327,7 @@ void DumpCustom::init_style()
else if (buffer_flag == 1) write_choice = &DumpCustom::write_string; else if (buffer_flag == 1) write_choice = &DumpCustom::write_string;
else write_choice = &DumpCustom::write_lines; else write_choice = &DumpCustom::write_lines;
// find current ptr for each compute,fix,variable // find current ptr for each compute,fix,variable and custom atom property
// check that fix frequency is acceptable // check that fix frequency is acceptable
int icompute; int icompute;
@ -353,9 +356,15 @@ void DumpCustom::init_style()
int icustom; int icustom;
for (int i = 0; i < ncustom; i++) { for (int i = 0; i < ncustom; i++) {
icustom = atom->find_custom(id_custom[i],flag_custom[i]); int flag,cols;
icustom = atom->find_custom(id_custom[i],flag,cols);
if (icustom < 0) if (icustom < 0)
error->all(FLERR,"Could not find custom per-atom property ID"); error->all(FLERR,"Could not find dump custom atom property name");
custom[i] = icustom;
if (!flag && !cols) custom_flag[i] = IVEC;
else if (flag && !cols) custom_flag[i] = DVEC;
else if (!flag && cols) custom_flag[i] = IARRAY;
else if (flag && cols) custom_flag[i] = DARRAY;
} }
// set index and check validity of region // set index and check validity of region
@ -1036,23 +1045,37 @@ int DumpCustom::count()
ptr = vbuf[field2index[i]]; ptr = vbuf[field2index[i]];
nstride = 1; nstride = 1;
} else if (thresh_array[ithresh] == DNAME) { } else if (thresh_array[ithresh] == IVEC) {
int iwhich,tmp;
i = nfield + ithresh; i = nfield + ithresh;
iwhich = atom->find_custom(id_custom[field2index[i]],tmp); int iwhich = custom[field2index[i]];
ptr = atom->dvector[iwhich];
nstride = 1;
} else if (thresh_array[ithresh] == INAME) {
int iwhich,tmp;
i = nfield + ithresh;
iwhich = atom->find_custom(id_custom[field2index[i]],tmp);
int *ivector = atom->ivector[iwhich]; int *ivector = atom->ivector[iwhich];
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
dchoose[i] = ivector[i]; dchoose[i] = ivector[i];
ptr = dchoose; ptr = dchoose;
nstride = 1; nstride = 1;
} else if (thresh_array[ithresh] == DVEC) {
i = nfield + ithresh;
int iwhich = custom[field2index[i]];
ptr = atom->dvector[iwhich];
nstride = 1;
} else if (thresh_array[ithresh] == IARRAY) {
i = nfield + ithresh;
int iwhich = custom[field2index[i]];
int **iarray = atom->iarray[iwhich];
int icol = argindex[i] - 1;
for (i = 0; i < nlocal; i++)
dchoose[i] = iarray[i][icol];
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == DARRAY) {
i = nfield + ithresh;
int iwhich = custom[field2index[i]];
double **darray = atom->darray[iwhich];
ptr = &darray[0][argindex[i]-1];
nstride = atom->dcols[iwhich];
} }
// unselect atoms that don't meet threshold criterion // unselect atoms that don't meet threshold criterion
@ -1257,7 +1280,7 @@ int DumpCustom::parse_fields(int narg, char **arg)
int i; int i;
for (int iarg = 0; iarg < narg; iarg++) { for (int iarg = 0; iarg < narg; iarg++) {
i = iarg; i = iarg;
if (strcmp(arg[iarg],"id") == 0) { if (strcmp(arg[iarg],"id") == 0) {
pack_choice[i] = &DumpCustom::pack_id; pack_choice[i] = &DumpCustom::pack_id;
if (sizeof(tagint) == sizeof(smallint)) vtype[i] = Dump::INT; if (sizeof(tagint) == sizeof(smallint)) vtype[i] = Dump::INT;
@ -1455,7 +1478,8 @@ int DumpCustom::parse_fields(int narg, char **arg)
if (ptr) { if (ptr) {
if (suffix[strlen(suffix)-1] != ']') if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in dump custom command"); error->all(FLERR,"Invalid attribute in dump custom command");
argindex[i] = atoi(ptr+1); suffix[strlen(suffix)-1] = '\0';
argindex[i] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0'; *ptr = '\0';
} else argindex[i] = 0; } else argindex[i] = 0;
@ -1491,7 +1515,8 @@ int DumpCustom::parse_fields(int narg, char **arg)
if (ptr) { if (ptr) {
if (suffix[strlen(suffix)-1] != ']') if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in dump custom command"); error->all(FLERR,"Invalid attribute in dump custom command");
argindex[i] = atoi(ptr+1); suffix[strlen(suffix)-1] = '\0';
argindex[i] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0'; *ptr = '\0';
} else argindex[i] = 0; } else argindex[i] = 0;
@ -1530,48 +1555,66 @@ int DumpCustom::parse_fields(int narg, char **arg)
field2index[i] = add_variable(suffix); field2index[i] = add_variable(suffix);
delete [] suffix; delete [] suffix;
// custom per-atom floating point value = d_ID // custom per-atom vector = i_ID or d_ID
} else if (strncmp(arg[iarg],"d_",2) == 0) { } else if (strncmp(arg[iarg],"i_",2) == 0 ||
strncmp(arg[iarg],"d_",2) == 0) {
int which = 0;
if (arg[iarg][0] == 'd') which = 1;
pack_choice[i] = &DumpCustom::pack_custom; pack_choice[i] = &DumpCustom::pack_custom;
vtype[i] = Dump::DOUBLE; if (!which) vtype[i] = Dump::INT;
else vtype[i] = Dump::DOUBLE;
int n = strlen(arg[iarg]); int n = strlen(arg[iarg]);
char *suffix = new char[n]; char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]); strcpy(suffix,&arg[iarg][2]);
argindex[i] = 0; argindex[i] = 0;
int tmp = -1; int flag,cols;
n = atom->find_custom(suffix,tmp); n = atom->find_custom(suffix,flag,cols);
if (n < 0) if ((!which && (n < 0 || flag || cols)) ||
error->all(FLERR,"Could not find custom per-atom property ID"); (which && (n < 0 || !flag || cols)))
error->all(FLERR,"Dump custom per-atom custom vector does not exist");
if (tmp != 1) field2index[i] = add_custom(suffix);
error->all(FLERR,"Custom per-atom property ID is not floating point");
field2index[i] = add_custom(suffix,1);
delete [] suffix; delete [] suffix;
// custom per-atom integer value = i_ID // custom per-atom array = i2_ID or d2_ID, must include bracketed index
} else if (strncmp(arg[iarg],"i2_",3) == 0 ||
strncmp(arg[iarg],"d2_",3) == 0) {
int which = 0;
if (arg[iarg][0] == 'd') which = 1;
} else if (strncmp(arg[iarg],"i_",2) == 0) {
pack_choice[i] = &DumpCustom::pack_custom; pack_choice[i] = &DumpCustom::pack_custom;
vtype[i] = Dump::INT; if (!which) vtype[i] = Dump::INT;
else vtype[i] = Dump::DOUBLE;
int n = strlen(arg[iarg]); int n = strlen(arg[iarg]);
char *suffix = new char[n]; char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]); strcpy(suffix,&arg[iarg][3]);
argindex[i] = 0;
int tmp = -1; char *ptr = strchr(suffix,'[');
n = atom->find_custom(suffix,tmp); if (ptr) {
if (n < 0) if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Could not find custom per-atom property ID"); error->all(FLERR,"Invalid attribute in dump custom command");
suffix[strlen(suffix)-1] = '\0';
argindex[i] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0';
} else error->all(FLERR,"Dump custom per-atom custom array is not indexed");
if (tmp != 0) int flag,cols;
error->all(FLERR,"Custom per-atom property ID is not integer"); n = atom->find_custom(suffix,flag,cols);
if ((!which && (n < 0 || flag || !cols)) ||
(which && (n < 0 || !flag || !cols)))
error->all(FLERR,"Dump custom per-atom custom array does not exist");
if (argindex[i] <= 0 || argindex[i] > cols)
error->all(FLERR,
"Dump custom per-atom custom array is accessed out-of-range");
field2index[i] = add_custom(suffix,0); field2index[i] = add_custom(suffix);
delete [] suffix; delete [] suffix;
} else return iarg; } else return iarg;
@ -1645,7 +1688,7 @@ int DumpCustom::add_variable(char *id)
id_variable = (char **) id_variable = (char **)
memory->srealloc(id_variable,(nvariable+1)*sizeof(char *), memory->srealloc(id_variable,(nvariable+1)*sizeof(char *),
"dump:id_variable"); "d<ump:id_variable");
delete [] variable; delete [] variable;
variable = new int[nvariable+1]; variable = new int[nvariable+1];
delete [] vbuf; delete [] vbuf;
@ -1661,27 +1704,28 @@ int DumpCustom::add_variable(char *id)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add custom atom property to list used by dump add custom atom property to list used by dump
return index of where this property is in list return index of where this property is in Atom class custom lists
if already in list, do not add, just return index, else add to list if already in list, do not add, just return index, else add to list
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int DumpCustom::add_custom(char *id, int flag) int DumpCustom::add_custom(char *id)
{ {
int icustom; int icustom;
for (icustom = 0; icustom < ncustom; icustom++) for (icustom = 0; icustom < ncustom; icustom++)
if ((strcmp(id,id_custom[icustom]) == 0) if (strcmp(id,id_custom[icustom]) == 0) break;
&& (flag == flag_custom[icustom])) break;
if (icustom < ncustom) return icustom; if (icustom < ncustom) return icustom;
id_custom = (char **) id_custom = (char **)
memory->srealloc(id_custom,(ncustom+1)*sizeof(char *),"dump:id_custom"); memory->srealloc(id_custom,(ncustom+1)*sizeof(char *),"dump:id_custom");
flag_custom = (int *)
memory->srealloc(flag_custom,(ncustom+1)*sizeof(int),"dump:flag_custom"); delete [] custom;
custom = new int[ncustom+1];
delete [] custom_flag;
custom_flag = new int[ncustom+1];
int n = strlen(id) + 1; int n = strlen(id) + 1;
id_custom[ncustom] = new char[n]; id_custom[ncustom] = new char[n];
strcpy(id_custom[ncustom],id); strcpy(id_custom[ncustom],id);
flag_custom[ncustom] = flag;
ncustom++; ncustom++;
return ncustom-1; return ncustom-1;
@ -1777,7 +1821,7 @@ int DumpCustom::modify_param(int narg, char **arg)
if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
if (strncmp(arg[1],"c_",2) != 0) if (strncmp(arg[1],"c_",2) != 0)
error->all(FLERR,"Illegal dump_modify command"); error->all(FLERR,"Illegal dump_modify command");
if (refreshflag) error->all(FLERR,"Dump modify can only have one refresh"); if (refreshflag) error->all(FLERR,"Dump_modify can only have one refresh");
refreshflag = 1; refreshflag = 1;
int n = strlen(arg[1]); int n = strlen(arg[1]);
@ -1914,28 +1958,29 @@ int DumpCustom::modify_param(int narg, char **arg)
char *ptr = strchr(suffix,'['); char *ptr = strchr(suffix,'[');
if (ptr) { if (ptr) {
if (suffix[strlen(suffix)-1] != ']') if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in dump modify command"); error->all(FLERR,"Invalid attribute in dump_modify command");
argindex[nfield+nthresh] = atoi(ptr+1); suffix[strlen(suffix)-1] = '\0';
argindex[nfield+nthresh] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0'; *ptr = '\0';
} else argindex[nfield+nthresh] = 0; } else argindex[nfield+nthresh] = 0;
n = modify->find_compute(suffix); n = modify->find_compute(suffix);
if (n < 0) error->all(FLERR,"Could not find dump modify compute ID"); if (n < 0) error->all(FLERR,"Could not find dump_modify compute ID");
if (modify->compute[n]->peratom_flag == 0) if (modify->compute[n]->peratom_flag == 0)
error->all(FLERR, error->all(FLERR,
"Dump modify compute ID does not compute per-atom info"); "Dump_modify compute ID does not compute per-atom info");
if (argindex[nfield+nthresh] == 0 && if (argindex[nfield+nthresh] == 0 &&
modify->compute[n]->size_peratom_cols > 0) modify->compute[n]->size_peratom_cols > 0)
error->all(FLERR, error->all(FLERR,
"Dump modify compute ID does not compute per-atom vector"); "Dump_modify compute ID does not compute per-atom vector");
if (argindex[nfield+nthresh] > 0 && if (argindex[nfield+nthresh] > 0 &&
modify->compute[n]->size_peratom_cols == 0) modify->compute[n]->size_peratom_cols == 0)
error->all(FLERR, error->all(FLERR,
"Dump modify compute ID does not compute per-atom array"); "Dump_modify compute ID does not compute per-atom array");
if (argindex[nfield+nthresh] > 0 && if (argindex[nfield+nthresh] > 0 &&
argindex[nfield+nthresh] > modify->compute[n]->size_peratom_cols) argindex[nfield+nthresh] > modify->compute[n]->size_peratom_cols)
error->all(FLERR,"Dump modify compute ID vector is not large enough"); error->all(FLERR,"Dump_modify compute ID vector is not large enough");
field2index[nfield+nthresh] = add_compute(suffix); field2index[nfield+nthresh] = add_compute(suffix);
delete [] suffix; delete [] suffix;
@ -1955,25 +2000,26 @@ int DumpCustom::modify_param(int narg, char **arg)
char *ptr = strchr(suffix,'['); char *ptr = strchr(suffix,'[');
if (ptr) { if (ptr) {
if (suffix[strlen(suffix)-1] != ']') if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in dump modify command"); error->all(FLERR,"Invalid attribute in dump_modify command");
argindex[nfield+nthresh] = atoi(ptr+1); suffix[strlen(suffix)-1] = '\0';
argindex[nfield+nthresh] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0'; *ptr = '\0';
} else argindex[nfield+nthresh] = 0; } else argindex[nfield+nthresh] = 0;
n = modify->find_fix(suffix); n = modify->find_fix(suffix);
if (n < 0) error->all(FLERR,"Could not find dump modify fix ID"); if (n < 0) error->all(FLERR,"Could not find dump_modify fix ID");
if (modify->fix[n]->peratom_flag == 0) if (modify->fix[n]->peratom_flag == 0)
error->all(FLERR,"Dump modify fix ID does not compute per-atom info"); error->all(FLERR,"Dump_modify fix ID does not compute per-atom info");
if (argindex[nfield+nthresh] == 0 && if (argindex[nfield+nthresh] == 0 &&
modify->fix[n]->size_peratom_cols > 0) modify->fix[n]->size_peratom_cols > 0)
error->all(FLERR,"Dump modify fix ID does not compute per-atom vector"); error->all(FLERR,"Dump_modify fix ID does not compute per-atom vector");
if (argindex[nfield+nthresh] > 0 && if (argindex[nfield+nthresh] > 0 &&
modify->fix[n]->size_peratom_cols == 0) modify->fix[n]->size_peratom_cols == 0)
error->all(FLERR,"Dump modify fix ID does not compute per-atom array"); error->all(FLERR,"Dump_modify fix ID does not compute per-atom array");
if (argindex[nfield+nthresh] > 0 && if (argindex[nfield+nthresh] > 0 &&
argindex[nfield+nthresh] > modify->fix[n]->size_peratom_cols) argindex[nfield+nthresh] > modify->fix[n]->size_peratom_cols)
error->all(FLERR,"Dump modify fix ID vector is not large enough"); error->all(FLERR,"Dump_modify fix ID vector is not large enough");
field2index[nfield+nthresh] = add_fix(suffix); field2index[nfield+nthresh] = add_fix(suffix);
delete [] suffix; delete [] suffix;
@ -1992,18 +2038,23 @@ int DumpCustom::modify_param(int narg, char **arg)
argindex[nfield+nthresh] = 0; argindex[nfield+nthresh] = 0;
n = input->variable->find(suffix); n = input->variable->find(suffix);
if (n < 0) error->all(FLERR,"Could not find dump modify variable name"); if (n < 0) error->all(FLERR,"Could not find dump_modify variable name");
if (input->variable->atomstyle(n) == 0) if (input->variable->atomstyle(n) == 0)
error->all(FLERR,"Dump modify variable is not atom-style variable"); error->all(FLERR,"Dump_modify variable is not atom-style variable");
field2index[nfield+nthresh] = add_variable(suffix); field2index[nfield+nthresh] = add_variable(suffix);
delete [] suffix; delete [] suffix;
// custom per atom floating point value = d_ID // custom per-atom vector = i_ID or d_ID
// must grow field2index and argindex arrays, since access is beyond nfield // must grow field2index and argindex arrays, since access is beyond nfield
} else if (strncmp(arg[1],"d_",2) == 0) { } else if (strncmp(arg[1],"i_",2) == 0 ||
thresh_array[nthresh] = DNAME; strncmp(arg[1],"d_",2) == 0) {
int which = 0;
if (arg[1][0] == 'd') which = 1;
if (!which) thresh_array[nthresh] = IVEC;
else thresh_array[nthresh] = DVEC;
memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); memory->grow(field2index,nfield+nthresh+1,"dump:field2index");
memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); memory->grow(argindex,nfield+nthresh+1,"dump:argindex");
int n = strlen(arg[1]); int n = strlen(arg[1]);
@ -2011,34 +2062,51 @@ int DumpCustom::modify_param(int narg, char **arg)
strcpy(suffix,&arg[1][2]); strcpy(suffix,&arg[1][2]);
argindex[nfield+nthresh] = 0; argindex[nfield+nthresh] = 0;
int tmp = -1; int flag,cols;
n = atom->find_custom(suffix,tmp); n = atom->find_custom(suffix,flag,cols);
if ((n < 0) || (tmp != 1)) if ((!which && (n < 0 || flag || cols)) ||
error->all(FLERR,"Could not find dump modify " (which && (n < 0 || !flag || cols)))
"custom atom floating point property ID"); error->all(FLERR,"Could not find dump_modify per-atom custom vector");
field2index[nfield+nthresh] = add_custom(suffix);
field2index[nfield+nthresh] = add_custom(suffix,1);
delete [] suffix; delete [] suffix;
// custom per atom integer value = i_ID // custom per-atom array = i2_ID or d2_ID, must include bracketed index
// must grow field2index and argindex arrays, since access is beyond nfield // must grow field2index and argindex arrays, since access is beyond nfield
} else if (strncmp(arg[1],"i_",2) == 0) { } else if (strncmp(arg[1],"i2_",3) == 0 ||
thresh_array[nthresh] = INAME; strncmp(arg[1],"d2_",3) == 0) {
int which = 0;
if (arg[1][0] == 'd') which = 1;
if (!which) thresh_array[nthresh] = IARRAY;
else thresh_array[nthresh] = DARRAY;
memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); memory->grow(field2index,nfield+nthresh+1,"dump:field2index");
memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); memory->grow(argindex,nfield+nthresh+1,"dump:argindex");
int n = strlen(arg[1]); int n = strlen(arg[1]);
char *suffix = new char[n]; char *suffix = new char[n];
strcpy(suffix,&arg[1][2]); strcpy(suffix,&arg[1][3]);
argindex[nfield+nthresh] = 0;
int tmp = -1; char *ptr = strchr(suffix,'[');
n = atom->find_custom(suffix,tmp); if (ptr) {
if ((n < 0) || (tmp != 0)) if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Could not find dump modify " error->all(FLERR,"Invalid attribute in dump custom command");
"custom atom integer property ID"); suffix[strlen(suffix)-1] = '\0';
argindex[nfield+nthresh] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0';
} else error->all(FLERR,"Dump_modify per-atom custom array is not indexed");
field2index[nfield+nthresh] = add_custom(suffix,0); int flag,cols;
n = atom->find_custom(suffix,flag,cols);
if ((!which && (n < 0 || flag || !cols)) ||
(which && (n < 0 || !flag || !cols)))
error->all(FLERR,"Could not find dump_modify per-atom custom array");
if (argindex[nfield+nthresh] <= 0 || argindex[nfield+nthresh] > cols)
error->all(FLERR,
"Dump_modify per-atom custom array is accessed out-of-range");
field2index[nfield+nthresh] = add_custom(suffix);
delete [] suffix; delete [] suffix;
} else error->all(FLERR,"Invalid dump_modify thresh attribute"); } else error->all(FLERR,"Invalid dump_modify thresh attribute");
@ -2165,27 +2233,36 @@ void DumpCustom::pack_variable(int n)
void DumpCustom::pack_custom(int n) void DumpCustom::pack_custom(int n)
{ {
int flag = custom_flag[field2index[n]];
int index = field2index[n]; int iwhich = custom[field2index[n]];
int index = argindex[n];
if (flag_custom[index] == 0) { // integer
int iwhich,tmp; if (flag == IVEC) {
iwhich = atom->find_custom(id_custom[index],tmp);
int *ivector = atom->ivector[iwhich]; int *ivector = atom->ivector[iwhich];
for (int i = 0; i < nchoose; i++) { for (int i = 0; i < nchoose; i++) {
buf[n] = ivector[clist[i]]; buf[n] = ivector[clist[i]];
n += size_one; n += size_one;
} }
} else if (flag_custom[index] == 1) { // double } else if (flag == DVEC) {
int iwhich,tmp;
iwhich = atom->find_custom(id_custom[index],tmp);
double *dvector = atom->dvector[iwhich]; double *dvector = atom->dvector[iwhich];
for (int i = 0; i < nchoose; i++) { for (int i = 0; i < nchoose; i++) {
buf[n] = dvector[clist[i]]; buf[n] = dvector[clist[i]];
n += size_one; n += size_one;
} }
} else if (flag == IARRAY) {
index--;
int **iarray = atom->iarray[iwhich];
for (int i = 0; i < nchoose; i++) {
buf[n] = iarray[clist[i]][index];
n += size_one;
}
} else if (flag == DARRAY) {
index--;
double **darray = atom->darray[iwhich];
for (int i = 0; i < nchoose; i++) {
buf[n] = darray[clist[i]][index];
n += size_one;
}
} }
} }

View File

@ -70,29 +70,30 @@ class DumpCustom : public Dump {
int nfield; // # of keywords listed by user int nfield; // # of keywords listed by user
int ioptional; // index of start of optional args int ioptional; // index of start of optional args
int *field2index; // which compute,fix,variable calcs this field int *field2index; // which compute,fix,variable,custom calcs this field
int *argindex; // index into compute,fix scalar_atom,vector_atom int *argindex; // index into compute,fix,custom per-atom data
// 0 for scalar_atom, 1-N for vector_atom values // 0 for per-atom vector, 1-N for cols of per-atom array
int ncompute; // # of Compute objects used by dump int ncompute; // # of Computes accessed by dump
char **id_compute; // their IDs char **id_compute; // their IDs
class Compute **compute; // list of ptrs to the Compute objects class Compute **compute; // list of ptrs to the Computes
int nfix; // # of Fix objects used by dump int nfix; // # of Fixes used by dump
char **id_fix; // their IDs char **id_fix; // their IDs
class Fix **fix; // list of ptrs to the Fix objects class Fix **fix; // list of ptrs to the Fixes
int nvariable; // # of Variables used by dump int nvariable; // # of Variables used by dump
char **id_variable; // their names char **id_variable; // their names
int *variable; // list of indices for the Variables int *variable; // list of Variable indices in Variable class
double **vbuf; // local storage for variable evaluation double **vbuf; // local storage for variable evaluation
int ncustom; // # of custom atom properties int ncustom; // # of Custom atom properties used by dump
char **id_custom; // their names char **id_custom; // their names
int *flag_custom; // their data type int *custom; // list of Custom indices in Atom class
int *custom_flag; // list of IVEC,DVEC,IARRAY,DARRAY styles
int ntypes; // # of atom types int ntypes; // # of atom types
char **typenames; // array of element names for each type char **typenames; // array of element names for each type
// private methods // private methods
@ -108,7 +109,7 @@ class DumpCustom : public Dump {
int add_compute(char *); int add_compute(char *);
int add_fix(char *); int add_fix(char *);
int add_variable(char *); int add_variable(char *);
int add_custom(char *, int); int add_custom(char *);
virtual int modify_param(int, char **); virtual int modify_param(int, char **);
void header_format_binary(); void header_format_binary();

View File

@ -64,6 +64,7 @@ idregion(NULL), idvar(NULL), idprop(NULL)
idregion = new char[n]; idregion = new char[n];
strcpy(idregion,arg[iarg+1]); strcpy(idregion,arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"var") == 0) { } else if (strcmp(arg[iarg],"var") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal group command"); if (iarg+2 > narg) error->all(FLERR,"Illegal group command");
if (input->variable->find(arg[iarg+1]) < 0) if (input->variable->find(arg[iarg+1]) < 0)
@ -74,16 +75,21 @@ idregion(NULL), idvar(NULL), idprop(NULL)
idvar = new char[n]; idvar = new char[n];
strcpy(idvar,arg[iarg+1]); strcpy(idvar,arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"property") == 0) { } else if (strcmp(arg[iarg],"property") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal group command"); if (iarg+2 > narg) error->all(FLERR,"Illegal group command");
if (atom->find_custom(arg[iarg+1],typeflag) < 0) int flag,cols;
error->all(FLERR,"Per atom property for group dynamic does not exist"); iprop = atom->find_custom(arg[iarg+1],flag,cols);
if (iprop < 1 || cols)
error->all(FLERR,"Custom per-atom vector for group dynamic "
"does not exist");
propflag = 1; propflag = 1;
delete [] idprop; delete [] idprop;
int n = strlen(arg[iarg+1]) + 1; int n = strlen(arg[iarg+1]) + 1;
idprop = new char[n]; idprop = new char[n];
strcpy(idprop,arg[iarg+1]); strcpy(idprop,arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"every") == 0) { } else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal group command"); if (iarg+2 > narg) error->all(FLERR,"Illegal group command");
nevery = force->inumeric(FLERR,arg[iarg+1]); nevery = force->inumeric(FLERR,arg[iarg+1]);
@ -125,7 +131,7 @@ void FixGroup::init()
if (strstr(update->integrate_style,"respa")) if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels; nlevels_respa = ((Respa *) update->integrate)->nlevels;
// set current indices for region and variable // set current indices for region and variable and custom property
if (regionflag) { if (regionflag) {
iregion = domain->find_region(idregion); iregion = domain->find_region(idregion);
@ -143,9 +149,10 @@ void FixGroup::init()
} }
if (propflag) { if (propflag) {
iprop = atom->find_custom(idprop,typeflag); int cols;
if (iprop < 0) iprop = atom->find_custom(idprop,proptype,cols);
error->all(FLERR,"Per-atom property for group dynamic does not exist"); if (iprop < 0 || cols)
error->all(FLERR,"Group dynamic command custom property vector does not exist");
} }
// warn if any FixGroup is not at tail end of all post_integrate fixes // warn if any FixGroup is not at tail end of all post_integrate fixes
@ -222,11 +229,10 @@ void FixGroup::set_group()
update->post_integrate = 0; update->post_integrate = 0;
} }
// invoke per-atom property if defined // set ptr to custom atom vector
if (propflag && !typeflag) ivector = atom->ivector[iprop]; //check nlocal > 0? if (propflag && !proptype) ivector = atom->ivector[iprop];
if (propflag && proptype) dvector = atom->dvector[iprop];
if (propflag && typeflag) dvector = atom->dvector[iprop]; //check nlocal > 0?
// update region in case it has a variable dependence or is dynamic // update region in case it has a variable dependence or is dynamic
@ -246,8 +252,10 @@ void FixGroup::set_group()
inflag = 1; inflag = 1;
if (regionflag && !region->match(x[i][0],x[i][1],x[i][2])) inflag = 0; if (regionflag && !region->match(x[i][0],x[i][1],x[i][2])) inflag = 0;
if (varflag && var[i] == 0.0) inflag = 0; if (varflag && var[i] == 0.0) inflag = 0;
if (propflag && !typeflag && ivector[i] == 0) inflag = 0; if (propflag) {
if (propflag && typeflag && dvector[i] == 0) inflag = 0; if (!proptype && ivector[i] == 0) inflag = 0;
if (proptype && dvector[i] == 0.0) inflag = 0;
}
} else inflag = 0; } else inflag = 0;
if (inflag) mask[i] |= gbit; if (inflag) mask[i] |= gbit;

View File

@ -37,7 +37,7 @@ class FixGroup : public Fix {
private: private:
int gbit,gbitinverse; int gbit,gbitinverse;
int regionflag,varflag,propflag,typeflag; int regionflag,varflag,propflag,proptype;
int iregion,ivar,iprop; int iregion,ivar,iprop;
char *idregion,*idvar,*idprop; char *idregion,*idvar,*idprop;
class Region *region; class Region *region;

View File

@ -24,7 +24,7 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{MOLECULE,CHARGE,RMASS,INTEGER,DOUBLE}; enum{MOLECULE,CHARGE,RMASS,IVEC,DVEC,IARRAY,DARRAY};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -40,6 +40,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3; int iarg = 3;
nvalue = narg-iarg; nvalue = narg-iarg;
style = new int[nvalue]; style = new int[nvalue];
cols = new int[nvalue];
index = new int[nvalue]; index = new int[nvalue];
molecule_flag = 0; molecule_flag = 0;
@ -47,6 +48,8 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
rmass_flag = 0; rmass_flag = 0;
nvalue = 0; nvalue = 0;
values_peratom = 0;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"mol") == 0) { if (strcmp(arg[iarg],"mol") == 0) {
if (atom->molecule_flag) if (atom->molecule_flag)
@ -55,8 +58,11 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (molecule_flag) if (molecule_flag)
error->all(FLERR,"Fix property/atom cannot specify mol twice"); error->all(FLERR,"Fix property/atom cannot specify mol twice");
style[nvalue] = MOLECULE; style[nvalue] = MOLECULE;
cols[nvalue] = 0;
atom->molecule_flag = molecule_flag = 1; atom->molecule_flag = molecule_flag = 1;
values_peratom++;
nvalue++; nvalue++;
iarg++;
} else if (strcmp(arg[iarg],"q") == 0) { } else if (strcmp(arg[iarg],"q") == 0) {
if (atom->q_flag) if (atom->q_flag)
error->all(FLERR,"Fix property/atom q when atom_style " error->all(FLERR,"Fix property/atom q when atom_style "
@ -64,8 +70,11 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (q_flag) if (q_flag)
error->all(FLERR,"Fix property/atom cannot specify q twice"); error->all(FLERR,"Fix property/atom cannot specify q twice");
style[nvalue] = CHARGE; style[nvalue] = CHARGE;
cols[nvalue] = 0;
atom->q_flag = q_flag = 1; atom->q_flag = q_flag = 1;
values_peratom++;
nvalue++; nvalue++;
iarg++;
} else if (strcmp(arg[iarg],"rmass") == 0) { } else if (strcmp(arg[iarg],"rmass") == 0) {
if (atom->rmass_flag) if (atom->rmass_flag)
error->all(FLERR,"Fix property/atom rmass when atom_style " error->all(FLERR,"Fix property/atom rmass when atom_style "
@ -73,27 +82,45 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (rmass_flag) if (rmass_flag)
error->all(FLERR,"Fix property/atom cannot specify rmass twice"); error->all(FLERR,"Fix property/atom cannot specify rmass twice");
style[nvalue] = RMASS; style[nvalue] = RMASS;
cols[nvalue] = 0;
atom->rmass_flag = rmass_flag = 1; atom->rmass_flag = rmass_flag = 1;
values_peratom++;
nvalue++; nvalue++;
} else if (strstr(arg[iarg],"i_") == arg[iarg]) { iarg++;
style[nvalue] = INTEGER; } else if (strstr(arg[iarg],"i_") == arg[iarg] ||
int tmp; strstr(arg[iarg],"d_") == arg[iarg]) {
index[nvalue] = atom->find_custom(&arg[iarg][2],tmp); int which = 0;
if (arg[iarg][0] == 'd') which = 1;
if (which == 0) style[nvalue] = IVEC;
else style[nvalue] = DVEC;
int tmp1,tmp2;
index[nvalue] = atom->find_custom(&arg[iarg][2],tmp1,tmp2);
if (index[nvalue] >= 0) if (index[nvalue] >= 0)
error->all(FLERR,"Fix property/atom vector name already exists"); error->all(FLERR,"Fix property/atom vector name already exists");
index[nvalue] = atom->add_custom(&arg[iarg][2],0); cols[nvalue] = 0;
index[nvalue] = atom->add_custom(&arg[iarg][2],which,cols[nvalue]);
values_peratom++;
nvalue++; nvalue++;
} else if (strstr(arg[iarg],"d_") == arg[iarg]) { iarg++;
style[nvalue] = DOUBLE; } else if (strstr(arg[iarg],"i2_") == arg[iarg] ||
int tmp; strstr(arg[iarg],"d2_") == arg[iarg]) {
index[nvalue] = atom->find_custom(&arg[iarg][2],tmp); if (iarg+2 > narg) error->all(FLERR,"Illegal fix property/atom command");
int which = 0;
if (arg[iarg][0] == 'd') which = 1;
if (which == 0) style[nvalue] = IARRAY;
else style[nvalue] = DARRAY;
int tmp1,tmp2;
index[nvalue] = atom->find_custom(&arg[iarg][3],tmp1,tmp2);
if (index[nvalue] >= 0) if (index[nvalue] >= 0)
error->all(FLERR,"Fix property/atom vector name already exists"); error->all(FLERR,"Fix property/atom array name already exists");
index[nvalue] = atom->add_custom(&arg[iarg][2],1); cols[nvalue] = utils::inumeric(FLERR,arg[iarg+1],true,lmp);
if (cols[nvalue] < 1)
error->all(FLERR,"Invalid array columns in fix property/atom");
index[nvalue] = atom->add_custom(&arg[iarg][3],which,cols[nvalue]);
values_peratom += cols[nvalue];
nvalue++; nvalue++;
iarg += 2;
} else break; } else break;
iarg++;
} }
// optional args // optional args
@ -109,7 +136,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal fix property/atom command"); } else error->all(FLERR,"Illegal fix property/atom command");
} }
if (border) comm_border = nvalue; if (border) comm_border = values_peratom;
// warn if mol or charge keyword used without ghost yes // warn if mol or charge keyword used without ghost yes
@ -134,8 +161,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
// register with Atom class // register with Atom class
nmax_old = 0; nmax_old = 0;
if (!lmp->kokkos) if (!lmp->kokkos) grow_arrays(atom->nmax);
grow_arrays(atom->nmax);
atom->add_callback(0); atom->add_callback(0);
atom->add_callback(1); atom->add_callback(1);
if (border) atom->add_callback(2); if (border) atom->add_callback(2);
@ -154,27 +180,32 @@ FixPropertyAtom::~FixPropertyAtom()
// deallocate per-atom vectors in Atom class // deallocate per-atom vectors in Atom class
// set ptrs to NULL, so they no longer exist for Atom class // set ptrs to NULL, so they no longer exist for Atom class
for (int m = 0; m < nvalue; m++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[m] == MOLECULE) { if (style[nv] == MOLECULE) {
atom->molecule_flag = 0; atom->molecule_flag = 0;
memory->destroy(atom->molecule); memory->destroy(atom->molecule);
atom->molecule = NULL; atom->molecule = NULL;
} else if (style[m] == CHARGE) { } else if (style[nv] == CHARGE) {
atom->q_flag = 0; atom->q_flag = 0;
memory->destroy(atom->q); memory->destroy(atom->q);
atom->q = NULL; atom->q = NULL;
} else if (style[m] == RMASS) { } else if (style[nv] == RMASS) {
atom->rmass_flag = 0; atom->rmass_flag = 0;
memory->destroy(atom->rmass); memory->destroy(atom->rmass);
atom->rmass = NULL; atom->rmass = NULL;
} else if (style[m] == INTEGER) { } else if (style[nv] == IVEC) {
atom->remove_custom(0,index[m]); atom->remove_custom(index[nv],0,cols[nv]);
} else if (style[m] == DOUBLE) { } else if (style[nv] == DVEC) {
atom->remove_custom(1,index[m]); atom->remove_custom(index[nv],1,cols[nv]);
} else if (style[nv] == IARRAY) {
atom->remove_custom(index[nv],0,cols[nv]);
} else if (style[nv] == DARRAY) {
atom->remove_custom(index[nv],1,cols[nv]);
} }
} }
delete [] style; delete [] style;
delete [] cols;
delete [] index; delete [] index;
delete [] astyle; delete [] astyle;
} }
@ -206,7 +237,7 @@ void FixPropertyAtom::init()
void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf, void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf,
tagint id_offset) tagint id_offset)
{ {
int j,m; int j,k,m,iword,ncol,nv;
tagint itag; tagint itag;
char *next; char *next;
@ -222,7 +253,7 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf,
int nwords = utils::trim_and_count_words(buf); int nwords = utils::trim_and_count_words(buf);
*next = '\n'; *next = '\n';
if (nwords != nvalue+1) if (nwords != values_peratom+1)
error->all(FLERR,fmt::format("Incorrect {} format in data file",keyword)); error->all(FLERR,fmt::format("Incorrect {} format in data file",keyword));
char **values = new char*[nwords]; char **values = new char*[nwords];
@ -254,16 +285,27 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf,
"data file",itag, keyword)); "data file",itag, keyword));
// assign words in line to per-atom vectors // assign words in line to per-atom vectors
// iword = position in vector of words
if ((m = atom->map(itag)) >= 0) { if ((m = atom->map(itag)) >= 0) {
for (j = 0; j < nvalue; j++) { iword = 1;
if (style[j] == MOLECULE) atom->molecule[m] = ATOTAGINT(values[j+1]); for (nv = 0; nv < nvalue; nv++) {
else if (style[j] == CHARGE) atom->q[m] = atof(values[j+1]); if (style[nv] == MOLECULE) atom->molecule[m] = ATOTAGINT(values[iword++]);
else if (style[j] == RMASS) atom->rmass[m] = atof(values[j+1]); else if (style[nv] == CHARGE) atom->q[m] = atof(values[iword++]);
else if (style[j] == INTEGER) else if (style[nv] == RMASS) atom->rmass[m] = atof(values[iword++]);
atom->ivector[index[j]][m] = atoi(values[j+1]); else if (style[nv] == IVEC)
else if (style[j] == DOUBLE) atom->ivector[index[nv]][m] = atoi(values[iword++]);
atom->dvector[index[j]][m] = atof(values[j+1]); else if (style[nv] == DVEC)
atom->dvector[index[nv]][m] = atof(values[iword++]);
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][m][k] = atoi(values[iword++]);
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][m][k] = atof(values[iword++]);
}
} }
} }
@ -291,13 +333,13 @@ bigint FixPropertyAtom::read_data_skip_lines(char * /*keyword*/)
return size I own for Mth data section return size I own for Mth data section
# of data sections = 1 for this fix # of data sections = 1 for this fix
nx = # of local atoms nx = # of local atoms
ny = columns = tag + nvalues ny = columns = tag + values_peratom
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixPropertyAtom::write_data_section_size(int /*mth*/, int &nx, int &ny) void FixPropertyAtom::write_data_section_size(int /*mth*/, int &nx, int &ny)
{ {
nx = atom->nlocal; nx = atom->nlocal;
ny = nvalue + 1; ny = values_peratom + 1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -307,7 +349,7 @@ void FixPropertyAtom::write_data_section_size(int /*mth*/, int &nx, int &ny)
void FixPropertyAtom::write_data_section_pack(int /*mth*/, double **buf) void FixPropertyAtom::write_data_section_pack(int /*mth*/, double **buf)
{ {
int i; int i,k,ncol;
// 1st column = atom tag // 1st column = atom tag
// rest of columns = per-atom values // rest of columns = per-atom values
@ -317,23 +359,42 @@ void FixPropertyAtom::write_data_section_pack(int /*mth*/, double **buf)
for (i = 0; i < nlocal; i++) buf[i][0] = ubuf(tag[i]).d; for (i = 0; i < nlocal; i++) buf[i][0] = ubuf(tag[i]).d;
for (int m = 0; m < nvalue; m++) { int icol = 1;
int mp1 = m+1; for (int nv = 0; nv < nvalue; nv++) {
if (style[m] == MOLECULE) { if (style[nv] == MOLECULE) {
tagint *molecule = atom->molecule; tagint *molecule = atom->molecule;
for (i = 0; i < nlocal; i++) buf[i][mp1] = ubuf(molecule[i]).d; for (i = 0; i < nlocal; i++) buf[i][icol] = ubuf(molecule[i]).d;
} else if (style[m] == CHARGE) { icol++;
} else if (style[nv] == CHARGE) {
double *q = atom->q; double *q = atom->q;
for (i = 0; i < nlocal; i++) buf[i][mp1] = q[i]; for (i = 0; i < nlocal; i++) buf[i][icol] = q[i];
} else if (style[m] == RMASS) { icol++;
} else if (style[nv] == RMASS) {
double *rmass = atom->rmass; double *rmass = atom->rmass;
for (i = 0; i < nlocal; i++) buf[i][mp1] = rmass[i]; for (i = 0; i < nlocal; i++) buf[i][icol] = rmass[i];
} else if (style[m] == INTEGER) { icol++;
int *ivec = atom->ivector[index[m]]; } else if (style[nv] == IVEC) {
for (i = 0; i < nlocal; i++) buf[i][mp1] = ubuf(ivec[i]).d; int *ivec = atom->ivector[index[nv]];
} else if (style[m] == DOUBLE) { for (i = 0; i < nlocal; i++) buf[i][icol] = ubuf(ivec[i]).d;
double *dvec = atom->dvector[index[m]]; icol++;
for (i = 0; i < nlocal; i++) buf[i][mp1] = dvec[i]; } else if (style[nv] == DVEC) {
double *dvec = atom->dvector[index[nv]];
for (i = 0; i < nlocal; i++) buf[i][icol] = dvec[i];
icol++;
} else if (style[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
for (i = 0; i < nlocal; i++)
for (k = 0; k < ncol; k++)
buf[i][icol+k] = ubuf(iarray[i][k]).d;
icol += ncol;
} else if (style[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
for (i = 0; i < nlocal; i++)
for (k = 0; k < ncol; k++)
buf[i][icol+k] = ubuf(darray[i][k]).d;
icol += ncol;
} }
} }
} }
@ -361,16 +422,33 @@ void FixPropertyAtom::write_data_section_keyword(int /*mth*/, FILE *fp)
void FixPropertyAtom::write_data_section(int /*mth*/, FILE *fp, void FixPropertyAtom::write_data_section(int /*mth*/, FILE *fp,
int n, double **buf, int /*index*/) int n, double **buf, int /*index*/)
{ {
int m; int k,icol,ncol,nv;
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
fprintf(fp,TAGINT_FORMAT,(tagint) ubuf(buf[i][0]).i); fprintf(fp,TAGINT_FORMAT,(tagint) ubuf(buf[i][0]).i);
for (m = 0; m < nvalue; m++) { icol = 1;
if (style[m] == MOLECULE) for (nv = 0; nv < nvalue; nv++) {
fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[i][m+1]).i); if (style[nv] == MOLECULE)
else if (style[m] == INTEGER) fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[i][icol++]).i);
fprintf(fp," %d",(int) ubuf(buf[i][m+1]).i); else if (style[nv] == CHARGE)
else fprintf(fp," %g",buf[i][m+1]); fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == RMASS)
fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == IVEC)
fprintf(fp," %d",(int) ubuf(buf[i][icol++]).i);
else if (style[nv] == DVEC)
fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
fprintf(fp," %d",(int) ubuf(buf[i][icol+k]).i);
icol += ncol;
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
fprintf(fp," %g",buf[i][icol+k]);
icol += ncol;
}
} }
fprintf(fp,"\n"); fprintf(fp,"\n");
} }
@ -387,42 +465,52 @@ double FixPropertyAtom::memory_usage()
if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(tagint); if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(tagint);
else if (style[m] == CHARGE) bytes = atom->nmax * sizeof(double); else if (style[m] == CHARGE) bytes = atom->nmax * sizeof(double);
else if (style[m] == RMASS) bytes = atom->nmax * sizeof(double); else if (style[m] == RMASS) bytes = atom->nmax * sizeof(double);
else if (style[m] == INTEGER) bytes = atom->nmax * sizeof(int); else if (style[m] == IVEC) bytes = atom->nmax * sizeof(int);
else if (style[m] == DOUBLE) bytes = atom->nmax * sizeof(double); else if (style[m] == DVEC) bytes = atom->nmax * sizeof(double);
else if (style[m] == IARRAY) bytes = atom->nmax * cols[m] * sizeof(int);
else if (style[m] == DARRAY) bytes = atom->nmax * cols[m] * sizeof(double);
} }
return bytes; return bytes;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
allocate atom-based arrays allocate atom-based arrays
initialize new values to 0, also initialize new values to 0
since AtomVec class won't do it as atoms are added, since AtomVec class won't do it as atoms are added,
e.g. in create_atom() or data_atom() e.g. in create_atom() or data_atom()
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixPropertyAtom::grow_arrays(int nmax) void FixPropertyAtom::grow_arrays(int nmax)
{ {
for (int m = 0; m < nvalue; m++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[m] == MOLECULE) { if (style[nv] == MOLECULE) {
memory->grow(atom->molecule,nmax,"atom:molecule"); memory->grow(atom->molecule,nmax,"atom:molecule");
size_t nbytes = (nmax-nmax_old) * sizeof(tagint); size_t nbytes = (nmax-nmax_old) * sizeof(tagint);
memset(&atom->molecule[nmax_old],0,nbytes); memset(&atom->molecule[nmax_old],0,nbytes);
} else if (style[m] == CHARGE) { } else if (style[nv] == CHARGE) {
memory->grow(atom->q,nmax,"atom:q"); memory->grow(atom->q,nmax,"atom:q");
size_t nbytes = (nmax-nmax_old) * sizeof(double); size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->q[nmax_old],0,nbytes); memset(&atom->q[nmax_old],0,nbytes);
} else if (style[m] == RMASS) { } else if (style[nv] == RMASS) {
memory->grow(atom->rmass,nmax,"atom:rmass"); memory->grow(atom->rmass,nmax,"atom:rmass");
size_t nbytes = (nmax-nmax_old) * sizeof(double); size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->rmass[nmax_old],0,nbytes); memset(&atom->rmass[nmax_old],0,nbytes);
} else if (style[m] == INTEGER) { } else if (style[nv] == IVEC) {
memory->grow(atom->ivector[index[m]],nmax,"atom:ivector"); memory->grow(atom->ivector[index[nv]],nmax,"atom:ivector");
size_t nbytes = (nmax-nmax_old) * sizeof(int); size_t nbytes = (nmax-nmax_old) * sizeof(int);
memset(&atom->ivector[index[m]][nmax_old],0,nbytes); memset(&atom->ivector[index[nv]][nmax_old],0,nbytes);
} else if (style[m] == DOUBLE) { } else if (style[nv] == DVEC) {
memory->grow(atom->dvector[index[m]],nmax,"atom:dvector"); memory->grow(atom->dvector[index[nv]],nmax,"atom:dvector");
size_t nbytes = (nmax-nmax_old) * sizeof(double); size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->dvector[index[m]][nmax_old],0,nbytes); memset(&atom->dvector[index[nv]][nmax_old],0,nbytes);
} else if (style[nv] == IARRAY) {
memory->grow(atom->iarray[index[nv]],nmax,cols[nv],"atom:iarray");
size_t nbytes = (nmax-nmax_old) * cols[nv] * sizeof(int);
if (nbytes) memset(&atom->iarray[index[nv]][nmax_old][0],0,nbytes);
} else if (style[nv] == DARRAY) {
memory->grow(atom->darray[index[nv]],nmax,cols[nv],"atom:darray");
size_t nbytes = (nmax-nmax_old) * cols[nv] * sizeof(double);
if (nbytes) memset(&atom->darray[index[nv]][nmax_old][0],0,nbytes);
} }
} }
@ -435,17 +523,28 @@ void FixPropertyAtom::grow_arrays(int nmax)
void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/) void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
{ {
for (int m = 0; m < nvalue; m++) { int k,ncol;
if (style[m] == MOLECULE)
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE)
atom->molecule[j] = atom->molecule[i]; atom->molecule[j] = atom->molecule[i];
else if (style[m] == CHARGE) else if (style[nv] == CHARGE)
atom->q[j] = atom->q[i]; atom->q[j] = atom->q[i];
else if (style[m] == RMASS) else if (style[nv] == RMASS)
atom->rmass[j] = atom->rmass[i]; atom->rmass[j] = atom->rmass[i];
else if (style[m] == INTEGER) else if (style[nv] == IVEC)
atom->ivector[index[m]][j] = atom->ivector[index[m]][i]; atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i];
else if (style[m] == DOUBLE) else if (style[nv] == DVEC)
atom->dvector[index[m]][j] = atom->dvector[index[m]][i]; atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
}
} }
} }
@ -455,40 +554,56 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
int FixPropertyAtom::pack_border(int n, int *list, double *buf) int FixPropertyAtom::pack_border(int n, int *list, double *buf)
{ {
int i,j,k; int i,j,k,ncol;
int m = 0; int m = 0;
for (k = 0; k < nvalue; k++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[k] == MOLECULE) { if (style[nv] == MOLECULE) {
tagint *molecule = atom->molecule; tagint *molecule = atom->molecule;
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
} }
} else if (style[k] == CHARGE) { } else if (style[nv] == CHARGE) {
double *q = atom->q; double *q = atom->q;
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = q[j]; buf[m++] = q[j];
} }
} else if (style[k] == RMASS) { } else if (style[nv] == RMASS) {
double *rmass = atom->rmass; double *rmass = atom->rmass;
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = rmass[j]; buf[m++] = rmass[j];
} }
} else if (style[k] == INTEGER) { } else if (style[nv] == IVEC) {
int *ivector = atom->ivector[index[k]]; int *ivector = atom->ivector[index[nv]];
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = ubuf(ivector[j]).d; buf[m++] = ubuf(ivector[j]).d;
} }
} else if (style[k] == DOUBLE) { } else if (style[nv] == DVEC) {
double *dvector = atom->dvector[index[k]]; double *dvector = atom->dvector[index[nv]];
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = dvector[j]; buf[m++] = dvector[j];
} }
} else if (style[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
for (i = 0; i < n; i++) {
j = list[i];
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(iarray[j][k]).d;
}
} else if (style[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
for (i = 0; i < n; i++) {
j = list[i];
for (k = 0; k < ncol; k++)
buf[m++] = darray[j][k];
}
} }
} }
@ -501,35 +616,49 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf)
int FixPropertyAtom::unpack_border(int n, int first, double *buf) int FixPropertyAtom::unpack_border(int n, int first, double *buf)
{ {
int i,k,last; int i,k,last,ncol;
int m = 0; int m = 0;
for (k = 0; k < nvalue; k++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[k] == MOLECULE) { if (style[nv] == MOLECULE) {
tagint *molecule = atom->molecule; tagint *molecule = atom->molecule;
last = first + n; last = first + n;
for (i = first; i < last; i++) for (i = first; i < last; i++)
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
} else if (style[k] == CHARGE) { } else if (style[nv] == CHARGE) {
double *q = atom->q; double *q = atom->q;
last = first + n; last = first + n;
for (i = first; i < last; i++) for (i = first; i < last; i++)
q[i] = buf[m++]; q[i] = buf[m++];
} else if (style[k] == RMASS) { } else if (style[nv] == RMASS) {
double *rmass = atom->rmass; double *rmass = atom->rmass;
last = first + n; last = first + n;
for (i = first; i < last; i++) for (i = first; i < last; i++)
rmass[i] = buf[m++]; rmass[i] = buf[m++];
} else if (style[k] == INTEGER) { } else if (style[nv] == IVEC) {
int *ivector = atom->ivector[index[k]]; int *ivector = atom->ivector[index[nv]];
last = first + n; last = first + n;
for (i = first; i < last; i++) for (i = first; i < last; i++)
ivector[i] = (int) ubuf(buf[m++]).i; ivector[i] = (int) ubuf(buf[m++]).i;
} else if (style[k] == DOUBLE) { } else if (style[nv] == DVEC) {
double *dvector = atom->dvector[index[k]]; double *dvector = atom->dvector[index[nv]];
last = first + n; last = first + n;
for (i = first; i < last; i++) for (i = first; i < last; i++)
dvector[i] = buf[m++]; dvector[i] = buf[m++];
} else if (style[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
last = first + n;
for (i = first; i < last; i++)
for (k = 0; k < ncol; k++)
iarray[i][k] = (int) ubuf(buf[m++]).i;
} else if (style[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
last = first + n;
for (i = first; i < last; i++)
for (k = 0; k < ncol; k++)
darray[i][k] = buf[m++];
} }
} }
@ -542,14 +671,27 @@ int FixPropertyAtom::unpack_border(int n, int first, double *buf)
int FixPropertyAtom::pack_exchange(int i, double *buf) int FixPropertyAtom::pack_exchange(int i, double *buf)
{ {
for (int m = 0; m < nvalue; m++) { int k,ncol;
if (style[m] == MOLECULE) buf[m] = ubuf(atom->molecule[i]).d;
else if (style[m] == CHARGE) buf[m] = atom->q[i]; int m = 0;
else if (style[m] == RMASS) buf[m] = atom->rmass[i]; for (int nv = 0; nv < nvalue; nv++) {
else if (style[m] == INTEGER) buf[m] = ubuf(atom->ivector[index[m]][i]).d; if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (style[m] == DOUBLE) buf[m] = atom->dvector[index[m]][i]; else if (style[nv] == CHARGE) buf[m++] = atom->q[i];
else if (style[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = atom->darray[index[nv]][i][k];
}
} }
return nvalue;
return m;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -558,19 +700,32 @@ int FixPropertyAtom::pack_exchange(int i, double *buf)
int FixPropertyAtom::unpack_exchange(int nlocal, double *buf) int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
{ {
for (int m = 0; m < nvalue; m++) { int k,ncol;
if (style[m] == MOLECULE)
atom->molecule[nlocal] = (tagint) ubuf(buf[m]).i; int m = 0;
else if (style[m] == CHARGE) for (int nv = 0; nv < nvalue; nv++) {
atom->q[nlocal] = buf[m]; if (style[nv] == MOLECULE)
else if (style[m] == RMASS) atom->molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
atom->rmass[nlocal] = buf[m]; else if (style[nv] == CHARGE)
else if (style[m] == INTEGER) atom->q[nlocal] = buf[m++];
atom->ivector[index[m]][nlocal] = (int) ubuf(buf[m]).i; else if (style[nv] == RMASS)
else if (style[m] == DOUBLE) atom->rmass[nlocal] = buf[m++];
atom->dvector[index[m]][nlocal] = buf[m]; else if (style[nv] == IVEC)
atom->ivector[index[nv]][nlocal] = (int) ubuf(buf[m++]).i;
else if (style[nv] == DVEC)
atom->dvector[index[nv]][nlocal] = buf[m++];
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][nlocal][k] = (int) ubuf(buf[m++]).i;
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][nlocal][k] = buf[m++];
}
} }
return nvalue;
return m;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -579,19 +734,31 @@ int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
int FixPropertyAtom::pack_restart(int i, double *buf) int FixPropertyAtom::pack_restart(int i, double *buf)
{ {
int k,ncol;
// pack buf[0] this way because other fixes unpack it // pack buf[0] this way because other fixes unpack it
buf[0] = nvalue+1;
buf[0] = values_peratom+1;
int m = 1; int m = 1;
for (int j = 0; j < nvalue; j++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[j] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d; if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (style[j] == CHARGE) buf[m++] = atom->q[i]; else if (style[nv] == CHARGE) buf[m++] = atom->q[i];
else if (style[j] == RMASS) buf[m++] = atom->rmass[i]; else if (style[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (style[j] == INTEGER) buf[m++] = ubuf(atom->ivector[index[j]][i]).d; else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (style[j] == DOUBLE) buf[m++] = atom->dvector[index[j]][i]; else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = atom->darray[index[nv]][i][k];
}
} }
return nvalue+1; return values_peratom+1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -600,6 +767,7 @@ int FixPropertyAtom::pack_restart(int i, double *buf)
void FixPropertyAtom::unpack_restart(int nlocal, int nth) void FixPropertyAtom::unpack_restart(int nlocal, int nth)
{ {
int k,ncol;
double **extra = atom->extra; double **extra = atom->extra;
// skip to Nth set of extra values // skip to Nth set of extra values
@ -609,17 +777,26 @@ void FixPropertyAtom::unpack_restart(int nlocal, int nth)
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]); for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++; m++;
for (int i = 0; i < nvalue; i++) { for (int nv = 0; nv < nvalue; nv++) {
if (style[i] == MOLECULE) if (style[nv] == MOLECULE)
atom->molecule[nlocal] = (tagint) ubuf(extra[nlocal][m++]).i; atom->molecule[nlocal] = (tagint) ubuf(extra[nlocal][m++]).i;
else if (style[i] == CHARGE) else if (style[nv] == CHARGE)
atom->q[nlocal] = extra[nlocal][m++]; atom->q[nlocal] = extra[nlocal][m++];
else if (style[i] == RMASS) else if (style[nv] == RMASS)
atom->rmass[nlocal] = extra[nlocal][m++]; atom->rmass[nlocal] = extra[nlocal][m++];
else if (style[i] == INTEGER) else if (style[nv] == IVEC)
atom->ivector[index[i]][nlocal] = (int) ubuf(extra[nlocal][m++]).i; atom->ivector[index[nv]][nlocal] = (int) ubuf(extra[nlocal][m++]).i;
else if (style[i] == DOUBLE) else if (style[nv] == DVEC)
atom->dvector[index[i]][nlocal] = extra[nlocal][m++]; atom->dvector[index[nv]][nlocal] = extra[nlocal][m++];
else if (style[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][nlocal][k] = (int) ubuf(extra[nlocal][m++]).i;
} else if (style[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][nlocal][k] = extra[nlocal][m++];
}
} }
} }
@ -629,7 +806,7 @@ void FixPropertyAtom::unpack_restart(int nlocal, int nth)
int FixPropertyAtom::maxsize_restart() int FixPropertyAtom::maxsize_restart()
{ {
return nvalue+1; return values_peratom+1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -638,5 +815,5 @@ int FixPropertyAtom::maxsize_restart()
int FixPropertyAtom::size_restart(int /*nlocal*/) int FixPropertyAtom::size_restart(int /*nlocal*/)
{ {
return nvalue+1; return values_peratom+1;
} }

View File

@ -52,10 +52,13 @@ class FixPropertyAtom : public Fix {
protected: protected:
int nvalue,border; int nvalue,border;
int molecule_flag,q_flag,rmass_flag; int molecule_flag,q_flag,rmass_flag; // flags for specific fields
int *style,*index; int *style; // style of each value, see enum
char *astyle; int *index; // indices into atom custom data structs
int *cols; // columns per value, for arrays
char *astyle; // atom style at instantiation
int values_peratom; // # of values per atom, including multiple for arrays
int nmax_old; // length of peratom arrays the last time they grew int nmax_old; // length of peratom arrays the last time they grew
}; };

View File

@ -20,17 +20,18 @@
#include "group.h" #include "group.h"
#include "modify.h" #include "modify.h"
#include "compute.h" #include "compute.h"
#include "force.h"
#include "fix.h" #include "fix.h"
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "utils.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "force.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{KEYWORD,COMPUTE,FIX,VARIABLE,DNAME,INAME}; enum{KEYWORD,COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY};
#define INVOKED_PERATOM 8 #define INVOKED_PERATOM 8
@ -224,15 +225,11 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
pack_choice[nvalues++] = &FixStoreState::pack_tqz; pack_choice[nvalues++] = &FixStoreState::pack_tqz;
} else if (strncmp(arg[iarg],"c_",2) == 0 || } else if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"d_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 || strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"i_",2) == 0 ||
strncmp(arg[iarg],"v_",2) == 0) { strncmp(arg[iarg],"v_",2) == 0) {
cfv_any = 1; cfv_any = 1;
if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE; if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
else if (arg[iarg][0] == 'd') which[nvalues] = DNAME;
else if (arg[iarg][0] == 'f') which[nvalues] = FIX; else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
else if (arg[iarg][0] == 'i') which[nvalues] = INAME;
else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE; else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
int n = strlen(arg[iarg]); int n = strlen(arg[iarg]);
@ -241,11 +238,50 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
char *ptr = strchr(suffix,'['); char *ptr = strchr(suffix,'[');
if (ptr) { if (ptr) {
if (which[nvalues] == VARIABLE)
error->all(FLERR,"Illegal fix store/state command");
if (suffix[strlen(suffix)-1] != ']') if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal fix store/state command"); error->all(FLERR,"Illegal fix store/state command");
argindex[nvalues] = atoi(ptr+1); suffix[strlen(suffix)-1] = '\0';
argindex[nvalues] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0'; *ptr = '\0';
} else argindex[nvalues] = 0; } else argindex[nvalues] = 0;
n = strlen(suffix) + 1;
ids[nvalues] = new char[n];
strcpy(ids[nvalues],suffix);
nvalues++;
delete [] suffix;
} else if (strncmp(arg[iarg],"i_",2) == 0 ||
strncmp(arg[iarg],"d_",2) == 0 ||
strncmp(arg[iarg],"i2_",3) == 0 ||
strncmp(arg[iarg],"d2_",3) == 0) {
if (strncmp(arg[iarg],"i_",2) == 0) which[nvalues] = IVEC;
else if (strncmp(arg[iarg],"d_",2) == 0) which[nvalues] = DVEC;
else if (strncmp(arg[iarg],"i2_",3) == 0) which[nvalues] = IARRAY;
else if (strncmp(arg[iarg],"d2_",3) == 0) which[nvalues] = DARRAY;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
if (which[nvalues] == IVEC || which[nvalues] == DVEC)
strcpy(suffix,&arg[iarg][2]);
else strcpy(suffix,&arg[iarg][3]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (which[nvalues] == IVEC || which[nvalues] == DVEC)
error->all(FLERR,"Illegal fix store/state command");
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal fix store/state command");
suffix[strlen(suffix)-1] = '\0';
argindex[nvalues] = utils::inumeric(FLERR,ptr+1,true,lmp);
*ptr = '\0';
} else {
if (which[nvalues] == IARRAY || which[nvalues] == DARRAY)
error->all(FLERR,"Illegal fix store/state command");
argindex[nvalues] = 0;
}
n = strlen(suffix) + 1; n = strlen(suffix) + 1;
ids[nvalues] = new char[n]; ids[nvalues] = new char[n];
@ -274,68 +310,79 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
// error check // error check
for (int i = 0; i < nvalues; i++) { for (int m = 0; m < nvalues; m++) {
if (which[i] == COMPUTE) { if (which[m] == COMPUTE) {
int icompute = modify->find_compute(ids[i]); int icompute = modify->find_compute(ids[m]);
if (icompute < 0) if (icompute < 0)
error->all(FLERR,"Compute ID for fix store/state does not exist"); error->all(FLERR,"Compute ID for fix store/state does not exist");
if (modify->compute[icompute]->peratom_flag == 0) if (modify->compute[icompute]->peratom_flag == 0)
error->all(FLERR,"Fix store/state compute " error->all(FLERR,"Fix store/state compute "
"does not calculate per-atom values"); "does not calculate per-atom values");
if (argindex[i] == 0 && if (argindex[m] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0) modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Fix store/state compute does not " error->all(FLERR,"Fix store/state compute does not "
"calculate a per-atom vector"); "calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0) if (argindex[m] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR, error->all(FLERR,
"Fix store/state compute does not " "Fix store/state compute does not "
"calculate a per-atom array"); "calculate a per-atom array");
if (argindex[i] && if (argindex[m] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols) argindex[m] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR, error->all(FLERR,
"Fix store/state compute array is accessed out-of-range"); "Fix store/state compute array is accessed out-of-range");
} else if (which[i] == INAME) { } else if (which[m] == FIX) {
int icustom,iflag; int ifix = modify->find_fix(ids[m]);
icustom = atom->find_custom(ids[i],iflag);
if ((icustom < 0) || (iflag != 0))
error->all(FLERR,
"Custom integer vector for fix store/state does not exist");
} else if (which[i] == DNAME) {
int icustom,iflag;
icustom = atom->find_custom(ids[i],iflag);
if ((icustom < 0) || (iflag != 1))
error->all(FLERR,
"Custom floating point vector for fix store/state does not exist");
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0) if (ifix < 0)
error->all(FLERR, error->all(FLERR,
"Fix ID for fix store/state does not exist"); "Fix ID for fix store/state does not exist");
if (modify->fix[ifix]->peratom_flag == 0) if (modify->fix[ifix]->peratom_flag == 0)
error->all(FLERR, error->all(FLERR,
"Fix store/state fix does not calculate per-atom values"); "Fix store/state fix does not calculate per-atom values");
if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0) if (argindex[m] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR, error->all(FLERR,
"Fix store/state fix does not calculate a per-atom vector"); "Fix store/state fix does not calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0) if (argindex[m] && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR, error->all(FLERR,
"Fix store/state fix does not calculate a per-atom array"); "Fix store/state fix does not calculate a per-atom array");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols) if (argindex[m] && argindex[m] > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR, error->all(FLERR,
"Fix store/state fix array is accessed out-of-range"); "Fix store/state fix array is accessed out-of-range");
if (nevery % modify->fix[ifix]->peratom_freq) if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR, error->all(FLERR,
"Fix for fix store/state not computed at compatible time"); "Fix for fix store/state not computed at compatible time");
} else if (which[i] == VARIABLE) { } else if (which[m] == VARIABLE) {
int ivariable = input->variable->find(ids[i]); int ivariable = input->variable->find(ids[m]);
if (ivariable < 0) if (ivariable < 0)
error->all(FLERR,"Variable name for fix store/state does not exist"); error->all(FLERR,"Variable name for fix store/state does not exist");
if (input->variable->atomstyle(ivariable) == 0) if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Fix store/state variable is not atom-style variable"); error->all(FLERR,"Fix store/state variable is not atom-style variable");
} else if (which[m] == IVEC) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if ((icustom < 0) || flag || cols)
error->all(FLERR,
"Custom integer vector for fix store/state does not exist");
} else if (which[m] == DVEC) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if ((icustom < 0) || !flag || cols)
error->all(FLERR,
"Custom floating point vector for fix store/state does not exist");
} else if (which[m] == IARRAY) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if ((icustom < 0) || flag || !cols)
error->all(FLERR,
"Custom integer array for fix store/state does not exist");
} else if (which[m] == DARRAY) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if ((icustom < 0) || !flag || !cols)
error->all(FLERR,
"Custom floating point array for fix store/state does not exist");
} }
} }
@ -413,22 +460,6 @@ void FixStoreState::init()
error->all(FLERR,"Compute ID for fix store/state does not exist"); error->all(FLERR,"Compute ID for fix store/state does not exist");
value2index[m] = icompute; value2index[m] = icompute;
} else if (which[m] == INAME) {
int icustom,iflag;
icustom = atom->find_custom(ids[m],iflag);
if ((icustom < 0) || (iflag != 0))
error->all(FLERR,
"Custom integer vector for fix store/state does not exist");
value2index[m] = icustom;
} else if (which[m] == DNAME) {
int icustom,iflag;
icustom = atom->find_custom(ids[m],iflag);
if ((icustom < 0) || (iflag != 1))
error->all(FLERR,
"Custom floating point vector for fix store/state does not exist");
value2index[m] = icustom;
} else if (which[m] == FIX) { } else if (which[m] == FIX) {
int ifix = modify->find_fix(ids[m]); int ifix = modify->find_fix(ids[m]);
if (ifix < 0) if (ifix < 0)
@ -440,6 +471,35 @@ void FixStoreState::init()
if (ivariable < 0) if (ivariable < 0)
error->all(FLERR,"Variable name for fix store/state does not exist"); error->all(FLERR,"Variable name for fix store/state does not exist");
value2index[m] = ivariable; value2index[m] = ivariable;
} else if (which[m] == IVEC) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if (icustom < 0 || flag || cols)
error->all(FLERR,
"Custom integer vector for fix store/state does not exist");
value2index[m] = icustom;
} else if (which[m] == DVEC) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if (icustom < 0 || !flag || cols)
error->all(FLERR,
"Custom floating point vector for fix store/state does not exist");
value2index[m] = icustom;
} else if (which[m] == IARRAY) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if (icustom < 0 || flag || !cols)
error->all(FLERR,
"Custom integer array for fix store/state does not exist");
value2index[m] = icustom;
} else if (which[m] == DARRAY) {
int icustom,flag,cols;
icustom = atom->find_custom(ids[m],flag,cols);
if (icustom < 0 || !flag || !cols)
error->all(FLERR,
"Custom floating point array for fix store/state does not exist");
value2index[m] = icustom;
} }
} }
} }
@ -525,22 +585,33 @@ void FixStoreState::end_of_step()
if (mask[i] & groupbit) values[i][m] = fix_array[i][jm1]; if (mask[i] & groupbit) values[i][m] = fix_array[i][jm1];
} }
// access custom atom property fields
} else if (which[m] == INAME) {
int *ivector = atom->ivector[n];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) values[i][m] = ivector[i];
} else if (which[m] == DNAME) {
double *dvector = atom->dvector[n];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) values[i][m] = dvector[i];
// evaluate atom-style variable // evaluate atom-style variable
} else if (which[m] == VARIABLE) { } else if (which[m] == VARIABLE) {
input->variable->compute_atom(n,igroup,&values[0][m],nvalues,0); input->variable->compute_atom(n,igroup,&values[0][m],nvalues,0);
// access custom atom property fields
} else if (which[m] == IVEC) {
int *ivector = atom->ivector[n];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) values[i][m] = ivector[i];
} else if (which[m] == DVEC) {
double *dvector = atom->dvector[n];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) values[i][m] = dvector[i];
} else if (which[m] == IARRAY) {
int **iarray = atom->iarray[n];
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
values[i][m] = iarray[i][jm1];
} else if (which[m] == DARRAY) {
double **darray = atom->darray[n];
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
values[i][m] = darray[i][jm1];
} }
} }
} }

View File

@ -47,18 +47,17 @@ int ImbalanceStore::options(int narg, char **arg)
void ImbalanceStore::compute(double *weight) void ImbalanceStore::compute(double *weight)
{ {
int dflag = 0; int flag,cols;
int idx = atom->find_custom(name,dflag); int index = atom->find_custom(name,flag,cols);
// property does not exist // property does not exist
if (idx < 0 || dflag != 1) return; if (index < 0 || flag != 1 || cols)
error->all(FLERR,"Balance weight store vector does not exist");
double *prop = atom->dvector[idx]; double *prop = atom->dvector[index];
const int nlocal = atom->nlocal; const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; ++i) prop[i] = weight[i];
for (int i = 0; i < nlocal; ++i)
prop[i] = weight[i];
} }
/* -------------------------------------------------------------------- */ /* -------------------------------------------------------------------- */

View File

@ -600,8 +600,9 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
expand arg to earg, for arguments with syntax c_ID[*] or f_ID[*] expand arg to earg, for args with syntax c_ID[*], f_ID[*], i2_ID[*], d2_ID[*]
fields to consider in input arg range from iarg to narg fields to consider in input arg range from iarg to narg
mode = 0 for args = vectors, mode = 1 for args = arrays
return new expanded # of values, and copy them w/out "*" into earg return new expanded # of values, and copy them w/out "*" into earg
if any expansion occurs, earg is new allocation, must be freed by caller if any expansion occurs, earg is new allocation, must be freed by caller
if no expansion occurs, earg just points to arg, caller need not free if no expansion occurs, earg just points to arg, caller need not free
@ -609,7 +610,7 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
int Input::expand_args(int narg, char **arg, int mode, char **&earg) int Input::expand_args(int narg, char **arg, int mode, char **&earg)
{ {
int n,iarg,index,nlo,nhi,nmax,expandflag,icompute,ifix; int n,iarg,index,nlo,nhi,nmax,expandflag;
char *ptr1,*ptr2,*str; char *ptr1,*ptr2,*str;
ptr1 = NULL; ptr1 = NULL;
@ -633,7 +634,9 @@ int Input::expand_args(int narg, char **arg, int mode, char **&earg)
expandflag = 0; expandflag = 0;
if (strncmp(arg[iarg],"c_",2) == 0 || if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0) { strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"i2_",3) == 0 ||
strncmp(arg[iarg],"d2_",3) == 0) {
ptr1 = strchr(&arg[iarg][2],'['); ptr1 = strchr(&arg[iarg][2],'[');
if (ptr1) { if (ptr1) {
@ -641,9 +644,10 @@ int Input::expand_args(int narg, char **arg, int mode, char **&earg)
if (ptr2) { if (ptr2) {
*ptr2 = '\0'; *ptr2 = '\0';
if (strchr(ptr1,'*')) { if (strchr(ptr1,'*')) {
if (arg[iarg][0] == 'c') { if (arg[iarg][0] == 'c') {
*ptr1 = '\0'; *ptr1 = '\0';
icompute = modify->find_compute(&arg[iarg][2]); int icompute = modify->find_compute(&arg[iarg][2]);
*ptr1 = '['; *ptr1 = '[';
// check for global vector/array, peratom array, local array // check for global vector/array, peratom array, local array
@ -665,9 +669,10 @@ int Input::expand_args(int narg, char **arg, int mode, char **&earg)
expandflag = 1; expandflag = 1;
} }
} }
} else if (arg[iarg][0] == 'f') { } else if (arg[iarg][0] == 'f') {
*ptr1 = '\0'; *ptr1 = '\0';
ifix = modify->find_fix(&arg[iarg][2]); int ifix = modify->find_fix(&arg[iarg][2]);
*ptr1 = '['; *ptr1 = '[';
// check for global vector/array, peratom array, local array // check for global vector/array, peratom array, local array
@ -689,10 +694,40 @@ int Input::expand_args(int narg, char **arg, int mode, char **&earg)
expandflag = 1; expandflag = 1;
} }
} }
}
} } else if (arg[iarg][0] == 'i') {
*ptr2 = ']'; *ptr1 = '\0';
} int flag,cols;
int icustom = atom->find_custom(&arg[iarg][3],flag,cols);
*ptr1 = '[';
// check for custom per-atom integer array
if (icustom >= 0) {
if (mode == 1 && !flag && cols) {
nmax = cols;
expandflag = 1;
}
}
} else if (arg[iarg][0] == 'd') {
*ptr1 = '\0';
int flag,cols;
int icustom = atom->find_custom(&arg[iarg][3],flag,cols);
*ptr1 = '[';
// check for custom per-atom floating point array
if (icustom >= 0) {
if (mode == 1 && flag && cols) {
nmax = cols;
expandflag = 1;
}
}
}
}
*ptr2 = ']';
}
} }
} }

View File

@ -45,12 +45,12 @@ using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET, enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,
MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI, MOLECULE,X,Y,Z,VX,VY,VZ,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM, DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM,
THETA,THETA_RANDOM,ANGMOM,OMEGA, THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY, SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
SMD_CONTACT_RADIUS,DPDTHETA,INAME,DNAME,VX,VY,VZ}; SMD_CONTACT_RADIUS,DPDTHETA,IVEC,DVEC,IARRAY,DARRAY};
#define BIG INT_MAX #define BIG INT_MAX
@ -573,28 +573,63 @@ void Set::command(int narg, char **arg)
set(DPDTHETA); set(DPDTHETA);
iarg += 2; iarg += 2;
} else if (strstr(arg[iarg],"i_") == arg[iarg]) { // custom per-atom vector
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); } else if (strstr(arg[iarg],"i_") == arg[iarg] ||
else ivalue = force->inumeric(FLERR,arg[iarg+1]); strstr(arg[iarg],"d_") == arg[iarg]) {
int flag;
index_custom = atom->find_custom(&arg[iarg][2],flag);
if (index_custom < 0 || flag != 0)
error->all(FLERR,"Set command integer vector does not exist");
set(INAME);
iarg += 2;
} else if (strstr(arg[iarg],"d_") == arg[iarg]) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
int which = 0;
if (arg[iarg][0] == 'd') which = 1;
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else if (!which) ivalue = force->inumeric(FLERR,arg[iarg+1]);
else dvalue = force->numeric(FLERR,arg[iarg+1]); else dvalue = force->numeric(FLERR,arg[iarg+1]);
int flag;
index_custom = atom->find_custom(&arg[iarg][2],flag); int flag,cols;
if (index_custom < 0 || flag != 1) index_custom = atom->find_custom(&arg[iarg][2],flag,cols);
error->all(FLERR,"Set command floating point vector does not exist"); if ((!which && (index_custom < 0 || flag || cols)) ||
set(DNAME); (which && (index_custom < 0 || !flag || cols)))
error->all(FLERR,"Set command per-atom custom vector does not exist");
if (!which) set(IVEC);
else set(DVEC);
iarg += 2; iarg += 2;
// custom per-atom array, must include bracketed index
} else if (strstr(arg[iarg],"i2_") == arg[iarg] ||
strstr(arg[iarg],"d2_") == arg[iarg]) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
int which = 0;
if (arg[iarg][0] == 'd') which = 1;
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else if (!which) ivalue = force->inumeric(FLERR,arg[iarg+1]);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][3]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Invalid attribute in set command");
icol_custom = atoi(ptr+1);
*ptr = '\0';
} else error->all(FLERR,"Set command per-atom custom array is not indexed");
int flag,cols;
index_custom = atom->find_custom(suffix,flag,cols);
delete [] suffix;
if ((!which && (index_custom < 0 || flag || !cols)) ||
(which && (index_custom < 0 || !flag || !cols)))
error->all(FLERR,"Set command per-atom custom array does not exist");
if (icol_custom <= 0 || icol_custom > cols)
error->all(FLERR,"Set command per-atom custom array is accessed out-of-range");
if (!which) set(IARRAY);
else set(DARRAY);
iarg += 2;
} else error->all(FLERR,"Illegal set command"); } else error->all(FLERR,"Illegal set command");
// statistics // statistics
@ -970,20 +1005,29 @@ void Set::set(int keyword)
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
} }
// set value for custom integer or double vector // set value for custom vector or array
else if (keyword == INAME) { else if (keyword == IVEC) {
atom->ivector[index_custom][i] = ivalue; atom->ivector[index_custom][i] = ivalue;
} }
else if (keyword == DNAME) { else if (keyword == DVEC) {
atom->dvector[index_custom][i] = dvalue; atom->dvector[index_custom][i] = dvalue;
} }
else if (keyword == IARRAY) {
atom->iarray[index_custom][i][icol_custom-1] = ivalue;
}
else if (keyword == DARRAY) {
atom->darray[index_custom][i][icol_custom-1] = dvalue;
}
count++; count++;
} }
// update bonus data numbers // update bonus data numbers
if (keyword == SHAPE) { if (keyword == SHAPE) {
bigint nlocal_bonus = avec_ellipsoid->nlocal_bonus; bigint nlocal_bonus = avec_ellipsoid->nlocal_bonus;
MPI_Allreduce(&nlocal_bonus,&atom->nellipsoids,1, MPI_Allreduce(&nlocal_bonus,&atom->nellipsoids,1,

View File

@ -32,7 +32,7 @@ class Set : protected Pointers {
private: private:
char *id; char *id;
int *select; int *select;
int style,ivalue,newtype,count,index_custom; int style,ivalue,newtype,count,index_custom,icol_custom;
int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag; int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag;
int cc_index; int cc_index;
bigint nsubset; bigint nsubset;