diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index f8c0937e9f..700abe0c2b 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -219,6 +219,7 @@ OPT. * :doc:`rigid/small (o) ` * :doc:`rx (k) ` * :doc:`saed/vtk ` + * :doc:`set ` * :doc:`setforce (k) ` * :doc:`setforce/spin ` * :doc:`sgcmc ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index d49850e7da..bb0545e406 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -398,6 +398,7 @@ accelerated styles exist. * :doc:`rigid/small ` - constrain many small clusters of atoms to move as a rigid body with NVE integration * :doc:`rx ` - solve reaction kinetic ODEs for a defined reaction set * :doc:`saed/vtk ` - time-average the intensities from :doc:`compute saed ` +* :doc:`set ` - reset an atom property via an atom-style variable every N steps * :doc:`setforce ` - set the force on each atom * :doc:`setforce/spin ` - set magnetic precession vectors on each atom * :doc:`sgcmc ` - fix for hybrid semi-grand canonical MD/MC simulations diff --git a/doc/src/fix_set.rst b/doc/src/fix_set.rst new file mode 100644 index 0000000000..cfa8231ccd --- /dev/null +++ b/doc/src/fix_set.rst @@ -0,0 +1,173 @@ +.. index:: fix set + +fix set command +=============== + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID set Nfreq rnflag set-args + +* ID, group-ID are documented in :doc:`fix ` command +* set = style name of this fix command +* Nfreq = reset per-atom properties every this many timesteps +* rnflag = 1 to reneighbor on next timestep, 0 to not +* set-args = identical to args for the :doc:`set ` command + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 10 all set 1 0 group all i_dump v_new + fix 10 all set 1 0 group all i_dump v_turnoff + +Description +""""""""""" + +Reset one or more properties of one or more atoms once every *Nfreq* +steps during a simulation. + +If the *rnflag* for reneighboring is set to 1, then a reneighboring +will be triggered on the next timestep (since the fix set operation +occurs at the end of the current timestep). This is important to do +if this command changes per-atom properties that need to be +communicated to ghost atoms. If this is not the case, an *rnflag* +setting of 0 can be used; reneighboring will only be triggered on +subsequent timesteps by the usual neighbor list criteria; see the +:doc:`neigh_modify command `. + +Here are two examples where an *rnflag* setting of 1 are needed. If a +custom per-atom property is changed and the :doc:`fix property/atom +` command to create the property used the *ghost +yes* keyword. Or if per-atom charges are changed, all pair styles +which compute Coulombic interactions require charge values for ghost +atoms. In both these examples, the re-neighboring will trigger the +changes in the owned atom properties to be immediately communicated to +ghost atoms. + +The arguments following *Nfreq* and *rnflag* are identical to those +allowed for the :doc:`set ` command, as in the examples above and +below. + +Note that the group-ID setting for this command is ignored. The +syntax for the :doc:`set ` arguments allows selection of which +atoms have their properties reset. + +This command can only be used to reset an atom property using a +per-atom variable. This option in allowed by many, but not all, of +the keyword/value pairs supported by the :doc:`set ` command. +The reason for this restriction is that if a per-atom variable is not +used, this command will typically not change atom properties during +the simulation. + +The :doc:`set ` command can be used with similar syntax to this +command to reset atom properties once before or between simulations. + +---------- + +Here is an example of input script commands which will output atoms +into a dump file only when their x-velocity crosses a threshold value +*vthresh* for the first time. Their position and x-velocity will then +be output every step for *twindow* timesteps. + +.. code-block:: LAMMPS + + variable vthresh equal 2 # threshold velocity + variable twindow equal 10 # dump for this many steps + # + # define custom property i_dump to store timestep threshold is crossed + # + fix 2 all property/atom i_dump + set group all i_dump -1 + # + # fix set command checks for threshold crossings every step + # resets i_dump from -1 to current timestep when crossing occurs + # + variable start atom "vx > v_vthresh && i_dump == -1" + variable new atom ternary(v_start,step,i_dump) + fix 3 all set 1 0 group all i_dump v_new + # + # dump command with thresh which enforces twindow + # + dump 1 all custom 1 tmp.dump id x y vx i_dump + variable dumpflag atom "i_dump >= 0 && (step-i_dump) < v_twindow" + dump_modify 1 thresh v_dumpflag == 1 + # + # run the simulation + # final dump with all atom IDs which crossed threshold on which timestep + # + run 1000 + write_dump all custom tmp.dump.final id i_dump modify thresh i_dump >= 0 + +The tmp.dump.final file lists which atoms crossed the velocity +threshold. This command will print the *twindow* timesteps when a +specific atom ID (104 in this case) was output in the tmp.dump file: + +.. code-block:: LAMMPS + + % grep "^104 " tmp.dump + +If these commands are used instead of the above, then an atom can +cross the velocity threshold multiple times, and will be output for +*twindow* timesteps each time. However the write_dump command is no +longer useful. + +.. code-block:: LAMMPS + + variable vthresh equal 2 # threshold velocity + variable twindow equal 10 # dump for this many steps + # + # define custom property i_dump to store timestep threshold is crossed + # + fix 2 all property/atom i_dump + set group all i_dump -1 + # + # fix set command checks for threshold crossings every step + # resets i_dump from -1 to current timestep when crossing occurs + # + variable start atom "vx > v_vthresh && i_dump == -1" + variable turnon atom ternary(v_start,step,i_dump) + variable stop atom "v_turnon >= 0 && (step-v_turnon) < v_twindow" + variable turnoff atom ternary(v_stop,v_turnon,-1) + fix 3 all set 1 0 group all i_dump v_turnoff + # + # dump command with thresh which enforces twindow + # + dump 1 all custom 1 tmp.dump id x y vx i_dump + variable dumpflag atom "i_dump >= 0 && (step-i_dump) < v_twindow" + dump_modify 1 thresh v_dumpflag == 1 + # + # run the simulation + # + run 1000 + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. No global or per-atom quantities are stored by +this fix for access by various :doc:`output commands `. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. + +Restrictions +"""""""""""" + +none + +Related commands +"""""""""""""""" + +:doc:`set ` + +Default +""""""" + +none diff --git a/doc/src/run.rst b/doc/src/run.rst index c3dd171f6b..95e2f037db 100644 --- a/doc/src/run.rst +++ b/doc/src/run.rst @@ -103,14 +103,16 @@ must be done. .. note:: - If your input script changes the system between 2 runs, then the - initial setup must be performed to ensure the change is recognized by - all parts of the code that are affected. Examples are adding a - :doc:`fix ` or :doc:`dump ` or :doc:`compute `, changing - a :doc:`neighbor ` list parameter, or writing restart file - which can migrate atoms between processors. LAMMPS has no easy way to - check if this has happened, but it is an error to use the *pre no* - option in this case. + If your input script "changes" the system between 2 runs, then the + initial setup typically needs to be performed to ensure the change + is recognized by all parts of the code that are affected. Examples + are adding a :doc:`fix ` or :doc:`dump ` or + :doc:`compute `, changing a :doc:`neighbor ` + list parameter, using the :doc:`set ` command, or writing a + restart file via the :doc:`write_restart ` command, + which can migrate atoms between processors. LAMMPS has no easy way + to check if this has happened, but it is an error to use the *pre + no* option in these cases. If *post* is specified as "no", the full timing summary is skipped; only a one-line summary timing is printed. diff --git a/doc/src/set.rst b/doc/src/set.rst index 6be590e9b6..d2c865dfbb 100644 --- a/doc/src/set.rst +++ b/doc/src/set.rst @@ -22,21 +22,110 @@ Syntax for style = *region*, ID = a region ID * 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 *vx* or *vy* or *vz* or *charge* or - *dipole* or *dipole/random* or *quat* or *spin/atom* or *spin/atom/random* or - *spin/electron* or *radius/electron* 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 *temperature* 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 *epsilon* or - *i_name* or *d_name* or *i2_name* or *d2_name* + +* keyword = *angle* or *angmom* or *bond* or *cc* or *charge* or + *density* or *density/disc* or *diameter* or *dihedral* or *dipole* + or *dipole/random* or *dpd/theta* or *edpd/cv* or *edpd/temp* or + *epsilon* or *image* or *improper* or *length* or *mass* or *mol* or + *omega* or *quat* or *quat/random* or *radius/electron* or *shape* or + *smd/contact/radius* or *smd/mass/density* or *sph/cv* or *sph/e* or + *sph/rho* or *spin/atom* or *spin/atom/random* or *spin/electron* or + *temperature* or *theta* or *theta/random* or *tri* or *type* or + *type/fraction* or *type/ratio* or *type/subset* or *volume* or *vx* + or *vy* or *vz* or *x* or *y* or *z* or *i_name* or *d_name* or + *i2_name* or *d2_name* .. parsed-literal:: + *angle* value = numeric angle type or angle type label, for all angles between selected atoms + *angmom* values = Lx Ly Lz + Lx,Ly,Lz = components of angular momentum vector (distance-mass-velocity units) + any of Lx,Ly,Lz can be an atom-style variable (see below) + *bond* value = numeric bond type or bond type label, for all bonds between selected atoms + *cc* values = index cc + index = index of a chemical species (1 to Nspecies) + cc = chemical concentration of tDPD particles for a species (mole/volume units) + cc can be an atom-style variable (see below) + *charge* value = atomic charge (charge units) + value can be an atom-style variable (see below) + *density* value = particle density for a sphere or ellipsoid (mass/distance\^3 units), or for a triangle (mass/distance\^2 units) or line (mass/distance units) particle + value can be an atom-style variable (see below) + *density/disc* value = particle density for a 2d disc or ellipse (mass/distance\^2 units) + value can be an atom-style variable (see below) + *diameter* value = diameter of spherical particle (distance units) + value can be an atom-style variable (see below) + *dihedral* value = numeric dihedral type or dihedral type label, for all dihedrals between selected atoms + *dipole* values = x y z + x,y,z = orientation of dipole moment vector + any of x,y,z can be an atom-style variable (see below) + *dipole/random* values = seed Dlen + seed = random # seed (positive integer) for dipole moment orientations + Dlen = magnitude of dipole moment (dipole units) + *dpd/theta* value = internal temperature of DPD particles (temperature units) + value can be an atom-style variable (see below) + value can be NULL which sets internal temp of each particle to KE temp + *edpd/cv* value = volumetric heat capacity of eDPD particles (energy/temperature/volume units) + value can be an atom-style variable (see below) + *edpd/temp* value = temperature of eDPD particles (temperature units) + value can be an atom-style variable (see below) + *epsilon* value = dielectric constant of the medium where the atoms reside + value can be an atom-style variable (see below) + *image* values = nx ny nz + nx,ny,nz = which periodic image of the simulation box the atom is in + any of nx,ny,nz can be an atom-style variable (see below) + *improper* value = numeric improper type or improper type label, for all impropers between selected atoms + *length* value = len + len = length of line segment (distance units) + len can be an atom-style variable (see below) + *mass* value = per-atom mass (mass units) + value can be an atom-style variable (see below) + *mol* value = molecule ID + the moleclue ID can be an atom-style variable (see below) + *omega* values = Wx Wy Wz + Wx,Wy,Wz = components of angular velocity vector (radians/time units) + any of Wx,Wy,Wz can be an atom-style variable (see below) + *quat* values = a b c theta + a,b,c = unit vector to rotate particle around via right-hand rule + theta = rotation angle (degrees) + any of a,b,c,theta values can be an atom-style variable (see below) + *quat/random* value = seed + seed = random # seed (positive integer) for quaternion orientations + *radius/electron* value = eradius + eradius = electron radius (or fixed-core radius) (distance units) + value can be an atom-style variable (see below) + *shape* values = Sx Sy Sz + Sx,Sy,Sz = 3 diameters of ellipsoid (distance units) + any of Sx,Sy,Sz can be an atom-style variable (see below) + *smd/contact/radius* value = radius for short range interactions, i.e. contact and friction + value can be an atom-style variable (see below) + *smd/mass/density* value = set particle mass based on volume by providing a mass density + value can be an atom-style variable (see below) + *sph/cv* value = heat capacity of SPH particles (need units) + value can be an atom-style variable (see below) + *sph/e* value = energy of SPH particles (need units) + value can be an atom-style variable (see below) + *sph/rho* value = density of SPH particles (need units) + value can be an atom-style variable (see below) + *spin/atom* values = g x y z + g = magnitude of magnetic spin vector (in Bohr magneton's unit) + x,y,z = orientation of magnetic spin vector + any of x,y,z can be an atom-style variable (see below) + *spin/atom/random* values = seed Dlen + seed = random # seed (positive integer) for magnetic spin orientations + Dlen = magnitude of magnetic spin vector (in Bohr magneton's unit) + *spin/electron* value = espin + espin = electron spin (+1/-1), 0 = nuclei, 2 = fixed-core, 3 = pseudo-cores (i.e. ECP) + value can be an atom-style variable (see below) + *temperature* value = temperature for finite-size particles (temperature units) + value can be an atom-style variable (see below) + *theta* value = angle (degrees) + angle = orientation of line segment with respect to x-axis + value can be an atom-style variable (see below) + *theta/random* value = seed + seed = random # seed (positive integer) for line segment orienations + *tri* value = side + side = side length of equilateral triangle (distance units) + value can be an atom-style variable (see below) *type* value = numeric atom type or type label value can be an atom-style variable (see below) *type/fraction* values = type fraction seed @@ -51,104 +140,22 @@ Syntax type = numeric atom type or type label Nsubset = exact number of selected atoms to set to new atom type seed = random # seed (positive integer) - *mol* value = molecule ID - value can be an atom-style variable (see below) - *x*,\ *y*,\ *z* value = atom coordinate (distance units) + *volume* value = particle volume for Peridynamic particle (distance\^3 units) value can be an atom-style variable (see below) *vx*,\ *vy*,\ *vz* value = atom velocity (velocity units) value can be an atom-style variable (see below) - *charge* value = atomic charge (charge units) + *x*,\ *y*,\ *z* value = atom coordinate (distance units) value can be an atom-style variable (see below) - *dipole* values = x y z - x,y,z = orientation of dipole moment vector - any of x,y,z can be an atom-style variable (see below) - *dipole/random* value = seed Dlen - seed = random # seed (positive integer) for dipole moment orientations - Dlen = magnitude of dipole moment (dipole units) - *spin/atom* values = g x y z - g = magnitude of magnetic spin vector (in Bohr magneton's unit) - x,y,z = orientation of magnetic spin vector - any of x,y,z can be an atom-style variable (see below) - *spin/atom/random* value = seed Dlen - seed = random # seed (positive integer) for magnetic spin orientations - Dlen = magnitude of magnetic spin vector (in Bohr magneton's unit) - *radius/electron* values = eradius - eradius = electron radius (or fixed-core radius) (distance units) - *spin/electron* value = espin - espin = electron spin (+1/-1), 0 = nuclei, 2 = fixed-core, 3 = pseudo-cores (i.e. ECP) - *quat* values = a b c theta - a,b,c = unit vector to rotate particle around via right-hand rule - theta = rotation angle (degrees) - any of a,b,c,theta can be an atom-style variable (see below) - *quat/random* value = seed - seed = random # seed (positive integer) for quaternion orientations - *diameter* value = diameter of spherical particle (distance units) - value can be an atom-style variable (see below) - *shape* value = Sx Sy Sz - Sx,Sy,Sz = 3 diameters of ellipsoid (distance units) - *length* value = len - len = length of line segment (distance units) - len can be an atom-style variable (see below) - *tri* value = side - side = side length of equilateral triangle (distance units) - side can be an atom-style variable (see below) - *theta* value = angle (degrees) - angle = orientation of line segment with respect to x-axis - angle can be an atom-style variable (see below) - *theta/random* value = seed - seed = random # seed (positive integer) for line segment orienations - *angmom* values = Lx Ly Lz - Lx,Ly,Lz = components of angular momentum vector (distance-mass-velocity units) - any of Lx,Ly,Lz can be an atom-style variable (see below) - *omega* values = Wx Wy Wz - Wx,Wy,Wz = components of angular velocity vector (radians/time units) - any of wx,wy,wz can be an atom-style variable (see below) - *mass* value = per-atom mass (mass units) - value can be an atom-style variable (see below) - *density* value = particle density for a sphere or ellipsoid (mass/distance\^3 units), or for a triangle (mass/distance\^2 units) or line (mass/distance units) particle - value can be an atom-style variable (see below) - *density/disc* value = particle density for a 2d disc or ellipse (mass/distance\^2 units) - value can be an atom-style variable (see below) - *temperature* value = temperature for finite-size particles (temperature units) - value can be an atom-style variable (see below) - *volume* value = particle volume for Peridynamic particle (distance\^3 units) - value can be an atom-style variable (see below) - *image* nx ny nz - nx,ny,nz = which periodic image of the simulation box the atom is in - any of nx,ny,nz can be an atom-style variable (see below) - *bond* value = numeric bond type or bond type label, for all bonds between selected atoms - *angle* value = numeric angle type or angle type label, for all angles between selected atoms - *dihedral* value = numeric dihedral type or dihedral type label, for all dihedrals between selected atoms - *improper* value = numeric improper type or improper type label, for all impropers between selected atoms - *rheo/rho* value = density of RHEO particles (mass/distance\^3) - *rheo/status* value = status or phase of RHEO particles (unitless) - *sph/e* value = energy of SPH particles (need units) - value can be an atom-style variable (see below) - *sph/cv* value = heat capacity of SPH particles (need units) - value can be an atom-style variable (see below) - *sph/rho* value = density of SPH particles (need units) - value can be an atom-style variable (see below) - *smd/contact/radius* = radius for short range interactions, i.e. contact and friction - value can be an atom-style variable (see below) - *smd/mass/density* = set particle mass based on volume by providing a mass density - value can be an atom-style variable (see below) - *dpd/theta* value = internal temperature of DPD particles (temperature units) - value can be an atom-style variable (see below) - value can be NULL which sets internal temp of each particle to KE temp - *edpd/temp* value = temperature of eDPD particles (temperature units) - value can be an atom-style variable (see below) - *edpd/cv* value = volumetric heat capacity of eDPD particles (energy/temperature/volume units) - value can be an atom-style variable (see below) - *cc* values = index cc - index = index of a chemical species (1 to Nspecies) - cc = chemical concentration of tDPD particles for a species (mole/volume units) - *epsilon* value = dielectric constant of the medium where the atoms reside *i_name* value = custom integer vector with name + value can be an atom-style variable (see below) *d_name* value = custom floating-point vector with name - *i2_name* value = column of a custom integer array with name + value can be an atom-style variable (see below) + *i2_name* value = 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 + value can be an atom-style variable (see below) + *d2_name* value = custom floating-point array with name column specified as d2_name[N] where N is 1 to Ncol + value can be an atom-style variable (see below) Examples """""""" @@ -177,22 +184,26 @@ Description Set one or more properties of one or more atoms. Since atom properties are initially assigned by the :doc:`read_data `, -:doc:`read_restart ` or :doc:`create_atoms ` -commands, this command changes those assignments. This can be useful -for overriding the default values assigned by the -:doc:`create_atoms ` command (e.g. charge = 0.0). It can -be useful for altering pairwise and molecular force interactions, +:doc:`read_restart ` or :doc:`create_atoms +` commands, this command changes those assignments. +This can be useful for overriding the default values assigned by the +:doc:`create_atoms ` command (e.g. charge = 0.0). It +can be useful for altering pairwise and molecular force interactions, since force-field coefficients are defined in terms of types. It can be used to change the labeling of atoms by atom type or molecule ID -when they are output in :doc:`dump ` files. It can also be useful -for debugging purposes; i.e. positioning an atom at a precise location -to compute subsequent forces or energy. +when they are output in :doc:`dump ` files. It can also be +useful for debugging purposes; i.e. positioning an atom at a precise +location to compute subsequent forces or energy. Note that the *style* and *ID* arguments determine which atoms have their properties reset. The remaining keywords specify which properties to reset and what the new values are. Some strings like *type* or *mol* can be used as a style and/or a keyword. +The :doc:`fix set ` command can be used with similar syntax +to this command to reset atom properties once every *N* steps during a +simulation using via atom-style variables. + ---------- This section describes how to select which atoms to change @@ -211,8 +222,8 @@ can be specified, e.g. "C". The style *mol* selects all the atoms in a range of molecule IDs. In each of the range cases, the range can be specified as a single -numeric value, or a wildcard asterisk can be used to specify a range -of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For +numeric value, or with a wildcard asterisk to specify a range of +values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For example, for the style *type*, if N = the number of atom types, then an asterisk with no numeric values means all types from 1 to N. A leading asterisk means all types from 1 to n (inclusive). A trailing @@ -222,25 +233,25 @@ means all types from m to n (inclusive). For all the styles except The style *group* selects all the atoms in the specified group. The style *region* selects all the atoms in the specified geometric -region. See the :doc:`group ` and :doc:`region ` commands -for details of how to specify a group or region. +region. See the :doc:`group ` and :doc:`region ` +commands for details of how to specify a group or region. ---------- -This section describes the keyword options for which properties to +The next section describes the keyword options for which properties to change, for the selected atoms. Note that except where explicitly prohibited below, all of the keywords allow an :doc:`atom-style or atomfile-style variable ` to be used as the specified value(s). If the value is a -variable, it should be specified as v_name, where name is the -variable name. In this case, the variable will be evaluated, and its -resulting per-atom value used to determine the value assigned to each -selected atom. Note that the per-atom value from the variable will be -ignored for atoms that are not selected via the *style* and *ID* -settings explained above. A simple way to use per-atom values from -the variable to reset a property for all atoms is to use style *atom* -with *ID* = "\*"; this selects all atom IDs. +variable, it should be specified as v_name, where name is the variable +name. In this case, the variable will be evaluated, and its resulting +per-atom value used to determine the value assigned to each selected +atom. Note that the per-atom value from the variable will be ignored +for atoms that are not selected via the *style* and *ID* settings +explained above. A simple way to use per-atom values from the +variable to reset a property for all atoms is to use style *atom* with +*ID* = "\*"; this selects all atom IDs. Atom-style variables can specify formulas with various mathematical functions, and include :doc:`thermo_style ` command @@ -256,52 +267,110 @@ from a file. .. note:: Atom-style and atomfile-style variables return floating point - per-atom values. If the values are assigned to an integer variable, - such as the molecule ID, then the floating point value is truncated to - its integer portion, e.g. a value of 2.6 would become 2. + per-atom values. If the values are assigned to an integer + variable, such as the molecule ID, then the floating point value is + truncated to its integer portion, e.g. a value of 2.6 would + become 2. + +---------- .. versionchanged:: 28Mar2023 - Support for type labels was added for setting atom, bond, angle, - dihedral, and improper types + Support for type labels was added for setting angle types -Keyword *type* sets the atom type for all selected atoms. A specified -value can be either a numeric atom type or an atom type label. When -using a numeric type, the specified value must be from 1 to ntypes, -where ntypes was set by the :doc:`create_box ` command or -the *atom types* field in the header of the data file read by the -:doc:`read_data ` command. When using a type label it must -have been defined previously. See the :doc:`Howto type labels -` doc page for the allowed syntax of type labels -and a general discussion of how type labels can be used. +Keyword *angle* sets the angle type of all angles of selected atoms to +the specified value. The value can be a numeric type from 1 to +nangletypes. Or it can be a angle type label. See the :doc:`Howto +type labels ` doc page for the allowed syntax of +type labels and a general discussion of how type labels can be used. +All atoms in a particular angle must be selected atoms in order for +the change to be made. The value of nangletypes was set by the *angle +types* field in the header of the data file read by the +:doc:`read_data ` command. This keyword does NOT allow use +of an atom-style variable. -Keyword *type/fraction* sets the atom type for a fraction of the selected -atoms. The actual number of atoms changed is not guaranteed -to be exactly the specified fraction (0 <= *fraction* <= 1), but -should be statistically close. Random numbers are used in such a way -that a particular atom is changed or not changed, regardless of how -many processors are being used. This keyword does not allow use of an -atom-style variable. +Keyword *angmom* sets the angular momentum of selected atoms. The +particles must be ellipsoids as defined by the :doc:`atom_style +ellipsoid ` command or triangles as defined by the +:doc:`atom_style tri ` command. The angular momentum +vector of the particles is set to the 3 specified components. -Keywords *type/ratio* and *type/subset* also set the atom type for a -fraction of the selected atoms. The actual number of atoms changed -will be exactly the requested number. For *type/ratio* the specified -fraction (0 <= *fraction* <= 1) determines the number. For -*type/subset*, the specified *Nsubset* is the number. An iterative -algorithm is used which ensures the correct number of atoms are -selected, in a perfectly random fashion. Which atoms are selected -will change with the number of processors used. These keywords do not -allow use of an atom-style variable. +.. versionchanged:: 28Mar2023 -Keyword *mol* sets the molecule ID for all selected atoms. The -:doc:`atom style ` being used must support the use of -molecule IDs. + Support for type labels was added for setting bond types -Keywords *x*, *y*, *z*, and *charge* set the coordinates or -charge of all selected atoms. For *charge*, the :doc:`atom style -` being used must support the use of atomic -charge. Keywords *vx*, *vy*, and *vz* set the velocities of all -selected atoms. +Keyword *bond* sets the bond type of all bonds of selected atoms to +the specified value. The value can be a numeric type from 1 to +nbondtypes. Or it can be a bond type label. See the :doc:`Howto type +labels ` doc page for the allowed syntax of type +labels and a general discussion of how type labels can be used. All +atoms in a particular bond must be selected atoms in order for the +change to be made. The value of nbondtypes was set by the *bond +types* field in the header of the data file read by the +:doc:`read_data ` command. This keyword does NOT allow use +of an atom-style variable. + +Keyword *cc* sets the chemical concentration of a tDPD particle for a +specified species as defined by the DPD-MESO package. Currently, only +:doc:`atom_style tdpd ` defines particles with this +attribute. An integer for "index" selects a chemical species (1 to +Nspecies) where Nspecies is set by the atom_style command. The value +for the chemical concentration must be >= 0.0. + +Keyword *charge* set the charge of all selected atoms. The :doc:`atom +style ` being used must support the use of atomic charge. + +Keyword *density* or *density/disc* also sets the mass of all selected +particles, but in a different way. The particles must have a per-atom +mass attribute, as defined by the :doc:`atom_style ` +command. If the atom has a radius attribute (see :doc:`atom_style +sphere `) and its radius is non-zero, its mass is set from +the density and particle volume for 3d systems (the input density is +assumed to be in mass/distance\^3 units). For 2d, the default is for +LAMMPS to model particles with a radius attribute as spheres. +However, if the *density/disc* keyword is used, then they can be +modeled as 2d discs (circles). Their mass is set from the density and +particle area (the input density is assumed to be in mass/distance\^2 +units). + +If the atom has a shape attribute (see :doc:`atom_style ellipsoid +`) and its 3 shape parameters are non-zero, then its mass +is set from the density and particle volume (the input density is +assumed to be in mass/distance\^3 units). The *density/disc* keyword +has no effect; it does not (yet) treat 3d ellipsoids as 2d ellipses. + +If the atom has a length attribute (see :doc:`atom_style line +`) and its length is non-zero, then its mass is set from +the density and line segment length (the input density is assumed to +be in mass/distance units). If the atom has an area attribute (see +:doc:`atom_style tri `) and its area is non-zero, then its +mass is set from the density and triangle area (the input density is +assumed to be in mass/distance\^2 units). + +If none of these cases are valid, then the mass is set to the density +value directly (the input density is assumed to be in mass units). + +Keyword *diameter* sets the size of the selected atoms. The particles +must be finite-size spheres as defined by the :doc:`atom_style sphere +` command. The diameter of a particle can be set to 0.0, +which means they will be treated as point particles. Note that this +command does not adjust the particle mass, even if it was defined with +a density, e.g. via the :doc:`read_data ` command. + +.. versionchanged:: 28Mar2023 + + Support for type labels was added for setting dihedral types + +Keyword *dihedral* sets the dihedral type of all dihedrals of selected +atoms to the specified value. The value can be a numeric type from 1 +to ndihedraltypes. Or it can be a dihedral type label. See the +:doc:`Howto type labels ` doc page for the allowed +syntax of type labels and a general discussion of how type labels can +be used. All atoms in a particular dihedral must be selected atoms in +order for the change to be made. The value of ndihedraltypes was set +by the *dihedral types* field in the header of the data file read by +the :doc:`read_data ` command. This keyword does NOT allow +use of an atom-style variable. Keyword *dipole* uses the specified x,y,z values as components of a vector to set as the orientation of the dipole moment vectors of the @@ -313,40 +382,106 @@ moment vectors for the selected atoms and sets the magnitude of each to the specified *Dlen* value. For 2d systems, the z component of the orientation is set to 0.0. Random numbers are used in such a way that the orientation of a particular atom is the same, regardless of how -many processors are being used. This keyword does not allow use of an +many processors are being used. This keyword does NOT allow use of an atom-style variable. -.. versionchanged:: 15Sep2022 +Keyword *dpd/theta* sets the internal temperature of a DPD particle as +defined by the DPD-REACT package. If the specified value is a number +it must be >= 0.0. If the specified value is NULL, then the kinetic +temperature Tkin of each particle is computed as 3/2 k Tkin = KE = 1/2 +m v\^2 = 1/2 m (vx\*vx+vy\*vy+vz\*vz). Each particle's internal +temperature is set to Tkin. If the specified value is an atom-style +variable, then the variable is evaluated for each particle. If a +value >= 0.0, the internal temperature is set to that value. If it is +< 0.0, the computation of Tkin is performed and the internal +temperature is set to that value. -Keyword *spin/atom* uses the specified g value to set the magnitude of the -magnetic spin vectors, and the x,y,z values as components of a vector -to set as the orientation of the magnetic spin vectors of the selected -atoms. This keyword was previously called *spin*. +Keywords *edpd/temp* and *edpd/cv* set the temperature and volumetric +heat capacity of an eDPD particle as defined by the DPD-MESO package. +Currently, only :doc:`atom_style edpd ` defines particles +with these attributes. The values for the temperature and heat +capacity must be positive. -.. versionchanged:: 15Sep2022 +Keyword *epsilon* sets the dielectric constant of a particle to be +that of the medium where the particle resides as defined by the +DIELECTRIC package. Currently, only :doc:`atom_style dielectric +` defines particles with this attribute. The value for the +dielectric constant must be >= 0.0. Note that the set command with +this keyword will rescale the particle charge accordingly so that the +real charge (e.g., as read from a data file) stays intact. To change +the real charges, one needs to use the set command with the *charge* +keyword. Care must be taken to ensure that the real and scaled charges +and the dielectric constants are consistent. -Keyword *spin/atom/random* randomizes the orientation of the magnetic spin -vectors for the selected atoms and sets the magnitude of each to the -specified *Dlen* value. This keyword was previously called *spin/random*. +Keyword *image* sets which image of the simulation box the atom is +considered to be in. An image of 0 means it is inside the box as +defined. A value of 2 means add 2 box lengths to get the true value. +A value of -1 means subtract 1 box length to get the true value. +LAMMPS updates these flags as atoms cross periodic boundaries during +the simulation. The flags can be output with atom snapshots via the +:doc:`dump ` command. If a value of NULL is specified for any +of nx,ny,nz, then the current image value for that dimension is +unchanged. For non-periodic dimensions only a value of 0 can be +specified. This command can be useful after a system has been +equilibrated and atoms have diffused one or more box lengths in +various directions. This command can then reset the image values for +atoms so that they are effectively inside the simulation box, e.g if a +diffusion coefficient is about to be measured via the :doc:`compute +msd ` command. Care should be taken not to reset the +image flags of two atoms in a bond to the same value if the bond +straddles a periodic boundary (rather they should be different by +/- +1). This will not affect the dynamics of a simulation, but may mess +up analysis of the trajectories if a LAMMPS diagnostic or your own +analysis relies on the image flags to unwrap a molecule which +straddles the periodic box. -.. versionadded:: 15Sep2022 +.. versionchanged:: 28Mar2023 -Keyword *radius/electron* uses the specified value to set the radius of -electrons or fixed cores. + Support for type labels was added for setting improper types -.. versionadded:: 15Sep2022 +Keyword *improper* sets the improper type of all impropers of selected +atoms to the specified value. The value can be a numeric type from 1 +to nimpropertypes. Or it can be a improper type label. See the +:doc:`Howto type labels ` doc page for the allowed +syntax of type labels and a general discussion of how type labels can +be used. All atoms in a particular improper must be selected atoms in +order for the change to be made. The value of nimpropertypes was set +by the *improper types* field in the header of the data file read by +the :doc:`read_data ` command. This keyword does NOT allow +use of an atom-style variable. -Keyword *spin/electron* sets the spin of an electron (+/- 1) or indicates -nuclei (=0), fixed-cores (=2), or pseudo-cores (= 3). +Keyword *length* sets the length of selected atoms. The particles +must be line segments as defined by the :doc:`atom_style line +` command. If the specified value is non-zero the line +segment is (re)set to a length = the specified value, centered around +the particle position, with an orientation along the x-axis. If the +specified value is 0.0, the particle will become a point particle. +Note that this command does not adjust the particle mass, even if it +was defined with a density, e.g. via the :doc:`read_data ` +command. + +Keyword *mass* sets the mass of all selected particles. The particles +must have a per-atom mass attribute, as defined by the +:doc:`atom_style ` command. See the "mass" command for +how to set mass values on a per-type basis. + +Keyword *mol* sets the molecule ID for all selected atoms. The +:doc:`atom style ` being used must support the use of +molecule IDs. + +Keyword *omega* sets the angular velocity of selected atoms. The +particles must be spheres as defined by the :doc:`atom_style sphere +` command. The angular velocity vector of the particles +is set to the 3 specified components. Keyword *quat* uses the specified values to create a quaternion (4-vector) that represents the orientation of the selected atoms. The particles must define a quaternion for their orientation (e.g. ellipsoids, triangles, body particles) as defined by the -:doc:`atom_style ` command. Note that particles defined by -:doc:`atom_style ellipsoid ` have 3 shape parameters. The 3 -values must be non-zero for each particle set by this command. They -are used to specify the aspect ratios of an ellipsoidal particle, +:doc:`atom_style ` command. Note that particles defined +by :doc:`atom_style ellipsoid ` have 3 shape parameters. +The 3 values must be non-zero for each particle set by this command. +They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box's x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector @@ -360,51 +495,77 @@ ignored, since a rotation vector of (0,0,1) is the only valid choice. Keyword *quat/random* randomizes the orientation of the quaternion for the selected atoms. The particles must define a quaternion for their orientation (e.g. ellipsoids, triangles, body particles) as defined by -the :doc:`atom_style ` command. Random numbers are used in -such a way that the orientation of a particular atom is the same, +the :doc:`atom_style ` command. Random numbers are used +in such a way that the orientation of a particular atom is the same, regardless of how many processors are being used. For 2d systems, only orientations in the xy plane are generated. As with keyword *quat*, for ellipsoidal particles, the 3 shape values must be non-zero -for each particle set by this command. This keyword does not allow +for each particle set by this command. This keyword does NOT allow use of an atom-style variable. -Keyword *diameter* sets the size of the selected atoms. The particles -must be finite-size spheres as defined by the :doc:`atom_style sphere -` command. The diameter of a particle can be set to 0.0, -which means they will be treated as point particles. Note that this -command does not adjust the particle mass, even if it was defined with -a density, e.g. via the :doc:`read_data ` command. +.. versionadded:: 15Sep2022 + +Keyword *radius/electron* uses the specified value to set the radius +of electrons or fixed cores. Keyword *shape* sets the size and shape of the selected atoms. The particles must be ellipsoids as defined by the :doc:`atom_style -ellipsoid ` command. The *Sx*, *Sy*, *Sz* settings -are the 3 diameters of the ellipsoid in each direction. All 3 can be -set to the same value, which means the ellipsoid is effectively a -sphere. They can also all be set to 0.0 which means the particle will -be treated as a point particle. Note that this command does not -adjust the particle mass, even if it was defined with a density, -e.g. via the :doc:`read_data ` command. +ellipsoid ` command. The *Sx*, *Sy*, *Sz* settings are +the 3 diameters of the ellipsoid in each direction. All 3 can be set +to the same value, which means the ellipsoid is effectively a sphere. +They can also all be set to 0.0 which means the particle will be +treated as a point particle. Note that this command does not adjust +the particle mass, even if it was defined with a density, e.g. via the +:doc:`read_data ` command. -Keyword *length* sets the length of selected atoms. The particles -must be line segments as defined by the :doc:`atom_style line -` command. If the specified value is non-zero the line -segment is (re)set to a length = the specified value, centered around -the particle position, with an orientation along the x-axis. If the -specified value is 0.0, the particle will become a point particle. -Note that this command does not adjust the particle mass, even if it -was defined with a density, e.g. via the :doc:`read_data ` -command. +Keyword *smd/contact/radius* only applies to simulations with the +Smooth Mach Dynamics package MACHDYN. Itsets an interaction radius +for computing short-range interactions, e.g. repulsive forces to +prevent different individual physical bodies from penetrating each +other. Note that the SPH smoothing kernel diameter used for computing +long range, nonlocal interactions, is set using the *diameter* +keyword. -Keyword *tri* sets the size of selected atoms. The particles must be -triangles as defined by the :doc:`atom_style tri ` command. -If the specified value is non-zero the triangle is (re)set to be an -equilateral triangle in the xy plane with side length = the specified -value, with a centroid at the particle position, with its base -parallel to the x axis, and the y-axis running from the center of the -base to the top point of the triangle. If the specified value is 0.0, -the particle will become a point particle. Note that this command -does not adjust the particle mass, even if it was defined with a -density, e.g. via the :doc:`read_data ` command. +Keyword *smd/mass/density* sets the mass of all selected particles, +but it is only applicable to the Smooth Mach Dynamics package MACHDYN. +It assumes that the particle volume has already been correctly set and +calculates particle mass from the provided mass density value. + +Keywords *sph/cv*, *sph/e*, and *sph/rho* set the heat capacity, +energy, and density of smoothed particle hydrodynamics (SPH) +particles. See `this PDF guide `_ to +using SPH in LAMMPS. + +.. note:: + + Note that the SPH PDF guide file has not been updated for many + years and thus does not reflect the current *syntax* of the SPH + package commands. For that, please refer to the LAMMPS manual. + +.. versionchanged:: 15Sep2022 + +Keyword *spin/atom* uses the specified g value to set the magnitude of +the magnetic spin vectors, and the x,y,z values as components of a +vector to set as the orientation of the magnetic spin vectors of the +selected atoms. This keyword was previously called *spin*. + +.. versionchanged:: 15Sep2022 + +Keyword *spin/atom/random* randomizes the orientation of the magnetic +spin vectors for the selected atoms and sets the magnitude of each to +the specified *Dlen* value. This keyword does NOT allow use of an +atom-style variable. This keyword was previously called +*spin/random*. + +.. versionadded:: 15Sep2022 + +Keyword *spin/electron* sets the spin of an electron (+/- 1) or +indicates nuclei (=0), fixed-cores (=2), or pseudo-cores (= 3). + +Keyword *temperature* sets the temperature of a finite-size particle. +Currently, only the GRANULAR package supports this attribute. The +temperature must be added using an instance of :doc:`fix property/atom +` The values for the temperature must be positive. Keyword *theta* sets the orientation of selected atoms. The particles must be line segments as defined by the :doc:`atom_style line @@ -413,169 +574,71 @@ orientation angle of the line segments with respect to the x axis. Keyword *theta/random* randomizes the orientation of theta for the selected atoms. The particles must be line segments as defined by the -:doc:`atom_style line ` command. Random numbers are used in -such a way that the orientation of a particular atom is the same, +:doc:`atom_style line ` command. Random numbers are used +in such a way that the orientation of a particular atom is the same, regardless of how many processors are being used. This keyword does -not allow use of an atom-style variable. +NOT allow use of an atom-style variable. -Keyword *angmom* sets the angular momentum of selected atoms. The -particles must be ellipsoids as defined by the :doc:`atom_style -ellipsoid ` command or triangles as defined by the -:doc:`atom_style tri ` command. The angular momentum -vector of the particles is set to the 3 specified components. +Keyword *tri* sets the size of selected atoms. The particles must be +triangles as defined by the :doc:`atom_style tri ` +command. If the specified value is non-zero the triangle is (re)set +to be an equilateral triangle in the xy plane with side length = the +specified value, with a centroid at the particle position, with its +base parallel to the x axis, and the y-axis running from the center of +the base to the top point of the triangle. If the specified value is +0.0, the particle will become a point particle. Note that this +command does not adjust the particle mass, even if it was defined with +a density, e.g. via the :doc:`read_data ` command. -Keyword *omega* sets the angular velocity of selected atoms. The -particles must be spheres as defined by the :doc:`atom_style sphere -` command. The angular velocity vector of the particles is -set to the 3 specified components. +.. versionchanged:: 28Mar2023 -Keyword *mass* sets the mass of all selected particles. The particles -must have a per-atom mass attribute, as defined by the :doc:`atom_style -` command. See the "mass" command for how to set mass -values on a per-type basis. + Support for type labels was added for setting atom types -Keyword *density* or *density/disc* also sets the mass of all selected -particles, but in a different way. The particles must have a per-atom -mass attribute, as defined by the :doc:`atom_style ` -command. If the atom has a radius attribute (see :doc:`atom_style -sphere `) and its radius is non-zero, its mass is set from -the density and particle volume for 3d systems (the input density is -assumed to be in mass/distance\^3 units). For 2d, the default is for -LAMMPS to model particles with a radius attribute as spheres. However, -if the *density/disc* keyword is used, then they can be modeled as 2d -discs (circles). Their mass is set from the density and particle area -(the input density is assumed to be in mass/distance\^2 units). - -If the atom has a shape attribute (see :doc:`atom_style ellipsoid -`) and its 3 shape parameters are non-zero, then its mass is -set from the density and particle volume (the input density is assumed -to be in mass/distance\^3 units). The *density/disc* keyword has no -effect; it does not (yet) treat 3d ellipsoids as 2d ellipses. - -If the atom has a length attribute (see :doc:`atom_style line -`) and its length is non-zero, then its mass is set from the -density and line segment length (the input density is assumed to be in -mass/distance units). If the atom has an area attribute (see -:doc:`atom_style tri `) and its area is non-zero, then its -mass is set from the density and triangle area (the input density is -assumed to be in mass/distance\^2 units). - -If none of these cases are valid, then the mass is set to the density -value directly (the input density is assumed to be in mass units). - -Keyword *temperature* sets the temperature of a finite-size particle. -Currently, only the GRANULAR package supports this attribute. The -temperature must be added using an instance of -:doc:`fix property/atom ` The values for the -temperature must be positive. - -Keyword *volume* sets the volume of all selected particles. Currently, -only the :doc:`atom_style peri ` command defines particles -with a volume attribute. Note that this command does not adjust the -particle mass. - -Keyword *image* sets which image of the simulation box the atom is -considered to be in. An image of 0 means it is inside the box as -defined. A value of 2 means add 2 box lengths to get the true value. A -value of -1 means subtract 1 box length to get the true value. LAMMPS -updates these flags as atoms cross periodic boundaries during the -simulation. The flags can be output with atom snapshots via the -:doc:`dump ` command. If a value of NULL is specified for any of -nx,ny,nz, then the current image value for that dimension is unchanged. -For non-periodic dimensions only a value of 0 can be specified. This -command can be useful after a system has been equilibrated and atoms -have diffused one or more box lengths in various directions. This -command can then reset the image values for atoms so that they are -effectively inside the simulation box, e.g if a diffusion coefficient is -about to be measured via the :doc:`compute msd ` command. -Care should be taken not to reset the image flags of two atoms in a bond -to the same value if the bond straddles a periodic boundary (rather they -should be different by +/- 1). This will not affect the dynamics of a -simulation, but may mess up analysis of the trajectories if a LAMMPS -diagnostic or your own analysis relies on the image flags to unwrap a -molecule which straddles the periodic box. - -Keywords *bond*, *angle*, *dihedral*, and *improper*, set the bond -type (angle type, etc) of all bonds (angles, etc) of selected atoms to -the specified value. The value can be a numeric type from 1 to -nbondtypes (nangletypes, etc). Or it can be a type label (bond type -label, angle type label, etc). See the :doc:`Howto type labels +Keyword *type* sets the atom type for all selected atoms. A specified +value can be either a numeric atom type or an atom type label. When +using a numeric type, the specified value must be from 1 to ntypes, +where ntypes was set by the :doc:`create_box ` command or +the *atom types* field in the header of the data file read by the +:doc:`read_data ` command. When using a type label it must +have been defined previously. See the :doc:`Howto type labels ` doc page for the allowed syntax of type labels -and a general discussion of how type labels can be used. All atoms in -a particular bond (angle, etc) must be selected atoms in order for the -change to be made. The value of nbondtypes (nangletypes, etc) was set -by the *bond types* (\ *angle types*, etc) field in the header of the -data file read by the :doc:`read_data ` command. These -keywords do not allow use of an atom-style variable. +and a general discussion of how type labels can be used. -Keywords *rheo/rho* and *rheo/status* set the density and the status of -rheo particles. In particular, one can only set the phase in the status -as described by the :doc:`RHEO howto page `. +Keyword *type/fraction* sets the atom type for a fraction of the +selected atoms. The actual number of atoms changed is not guaranteed +to be exactly the specified fraction (0 <= *fraction* <= 1), but +should be statistically close. Random numbers are used in such a way +that a particular atom is changed or not changed, regardless of how +many processors are being used. This keyword does NOT allow use of an +atom-style variable. -Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat capacity, -and density of smoothed particle hydrodynamics (SPH) particles. See -`this PDF guide `_ to using SPH in LAMMPS. +Keywords *type/ratio* and *type/subset* also set the atom type for a +fraction of the selected atoms. The actual number of atoms changed +will be exactly the requested number. For *type/ratio* the specified +fraction (0 <= *fraction* <= 1) determines the number. For +*type/subset*, the specified *Nsubset* is the number. An iterative +algorithm is used which ensures the correct number of atoms are +selected, in a perfectly random fashion. Which atoms are selected +will change with the number of processors used. These keywords do not +allow use of an atom-style variable. -.. note:: +Keyword *volume* sets the volume of all selected particles. +Currently, only the :doc:`atom_style peri ` command +defines particles with a volume attribute. Note that this command +does not adjust the particle mass. - Please note that the SPH PDF guide file has not been updated for - many years and thus does not reflect the current *syntax* of the - SPH package commands. For that please refer to the LAMMPS manual. +Keywords *vx*, *vy*, and *vz* set the velocities of all selected +atoms. -Keyword *smd/mass/density* sets the mass of all selected particles, but -it is only applicable to the Smooth Mach Dynamics package MACHDYN. It -assumes that the particle volume has already been correctly set and -calculates particle mass from the provided mass density value. - -Keyword *smd/contact/radius* only applies to simulations with the Smooth -Mach Dynamics package MACHDYN. Itsets an interaction radius for -computing short-range interactions, e.g. repulsive forces to prevent -different individual physical bodies from penetrating each other. Note -that the SPH smoothing kernel diameter used for computing long range, -nonlocal interactions, is set using the *diameter* keyword. - -Keyword *dpd/theta* sets the internal temperature of a DPD particle as -defined by the DPD-REACT package. If the specified value is a number it -must be >= 0.0. If the specified value is NULL, then the kinetic -temperature Tkin of each particle is computed as 3/2 k Tkin = KE = 1/2 m -v\^2 = 1/2 m (vx\*vx+vy\*vy+vz\*vz). Each particle's internal -temperature is set to Tkin. If the specified value is an atom-style -variable, then the variable is evaluated for each particle. If a value ->= 0.0, the internal temperature is set to that value. If it is < 0.0, -the computation of Tkin is performed and the internal temperature is set -to that value. - -Keywords *edpd/temp* and *edpd/cv* set the temperature and volumetric -heat capacity of an eDPD particle as defined by the DPD-MESO package. -Currently, only :doc:`atom_style edpd ` defines particles -with these attributes. The values for the temperature and heat capacity -must be positive. - -Keyword *cc* sets the chemical concentration of a tDPD particle for a -specified species as defined by the DPD-MESO package. Currently, only -:doc:`atom_style tdpd ` defines particles with this -attribute. An integer for "index" selects a chemical species (1 to -Nspecies) where Nspecies is set by the atom_style command. The value for -the chemical concentration must be >= 0.0. - -Keyword *epsilon* sets the dielectric constant of a particle, precisely -of the medium where the particle resides as defined by the DIELECTRIC -package. Currently, only :doc:`atom_style dielectric ` -defines particles with this attribute. The value for the dielectric -constant must be >= 0.0. Note that the set command with this keyword -will rescale the particle charge accordingly so that the real charge -(e.g., as read from a data file) stays intact. To change the real -charges, one needs to use the set command with the *charge* -keyword. Care must be taken to ensure that the real and scaled charges, -and dielectric constants are consistent. +Keywords *x*, *y*, *z* set the coordinates of all selected atoms. Keywords *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 ` 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]. +and *d2_name*, the column of the array to set must also be included +following the name in brackets: e.g. d2_xyz[2] or i2_mySpin[3]. Restrictions """""""""""" @@ -584,7 +647,7 @@ You cannot set an atom attribute (e.g. *mol* or *q* or *volume*\ ) if the :doc:`atom_style ` does not have that attribute. This command requires inter-processor communication to coordinate the -setting of bond types (angle types, etc). This means that your system +setting of bond types (angle types, etc). This means that the system must be ready to perform a simulation before using one of these keywords (force fields set, atom mass set, etc). This is not necessary for other keywords. @@ -599,7 +662,7 @@ Related commands """""""""""""""" :doc:`create_box `, :doc:`create_atoms `, -:doc:`read_data ` +:doc:`read_data `, :doc:`fix set ` Default """"""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 993956fce9..668ad913ba 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -3369,6 +3369,7 @@ Rmin RMS rmsd rnage +rnflag rng rNEMD ro diff --git a/src/fix_set.cpp b/src/fix_set.cpp new file mode 100644 index 0000000000..ac4400c1a9 --- /dev/null +++ b/src/fix_set.cpp @@ -0,0 +1,85 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_set.h" + +#include "atom.h" +#include "error.h" +#include "set.h" +#include "update.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{SETCOMMAND,FIXSET}; // also used in Set class + +/* ---------------------------------------------------------------------- */ + +FixSet::FixSet(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 8) error->all(FLERR, 1, "Illegal fix set command: need at least eight arguments"); + + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR, "Fix {} Nevery must be > 0", style); + + reneighbor = utils::inumeric(FLERR, arg[4], false, lmp); + if (reneighbor < 0 || reneighbor > 1) + error->all(FLERR, "Fix {} rnflag must be 0/1", style); + + // create instance of Set class + + set = new Set(lmp); + + // pass remaining args to Set class + // only keywords which use per-atom variables are currently allowed + // NOTE: could also allow when set style = region, + // since atoms may move in/out of regions + + set->process_args(FIXSET,narg-5,&arg[5]); +} + +/* ---------------------------------------------------------------------- */ + +FixSet::~FixSet() +{ + delete set; +} + +/* ---------------------------------------------------------------------- */ + +int FixSet::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- + use the Set instance to update per-atom properties + NOTE: could return count of updated atoms from Set for use as fix output +---------------------------------------------------------------------- */ + +void FixSet::end_of_step() +{ + // select which atoms to act on + + set->selection(atom->nlocal); + + // loop over list of actions to reset atom attributes + + set->invoke_actions(); + + // trigger reneighboring on next timestep if requested + + if (reneighbor) next_reneighbor = update->ntimestep + 1; +} diff --git a/src/fix_set.h b/src/fix_set.h new file mode 100644 index 0000000000..7328ebd144 --- /dev/null +++ b/src/fix_set.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(set,FixSet); +// clang-format on +#else + +#ifndef LMP_FIX_SET_H +#define LMP_FIX_SET_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixSet : public Fix { + public: + FixSet(class LAMMPS *, int, char **); + ~FixSet() override; + int setmask() override; + void end_of_step() override; + + private: + int reneighbor; + + class Set *set; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/set.cpp b/src/set.cpp index 7ae41e9246..f42dd175ab 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -42,739 +42,418 @@ using namespace LAMMPS_NS; using namespace MathConst; +enum{SETCOMMAND,FIXSET}; // also used in FixSet class + enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; -enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET, - MOLECULE,X,Y,Z,VX,VY,VZ,CHARGE,MASS,SHAPE,LENGTH,TRI, - DIPOLE,DIPOLE_RANDOM,SPIN_ATOM,SPIN_RANDOM,SPIN_ELECTRON,RADIUS_ELECTRON, - QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,TEMPERATURE, - DIAMETER,RADIUS_ATOM,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, - RHEO_STATUS,SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY, - SMD_CONTACT_RADIUS,DPDTHETA,EPSILON,IVEC,DVEC,IARRAY,DARRAY}; +enum{ANGLE,ANGMOM,BOND,CC,CHARGE,DENSITY,DIAMETER,DIHEDRAL,DIPOLE, + DIPOLE_RANDOM,DPD_THETA,EDPD_CV,EDPD_TEMP,EPSILON,IMAGE,IMPROPER,LENGTH, + MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,SHAPE, + SMD_CONTACT_RADIUS,SMD_MASS_DENSITY,SPH_CV,SPH_E,SPH_RHO, + SPIN_ATOM,SPIN_ATOM_RANDOM,SPIN_ELECTRON,TEMPERATURE,THETA,THETA_RANDOM, + TRI,TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,VOLUME,VX,VY,VZ,X,Y,Z, + IVEC,DVEC,IARRAY,DARRAY}; + +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +Set::Set(class LAMMPS *lmp) : + Command(lmp), id(nullptr), region(nullptr), actions(nullptr), invoke_choice(nullptr), + vec1(nullptr), vec2(nullptr), vec3(nullptr), vec4(nullptr), select(nullptr) +{ + maxselect = maxvariable = 0; +} + +/* ---------------------------------------------------------------------- */ + +Set::~Set() +{ + memory->sfree(actions); + memory->sfree(invoke_choice); + + memory->destroy(select); + memory->destroy(vec1); + memory->destroy(vec2); + memory->destroy(vec3); + memory->destroy(vec4); +} /* ---------------------------------------------------------------------- */ void Set::command(int narg, char **arg) { if (domain->box_exist == 0) - error->all(FLERR, Error::NOLASTLINE, "Set command before simulation box is defined" - + utils::errorurl(0)); + error->all(FLERR, Error::NOLASTLINE, + "Set command before simulation box is defined" + utils::errorurl(0)); if (atom->natoms == 0) error->all(FLERR, Error::NOLASTLINE, "Set command on system without atoms"); if (narg < 4) error->all(FLERR, 1, "Illegal set command: need at least four arguments"); - // style and ID info + process_args(SETCOMMAND, narg, arg); - if (strcmp(arg[0],"atom") == 0) style = ATOM_SELECT; - else if (strcmp(arg[0],"mol") == 0) style = MOL_SELECT; - else if (strcmp(arg[0],"type") == 0) style = TYPE_SELECT; - else if (strcmp(arg[0],"group") == 0) style = GROUP_SELECT; - else if (strcmp(arg[0],"region") == 0) style = REGION_SELECT; - else error->all(FLERR, Error::ARGZERO, "Unknown set command style: {}", arg[0]); + if (comm->me == 0) utils::logmesg(lmp, "Setting atom values ...\n"); + + // select which atoms to act on - id = utils::strdup(arg[1]); - select = nullptr; selection(atom->nlocal); - // loop over keyword/value pairs - // call appropriate routine to reset attributes + // loop over list of actions to reset atom attributes - if (comm->me == 0) utils::logmesg(lmp,"Setting atom values ...\n"); + invoke_actions(); - int allcount,origarg; + // print stats for each action + // for CC option, include species index + + bigint bcount, allcount; + + for (int i = 0; i < naction; i++) { + Action *action = &actions[i]; + int iarg = action->argindex; + + if (action->count_action < 0) { + bcount = action->count_select; + } else { + bcount = action->count_action; + } + MPI_Allreduce(&bcount, &allcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + + if (comm->me == 0) { + if (strcmp(arg[iarg], "cc") == 0) { + utils::logmesg(lmp, " {} settings made for {} index {}\n", allcount, arg[iarg], + arg[iarg + 1]); + } else { + utils::logmesg(lmp, " {} settings made for {}\n", allcount, arg[iarg]); + } + } + } +} + +/* ---------------------------------------------------------------------- + process args of set command or fix set command +------------------------------------------------------------------------- */ + +void Set::process_args(int caller_flag, int narg, char **arg) +{ + caller = caller_flag; + + // style and ID info + + id = utils::strdup(arg[1]); + + if (strcmp(arg[0], "atom") == 0) { + style = ATOM_SELECT; + if (atom->tag_enable == 0) error->all(FLERR, "Cannot use set atom with no atom IDs defined"); + utils::bounds(FLERR, id, 1, MAXTAGINT, nlobig, nhibig, error); + + } else if (strcmp(arg[0], "mol") == 0) { + style = MOL_SELECT; + if (atom->molecule_flag == 0) + error->all(FLERR, "Cannot use set mol with no molecule IDs defined"); + utils::bounds(FLERR, id, 1, MAXTAGINT, nlobig, nhibig, error); + + } else if (strcmp(arg[0], "type") == 0) { + style = TYPE_SELECT; + if (char *typestr = utils::expand_type(FLERR, id, Atom::ATOM, lmp)) { + delete[] id; + id = typestr; + } + utils::bounds(FLERR, id, 1, atom->ntypes, nlo, nhi, error); + + } else if (strcmp(arg[0], "group") == 0) { + style = GROUP_SELECT; + int igroup = group->find(id); + if (igroup == -1) error->all(FLERR, "Could not find set group ID {}", id); + groupbit = group->bitmask[igroup]; + + } else if (strcmp(arg[0], "region") == 0) { + style = REGION_SELECT; + region = domain->get_region_by_id(id); + if (!region) error->all(FLERR, "Set region {} does not exist", id); + + } else + error->all(FLERR, "Unknown set or fix set command style: {}", arg[0]); + + delete[] id; + + // loop over remaining keyword/value pairs to create list of actions + // one action = keyword/value pair + + naction = maxaction = 0; + actions = nullptr; + Action *action; int iarg = 2; while (iarg < narg) { - varflag = varflag1 = varflag2 = varflag3 = varflag4 = 0; - count = 0; - origarg = iarg; - if (strcmp(arg[iarg],"type") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set type", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1], 1); - else ivalue = utils::expand_type_int(FLERR, arg[iarg+1], Atom::ATOM, lmp); - set(TYPE); - iarg += 2; + // grow list of actions if needed - } else if (strcmp(arg[iarg],"type/fraction") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/fraction", error); - newtype = utils::expand_type_int(FLERR, arg[iarg+1], Atom::ATOM, lmp); - fraction = utils::numeric(FLERR, arg[iarg+2], false, lmp); - ivalue = utils::inumeric(FLERR, arg[iarg+3], false, lmp); - if (newtype <= 0 || newtype > atom->ntypes) - error->all(FLERR, iarg + 1, "Invalid type value {} in set type/fraction command", newtype); - if (fraction < 0.0 || fraction > 1.0) - error->all(FLERR, iarg + 2, "Invalid fraction value {} in set type/fraction command", - fraction); - if (ivalue <= 0) - error->all(FLERR, iarg + 3, "Invalid random number seed {} in set type/fraction command", - ivalue); - setrandom(TYPE_FRACTION); - iarg += 4; - - } else if (strcmp(arg[iarg],"type/ratio") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/ratio", error); - newtype = utils::expand_type_int(FLERR, arg[iarg+1], Atom::ATOM, lmp); - fraction = utils::numeric(FLERR, arg[iarg+2], false, lmp); - ivalue = utils::inumeric(FLERR, arg[iarg+3], false, lmp); - if (newtype <= 0 || newtype > atom->ntypes) - error->all(FLERR, iarg + 1, "Invalid type value {} in set type/ratio command", newtype); - if (fraction < 0.0 || fraction > 1.0) - error->all(FLERR, iarg + 2, "Invalid fraction value {} in set type/ratio command", - fraction); - if (ivalue <= 0) - error->all(FLERR, iarg + 3, "Invalid random number seed {} in set type/ratio command", - ivalue); - setrandom(TYPE_RATIO); - iarg += 4; - - } else if (strcmp(arg[iarg],"type/subset") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/subset", error); - newtype = utils::expand_type_int(FLERR, arg[iarg+1], Atom::ATOM, lmp); - nsubset = utils::bnumeric(FLERR, arg[iarg+2], false, lmp); - ivalue = utils::inumeric(FLERR, arg[iarg+3], false, lmp); - if (newtype <= 0 || newtype > atom->ntypes) - error->all(FLERR, iarg + 1, "Invalid type value {} in set type/subset command", newtype); - if (nsubset < 0) - error->all(FLERR, iarg + 2, "Invalid subset size {} in set type/subset command", nsubset); - if (ivalue <= 0) - error->all(FLERR, iarg + 3, "Invalid random number seed {} in set type/subset command", - ivalue); - setrandom(TYPE_SUBSET); - iarg += 4; - - } else if (strcmp(arg[iarg],"mol") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set mol", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->molecule_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(MOLECULE); - iarg += 2; - - } else if (strcmp(arg[iarg],"x") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set x", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(X); - iarg += 2; - - } else if (strcmp(arg[iarg],"y") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set y", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(Y); - iarg += 2; - - } else if (strcmp(arg[iarg],"z") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set z", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(Z); - iarg += 2; - - } else if (strcmp(arg[iarg],"vx") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vx", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(VX); - iarg += 2; - - } else if (strcmp(arg[iarg],"vy") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vy", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(VY); - iarg += 2; - - } else if (strcmp(arg[iarg],"vz") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vz", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - set(VZ); - iarg += 2; - - } else if (strcmp(arg[iarg],"charge") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set charge", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->q_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(CHARGE); - iarg += 2; - - } else if (strcmp(arg[iarg],"mass") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set mass", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->rmass_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(MASS); - iarg += 2; - - } else if (strcmp(arg[iarg],"shape") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set shape", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (!atom->ellipsoid_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SHAPE); - iarg += 4; - - } else if (strcmp(arg[iarg],"length") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set length", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->line_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(LENGTH); - iarg += 2; - - } else if (strcmp(arg[iarg],"tri") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set tri", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->tri_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(TRI); - iarg += 2; - - } else if (strcmp(arg[iarg],"dipole") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set dipole", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (!atom->mu_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(DIPOLE); - iarg += 4; - - } else if (strcmp(arg[iarg],"dipole/random") == 0) { - if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "set dipole/random", error); - ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - dvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (!atom->mu_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0) - error->all(FLERR, iarg + 1, "Invalid random number seed in set command"); - if (dvalue <= 0.0) - error->all(FLERR, iarg + 2, "Invalid dipole length in set command"); - setrandom(DIPOLE_RANDOM); - iarg += 3; - - } else if ((strcmp(arg[iarg],"spin") == 0) || (strcmp(arg[iarg],"spin/atom") == 0)) { - if ((strcmp(arg[iarg],"spin") == 0) && (comm->me == 0)) - error->warning(FLERR, "Set attribute spin is deprecated. Please use spin/atom instead."); - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "set spin/atom", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4); - else zvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - if ((xvalue == 0.0) && (yvalue == 0.0) && (zvalue == 0.0)) - error->all(FLERR, Error::NOPOINTER, "At least one spin vector component must be non-zero"); - if (!atom->sp_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (dvalue <= 0.0) - error->all(FLERR, iarg + 1, "Invalid spin magnitude {} in set {} command", dvalue, - arg[iarg]); - set(SPIN_ATOM); - iarg += 5; - - } else if ((strcmp(arg[iarg],"spin/random") == 0) || - (strcmp(arg[iarg],"spin/atom/random") == 0)) { - if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "set spin/atom/random", error); - ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - dvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if ((strcmp(arg[iarg],"spin/random") == 0) && (comm->me == 0)) - error->warning(FLERR, "Set attribute spin/random is deprecated. " - "Please use spin/atom/random instead."); - if (!atom->sp_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0) - error->all(FLERR, iarg + 1, "Invalid random number seed {} in set {} command", ivalue, - arg[iarg]); - if (dvalue <= 0.0) - error->all(FLERR, iarg + 2, "Invalid spin magnitude {} in set {} command", dvalue, - arg[iarg]); - setrandom(SPIN_RANDOM); - iarg += 3; - - } else if (strcmp(arg[iarg],"radius/electron") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set radius/electron", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->eradius_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(RADIUS_ELECTRON); - iarg += 2; - - } else if (strcmp(arg[iarg],"spin/electron") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set spin/electron", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->spin_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SPIN_ELECTRON); - iarg += 2; - - } else if (strcmp(arg[iarg],"quat") == 0) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "set quat", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4); - else wvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(QUAT); - iarg += 5; - - } else if (strcmp(arg[iarg],"quat/random") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set quat/random", error); - ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0) - error->all(FLERR, iarg + 1, "Invalid random number seed in set command"); - setrandom(QUAT_RANDOM); - iarg += 2; - - } else if (strcmp(arg[iarg],"theta") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set theta", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = DEG2RAD * utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->line_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(THETA); - iarg += 2; - - } else if (strcmp(arg[iarg],"theta/random") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set theta/random", error); - ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->line_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0) - error->all(FLERR, iarg + 1, "Invalid random number seed in set command"); - set(THETA_RANDOM); - iarg += 2; - - } else if (strcmp(arg[iarg],"angmom") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set angmom", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (!atom->angmom_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(ANGMOM); - iarg += 4; - - } else if (strcmp(arg[iarg],"omega") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set omega", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (!atom->omega_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(OMEGA); - iarg += 4; - - } else if (strcmp(arg[iarg],"radius/atom") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set radius/atom", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->radius_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(RADIUS_ATOM); - iarg += 2; - - } else if (strcmp(arg[iarg],"diameter") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set diameter", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->radius_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(DIAMETER); - iarg += 2; - - } else if (strcmp(arg[iarg],"density") == 0 || - (strcmp(arg[iarg],"density/disc") == 0)) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set density", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->rmass_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (dvalue <= 0.0) error->all(FLERR, iarg + 1, "Invalid density in set command"); - discflag = 0; - if (strcmp(arg[iarg],"density/disc") == 0) { - discflag = 1; - if (domain->dimension != 2) - error->all(FLERR, Error::NOLASTLINE, "Density/disc option requires 2d simulation"); - } - set(DENSITY); - iarg += 2; - - } else if (strcmp(arg[iarg],"temperature") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->temperature_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(TEMPERATURE); - iarg += 2; - - } else if (strcmp(arg[iarg],"volume") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set volume", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->vfrac_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (dvalue <= 0.0) error->all(FLERR, iarg + 1, "Invalid volume in set command"); - set(VOLUME); - iarg += 2; - - } else if (strcmp(arg[iarg],"image") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set image", error); - ximageflag = yimageflag = zimageflag = 0; - if (strcmp(arg[iarg+1],"NULL") != 0) { - ximageflag = 1; - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else ximage = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - } - if (strcmp(arg[iarg+2],"NULL") != 0) { - yimageflag = 1; - if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2); - else yimage = utils::inumeric(FLERR,arg[iarg+2],false,lmp); - } - if (strcmp(arg[iarg+3],"NULL") != 0) { - zimageflag = 1; - if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3); - else zimage = utils::inumeric(FLERR,arg[iarg+3],false,lmp); - } - if (ximageflag && ximage && !domain->xperiodic) - error->all(FLERR, iarg + 1,"Cannot set non-zero image flag for non-periodic dimension"); - if (yimageflag && yimage && !domain->yperiodic) - error->all(FLERR, iarg + 2, "Cannot set non-zero image flag for non-periodic dimension"); - if (zimageflag && zimage && !domain->zperiodic) - error->all(FLERR, iarg + 3, "Cannot set non-zero image flag for non-periodic dimension"); - set(IMAGE); - iarg += 4; - - } else if (strcmp(arg[iarg],"bond") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set bond", error); - ivalue = utils::expand_type_int(FLERR, arg[iarg+1], Atom::BOND, lmp); - if (atom->avec->bonds_allow == 0) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0 || ivalue > atom->nbondtypes) - error->all(FLERR, iarg + 1, "Invalid value {} in set bond command", ivalue); - topology(BOND); - iarg += 2; - - } else if (strcmp(arg[iarg],"angle") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set angle", error); - ivalue = utils::expand_type_int(FLERR, arg[iarg+1], Atom::ANGLE, lmp); - if (atom->avec->angles_allow == 0) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0 || ivalue > atom->nangletypes) - error->all(FLERR, iarg + 1, "Invalid value {} in set angle command", ivalue); - topology(ANGLE); - iarg += 2; - - } else if (strcmp(arg[iarg],"dihedral") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set dihedral", error); - ivalue = utils::expand_type_int(FLERR, arg[iarg+1], Atom::DIHEDRAL, lmp); - if (atom->avec->dihedrals_allow == 0) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0 || ivalue > atom->ndihedraltypes) - error->all(FLERR, iarg + 1, "Invalid value {} in set dihedral command", ivalue); - topology(DIHEDRAL); - iarg += 2; - - } else if (strcmp(arg[iarg],"improper") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set improper", error); - ivalue = utils::expand_type_int(FLERR, arg[iarg+1], Atom::IMPROPER, lmp); - if (atom->avec->impropers_allow == 0) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - if (ivalue <= 0 || ivalue > atom->nimpropertypes) - error->all(FLERR, iarg + 1, "Invalid value {} in set improper command", ivalue); - topology(IMPROPER); - iarg += 2; - - } else if (strcmp(arg[iarg],"rheo/rho") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/rho", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->rho_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SPH_RHO); - iarg += 2; - - } else if (strcmp(arg[iarg],"rheo/status") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->rheo_status_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(RHEO_STATUS); - iarg += 2; - - } else if (strcmp(arg[iarg],"sph/e") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/e", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->esph_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SPH_E); - iarg += 2; - - } else if (strcmp(arg[iarg],"sph/cv") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/cv", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->cv_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SPH_CV); - iarg += 2; - - } else if (strcmp(arg[iarg],"sph/rho") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/rho", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->rho_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SPH_RHO); - iarg += 2; - - } else if (strcmp(arg[iarg],"edpd/temp") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set edpd/temp", error); - if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; - else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else { - dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (dvalue < 0.0) - error->all(FLERR, iarg + 1, "Invalid value {} in set edpd/temp command", dvalue); - } - if (!atom->edpd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(EDPD_TEMP); - iarg += 2; - - } else if (strcmp(arg[iarg],"edpd/cv") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set edpd/cv", error); - if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; - else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else { - dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (dvalue < 0.0) - error->all(FLERR, iarg + 1, "Invalid value {} in set edpd/cv command", dvalue); - } - if (!atom->edpd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(EDPD_CV); - iarg += 2; - - } else if (strcmp(arg[iarg],"cc") == 0) { - if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "set cc", error); - if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; - else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else { - cc_index = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - dvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (cc_index < 1) error->all(FLERR,"Illegal set command"); - error->all(FLERR, iarg + 1, "Invalid index value {} in set cc command", cc_index); - } - if (!atom->tdpd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(CC); - iarg += 3; - - } else if (strcmp(arg[iarg],"smd/mass/density") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set smd/mass/density", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->smd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SMD_MASS_DENSITY); - iarg += 2; - - } else if (strcmp(arg[iarg],"smd/contact/radius") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set smd/contact/radius", error); - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->smd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(SMD_CONTACT_RADIUS); - iarg += 2; - - } else if (strcmp(arg[iarg],"dpd/theta") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set dpd/theta", error); - if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; - else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else { - dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (dvalue < 0.0) - error->all(FLERR, iarg + 1, "Invalid value {} in set dpd/theta command", dvalue); - } - if (!atom->dpd_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(DPDTHETA); - iarg += 2; - - } else if (strcmp(arg[iarg],"epsilon") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set epsilon", error); - if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; - else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else { - dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (dvalue < 0.0) - error->all(FLERR, iarg + 1, "Invalid value {} in set epsilon command", dvalue); - } - if (!atom->dielectric_flag) - error->all(FLERR, iarg, "Cannot set attribute {} for atom style {}", arg[iarg], - atom->get_style()); - set(EPSILON); - iarg += 2; - - } else { - - // set custom per-atom vector or array or error out - - int flag,cols; - ArgInfo argi(arg[iarg],ArgInfo::DNAME|ArgInfo::INAME); - const char *pname = argi.get_name(); - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set", error); - index_custom = atom->find_custom(argi.get_name(),flag,cols); - if (index_custom < 0) - error->all(FLERR, iarg, "Set keyword or custom property {} does not exist", pname); - - switch (argi.get_type()) { - - case ArgInfo::INAME: - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (flag != 0) - error->all(FLERR, iarg + 1, "Set command custom property {} is not integer", pname); - - if (argi.get_dim() == 0) { - if (cols > 0) - error->all(FLERR, iarg, "Set command custom integer property {} is not a vector", - pname); - set(IVEC); - } else if (argi.get_dim() == 1) { - if (cols == 0) - error->all(FLERR, iarg, "Set command custom integer property {} is not an array", - pname); - icol_custom = argi.get_index1(); - if (icol_custom <= 0 || icol_custom > cols) - error->all(FLERR, iarg, "Set command per-atom custom integer array {} is accessed " - "out-of-range{}", pname, utils::errorurl(20)); - set(IARRAY); - } else error->all(FLERR,"Illegal set command"); - break; - - case ArgInfo::DNAME: - if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); - else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (flag != 1) error->all(FLERR,"Custom property {} is not floating-point",argi.get_name()); - - if (argi.get_dim() == 0) { - if (cols > 0) - error->all(FLERR, iarg, "Set command custom double property {} is not a vector", pname); - set(DVEC); - } else if (argi.get_dim() == 1) { - if (cols == 0) - error->all(FLERR, iarg, "Set command custom double property {} is not an array", pname); - icol_custom = argi.get_index1(); - if (icol_custom <= 0 || icol_custom > cols) - error->all(FLERR, iarg, "Set command per-atom custom double array {} is accessed " - "out-of-range{}", pname, utils::errorurl(20)); - set(DARRAY); - } else error->all(FLERR,"Illegal set command"); - break; - - default: - error->all(FLERR,"Illegal set command"); - break; - } - iarg += 2; + if (naction == maxaction) { + maxaction += DELTA; + actions = (Action *) memory->srealloc(actions,maxaction*sizeof(Action),"set:actions"); + invoke_choice = (FnPtrPack *) + memory->srealloc(invoke_choice,maxaction*sizeof(FnPtrPack),"set:invoke_choice"); } - // statistics - // for CC option, include species index + action = &actions[naction]; + action->argindex = iarg; + action->varflag = 0; + action->varflag1 = action->varflag2 = action->varflag3 = action->varflag4 = 0; - MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world); + if (strcmp(arg[iarg],"angle") == 0) { + action->keyword = ANGLE; + process_angle(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_angle; + } else if (strcmp(arg[iarg],"angmom") == 0) { + action->keyword = ANGMOM; + process_angmom(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_angmom; + } else if (strcmp(arg[iarg],"bond") == 0) { + action->keyword = BOND; + process_bond(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_bond; + } else if (strcmp(arg[iarg],"cc") == 0) { + action->keyword = CC; + process_cc(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_cc; + } else if (strcmp(arg[iarg],"charge") == 0) { + action->keyword = CHARGE; + process_charge(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_charge; + } else if (strcmp(arg[iarg],"density") == 0 ||(strcmp(arg[iarg],"density/disc") == 0)) { + action->keyword = DENSITY; + process_density(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_density; + } else if (strcmp(arg[iarg],"diameter") == 0) { + action->keyword = DIAMETER; + process_diameter(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_diameter; + } else if (strcmp(arg[iarg],"dihedral") == 0) { + action->keyword = DIHEDRAL; + process_dihedral(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_dihedral; + } else if (strcmp(arg[iarg],"dipole") == 0) { + action->keyword = DIPOLE; + process_dipole(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_dipole; + } else if (strcmp(arg[iarg],"dipole/random") == 0) { + action->keyword = DIPOLE_RANDOM; + process_dipole_random(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_dipole_random; + } else if (strcmp(arg[iarg],"dpd/theta") == 0) { + action->keyword = DPD_THETA; + process_dpd_theta(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_dpd_theta; + } else if (strcmp(arg[iarg],"edpd/cv") == 0) { + action->keyword = EDPD_CV; + process_edpd_cv(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_edpd_cv; + } else if (strcmp(arg[iarg],"edpd/temp") == 0) { + action->keyword = EDPD_TEMP; + process_edpd_temp(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_edpd_temp; + } else if (strcmp(arg[iarg],"epsilon") == 0) { + action->keyword = EPSILON; + process_epsilon(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_epsilon; + } else if (strcmp(arg[iarg],"image") == 0) { + action->keyword = IMAGE; + process_image(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_image; + } else if (strcmp(arg[iarg],"improper") == 0) { + action->keyword = IMPROPER; + process_improper(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_improper; + } else if (strcmp(arg[iarg],"length") == 0) { + action->keyword = LENGTH; + process_length(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_length; + } else if (strcmp(arg[iarg],"mass") == 0) { + action->keyword = MASS; + process_mass(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_mass; + } else if (strcmp(arg[iarg],"mol") == 0) { + action->keyword = MOLECULE; + process_mol(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_mol; + } else if (strcmp(arg[iarg],"omega") == 0) { + action->keyword = OMEGA; + process_omega(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_omega; + } else if (strcmp(arg[iarg],"quat") == 0) { + action->keyword = QUAT; + process_quat(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_quat; + } else if (strcmp(arg[iarg],"quat/random") == 0) { + action->keyword = QUAT_RANDOM; + process_quat_random(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_quat_random; + } else if (strcmp(arg[iarg],"radius/electron") == 0) { + action->keyword = RADIUS_ELECTRON; + process_radius_election(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_radius_election; + } else if (strcmp(arg[iarg],"shape") == 0) { + action->keyword = SHAPE; + process_shape(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_shape; + } else if (strcmp(arg[iarg],"smd/contact/radius") == 0) { + action->keyword = SMD_CONTACT_RADIUS; + process_smd_contact_radius(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_smd_contact_radius; + } else if (strcmp(arg[iarg],"smd/mass/density") == 0) { + action->keyword = SMD_MASS_DENSITY; + process_smd_mass_density(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_smd_mass_density; + } else if (strcmp(arg[iarg],"sph/cv") == 0) { + action->keyword = SPH_CV; + process_sph_cv(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_sph_cv; + } else if (strcmp(arg[iarg],"sph/e") == 0) { + action->keyword = SPH_E; + process_sph_e(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_sph_e; + } else if (strcmp(arg[iarg],"sph/rho") == 0) { + action->keyword = SPH_RHO; + process_sph_rho(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_sph_rho; + } else if ((strcmp(arg[iarg],"spin/atom") == 0) || (strcmp(arg[iarg],"spin") == 0)) { + action->keyword = SPIN_ATOM; + process_spin_atom(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_spin_atom; + } else if ((strcmp(arg[iarg],"spin/atom/random") == 0) || (strcmp(arg[iarg],"spin/random") == 0)) { + action->keyword = SPIN_ATOM_RANDOM; + process_spin_atom_random(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_spin_atom_random; + } else if (strcmp(arg[iarg],"spin/electron") == 0) { + action->keyword = SPIN_ELECTRON; + process_spin_electron(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_spin_electron; + } else if (strcmp(arg[iarg],"temperature") == 0) { + action->keyword = TEMPERATURE; + process_temperature(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_temperature; + } else if (strcmp(arg[iarg],"theta") == 0) { + action->keyword = THETA; + process_theta(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_theta; + } else if (strcmp(arg[iarg],"theta/random") == 0) { + action->keyword = THETA_RANDOM; + process_theta_random(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_theta_random; + } else if (strcmp(arg[iarg],"tri") == 0) { + action->keyword = TRI; + process_tri(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_tri; + } else if (strcmp(arg[iarg],"type") == 0) { + action->keyword = TYPE; + process_type(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_type; + } else if (strcmp(arg[iarg],"type/fraction") == 0) { + action->keyword = TYPE_FRACTION; + process_type_fraction(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_type_fraction; + } else if (strcmp(arg[iarg],"type/ratio") == 0) { + action->keyword = TYPE_RATIO; + process_type_ratio(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_type_ratio; + } else if (strcmp(arg[iarg],"type/subset") == 0) { + action->keyword = TYPE_SUBSET; + process_type_subset(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_type_subset; + } else if (strcmp(arg[iarg],"volume") == 0) { + action->keyword = VOLUME; + process_volume(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_volume; + } else if (strcmp(arg[iarg],"vx") == 0) { + action->keyword = VX; + process_vx(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_vx; + } else if (strcmp(arg[iarg],"vy") == 0) { + action->keyword = VY; + process_vy(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_vy; + } else if (strcmp(arg[iarg],"vz") == 0) { + action->keyword = VZ; + process_vz(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_vz; + } else if (strcmp(arg[iarg],"x") == 0) { + action->keyword = X; + process_x(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_x; + } else if (strcmp(arg[iarg],"y") == 0) { + action->keyword = Y; + process_y(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_y; + } else if (strcmp(arg[iarg],"z") == 0) { + action->keyword = Z; + process_z(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_z; - if (comm->me == 0) { - if (strcmp(arg[origarg],"cc") == 0) - utils::logmesg(lmp," {} settings made for {} index {}\n", - allcount,arg[origarg],arg[origarg+1]); - else - utils::logmesg(lmp," {} settings made for {}\n", - allcount,arg[origarg]); + } else if (utils::strmatch(arg[iarg],"^[id]2?_")) { + process_custom(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_custom; + + // unrecognized keyword + + } else error->all(FLERR,"Unrecognized set or fix set command keyword {}",arg[iarg]); + } + + // varflag = 1 if any action uses a per-atom variable + // varflag1234 = 1 if any action uses a specific per-atom variable + + varflag = 0; + varflag1 = varflag2 = varflag3 = varflag4 = 0; + for (int i = 0; i < naction; i++) { + action = &actions[i]; + if (action->varflag) varflag = 1; + if (action->varflag1) varflag1 = 1; + if (action->varflag2) varflag2 = 1; + if (action->varflag3) varflag3 = 1; + if (action->varflag4) varflag4 = 1; + } + + // error if any action of fix set command does not use a per-atom variable + // b/c fix set is then effectivly a no-op + + if (caller == FIXSET) { + for (int i = 0; i < naction; i++) { + action = &actions[i]; + if (!action->varflag) + error->all(FLERR,"Fix set command keyword {} does not invoke a per-atom variable", + arg[action->argindex]); } } - // free local memory + // warn if a keyword sets properties for atoms in rigid bodies + // assume no conflict for properties not in this list of cases - delete[] id; - delete[] select; + for (int i = 0; i < naction; i++) { + switch (actions[i].keyword) { + case X: + case Y: + case Z: + case MOLECULE: + case MASS: + case ANGMOM: + case SHAPE: + case DIAMETER: + case DENSITY: + case TEMPERATURE: + case QUAT: + case IMAGE: + if (modify->check_rigid_list_overlap(select)) + error->warning(FLERR,"Setting a property of atoms in rigid bodies " + "that has no effect unless rigid bodies are re-initialized"); + break; + default: + break; + } + } + + // in future, could return index to next arg to process + // i.e. if fix set command appends its own options + + // return iarg; } /* ---------------------------------------------------------------------- @@ -784,468 +463,143 @@ void Set::command(int narg, char **arg) void Set::selection(int n) { - delete[] select; - select = new int[n]; - int nlo,nhi; + // reallocate select vector if needed + // this method could be called many times by fix set command + + if (n > maxselect) { + memory->destroy(select); + memory->create(select,n,"set:select"); + maxselect = n; + } if (style == ATOM_SELECT) { - if (atom->tag_enable == 0) - error->all(FLERR, Error::NOLASTLINE, "Cannot use set atom with no atom IDs defined"); - bigint nlobig,nhibig; - utils::bounds(FLERR,id,1,MAXTAGINT,nlobig,nhibig,error); - tagint *tag = atom->tag; for (int i = 0; i < n; i++) if (tag[i] >= nlobig && tag[i] <= nhibig) select[i] = 1; else select[i] = 0; } else if (style == MOL_SELECT) { - if (atom->molecule_flag == 0) - error->all(FLERR, Error::NOLASTLINE, "Cannot use set mol with no molecule IDs defined"); - bigint nlobig,nhibig; - utils::bounds(FLERR,id,1,MAXTAGINT,nlobig,nhibig,error); - tagint *molecule = atom->molecule; for (int i = 0; i < n; i++) if (molecule[i] >= nlobig && molecule[i] <= nhibig) select[i] = 1; else select[i] = 0; } else if (style == TYPE_SELECT) { - utils::bounds_typelabel(FLERR,id,1,atom->ntypes,nlo,nhi,lmp,Atom::ATOM); - int *type = atom->type; for (int i = 0; i < n; i++) if (type[i] >= nlo && type[i] <= nhi) select[i] = 1; else select[i] = 0; } else if (style == GROUP_SELECT) { - int igroup = group->find(id); - if (igroup == -1) error->all(FLERR, Error::NOLASTLINE, "Could not find set group ID {}", id); - int groupbit = group->bitmask[igroup]; - int *mask = atom->mask; for (int i = 0; i < n; i++) if (mask[i] & groupbit) select[i] = 1; else select[i] = 0; } else if (style == REGION_SELECT) { - auto region = domain->get_region_by_id(id); - if (!region) error->all(FLERR, Error::NOLASTLINE, "Set region {} does not exist", id); region->prematch(); - double **x = atom->x; for (int i = 0; i < n; i++) - if (region->match(x[i][0],x[i][1],x[i][2])) - select[i] = 1; + if (region->match(x[i][0],x[i][1],x[i][2])) select[i] = 1; else select[i] = 0; } + + // count_select = count of selected owned atoms + + count_select = 0; + for (int i = 0; i < n; i++) + if (select[i]) count_select++; } /* ---------------------------------------------------------------------- - set owned atom properties directly - either scalar or per-atom values from atom-style variable(s) + loop over list of actions + perform each action on all selected atoms via call to invoke_choice() method ------------------------------------------------------------------------- */ -void Set::set(int keyword) +void Set::invoke_actions() { - // evaluate atom-style variable(s) if necessary + // reallocate per-atom variable storage if needed - vec1 = vec2 = vec3 = vec4 = nullptr; - - if (varflag) { - int nlocal = atom->nlocal; + if (varflag && atom->nlocal > maxvariable) { + maxvariable = atom->nlocal; if (varflag1) { - memory->create(vec1,nlocal,"set:vec1"); - input->variable->compute_atom(ivar1,0,vec1,1,0); + memory->destroy(vec1); + memory->create(vec1,maxvariable,"set:var1"); } if (varflag2) { - memory->create(vec2,nlocal,"set:vec2"); - input->variable->compute_atom(ivar2,0,vec2,1,0); + memory->destroy(vec2); + memory->create(vec2,maxvariable,"set:var2"); } if (varflag3) { - memory->create(vec3,nlocal,"set:vec3"); - input->variable->compute_atom(ivar3,0,vec3,1,0); + memory->destroy(vec3); + memory->create(vec3,maxvariable,"set:var3"); } if (varflag4) { - memory->create(vec4,nlocal,"set:vec4"); - input->variable->compute_atom(ivar4,0,vec4,1,0); + memory->destroy(vec4); + memory->create(vec4,maxvariable,"set:var4"); } } - // check if properties of atoms in rigid bodies are updated - // that are cached as per-body data. - switch (keyword) { - case X: - case Y: - case Z: - case MOLECULE: - case MASS: - case ANGMOM: - case SHAPE: - case DIAMETER: - case DENSITY: - case TEMPERATURE: - case QUAT: - case IMAGE: - if (modify->check_rigid_list_overlap(select)) - error->warning(FLERR,"Changing a property of atoms in rigid bodies " - "that has no effect unless rigid bodies are rebuild"); - break; - default: // assume no conflict for all other properties - break; - } + // loop over actions - // loop over selected atoms + for (int i = 0; i < naction; i++) { - auto avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); - auto avec_line = dynamic_cast(atom->style_match("line")); - auto avec_tri = dynamic_cast(atom->style_match("tri")); - auto avec_body = dynamic_cast(atom->style_match("body")); + Action *action = &actions[i]; - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - if (!select[i]) continue; + // use count_action to optionally override count_select + // if stays -1, count_select is used by caller + // if overwritten by an invoke method, count_action is used + // only a handful of invoke methods tally their own count - // overwrite dvalue, ivalue, xyzw value if variables defined - // else the input script scalar value remains in place + count_action = -1; - if (varflag) { - if (varflag1) { - dvalue = xvalue = vec1[i]; - ivalue = static_cast (dvalue); + // evaluate atom-style variable(s) if necessary + + if (action->varflag) { + if (action->varflag1) { + input->variable->compute_atom(action->ivar1,0,vec1,1,0); } - if (varflag2) yvalue = vec2[i]; - if (varflag3) zvalue = vec3[i]; - if (varflag4) wvalue = vec4[i]; - } - - // set values in per-atom arrays - // error check here in case atom-style variables generated bogus value - - if (keyword == TYPE) { - if (ivalue <= 0 || ivalue > atom->ntypes) - error->one(FLERR, Error::NOLASTLINE, "Invalid value {} in set type command", ivalue); - atom->type[i] = ivalue; - } - else if (keyword == MOLECULE) atom->molecule[i] = ivalue; - else if (keyword == X) atom->x[i][0] = dvalue; - else if (keyword == Y) atom->x[i][1] = dvalue; - else if (keyword == Z) atom->x[i][2] = dvalue; - else if (keyword == VX) atom->v[i][0] = dvalue; - else if (keyword == VY) atom->v[i][1] = dvalue; - else if (keyword == VZ) atom->v[i][2] = dvalue; - else if (keyword == CHARGE) { - atom->q[i] = dvalue; - // ensure that scaled charges are consistent the new charge value - if (atom->epsilon) atom->q_scaled[i] = dvalue / atom->epsilon[i]; - } else if (keyword == MASS) { - if (dvalue <= 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid mass {} in set command", dvalue); - atom->rmass[i] = dvalue; - } - else if (keyword == DIAMETER) { - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid diameter {} in set command", dvalue); - atom->radius[i] = 0.5 * dvalue; - } - else if (keyword == VOLUME) { - if (dvalue <= 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid volume {} in set command", dvalue); - atom->vfrac[i] = dvalue; - } - - else if (keyword == RHEO_STATUS) { - if (ivalue != 0 && ivalue != 1) - error->one(FLERR, Error::NOLASTLINE, "Invalid value {} in set command for rheo/status", - ivalue); - atom->rheo_status[i] = ivalue; - } - - else if (keyword == SPH_E) atom->esph[i] = dvalue; - else if (keyword == SPH_CV) atom->cv[i] = dvalue; - else if (keyword == SPH_RHO) atom->rho[i] = dvalue; - - else if (keyword == EDPD_TEMP) atom->edpd_temp[i] = dvalue; - else if (keyword == EDPD_CV) atom->edpd_cv[i] = dvalue; - else if (keyword == CC) atom->cc[i][cc_index-1] = dvalue; - - else if (keyword == SMD_MASS_DENSITY) { - // set mass from volume and supplied mass density - atom->rmass[i] = atom->vfrac[i] * dvalue; - } - else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue; - - else if (keyword == DPDTHETA) { - if (dvalue >= 0.0) atom->dpdTheta[i] = dvalue; - else { - double onemass; - if (atom->rmass) onemass = atom->rmass[i]; - else onemass = atom->mass[atom->type[i]]; - double vx = atom->v[i][0]; - double vy = atom->v[i][1]; - double vz = atom->v[i][2]; - double tfactor = force->mvv2e / (domain->dimension * force->boltz); - atom->dpdTheta[i] = tfactor * onemass * (vx*vx + vy*vy + vz*vz); + if (action->varflag2) { + input->variable->compute_atom(action->ivar2,0,vec2,1,0); + } + if (action->varflag3) { + input->variable->compute_atom(action->ivar3,0,vec3,1,0); + } + if (action->varflag4) { + input->variable->compute_atom(action->ivar4,0,vec4,1,0); } } - // set shape of ellipsoidal particle + // invoke the action to reset per-atom or per-topology values - else if (keyword == SHAPE) { - if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid ellipsoid shape {} {} {} in set command", - xvalue, yvalue, zvalue); - if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) { - if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid ellipsoid shape {} {} {} in set command", - xvalue, yvalue, zvalue); - } - avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue); - } + (this->*invoke_choice[i])(action); - // set length of line particle - - else if (keyword == LENGTH) { - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid length {} in set command", dvalue); - avec_line->set_length(i,dvalue); - } - - // set corners of tri particle - - else if (keyword == TRI) { - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid length {} in set command", dvalue); - avec_tri->set_equilateral(i,dvalue); - } - - // set rmass via density - // if radius > 0.0, treat as sphere or disc - // if shape > 0.0, treat as ellipsoid (or ellipse, when uncomment below) - // if length > 0.0, treat as line - // if area > 0.0, treat as tri - // else set rmass to density directly - - else if (keyword == DENSITY) { - if (dvalue <= 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid density {} in set command", dvalue); - if (atom->radius_flag && atom->radius[i] > 0.0) - if (discflag) - atom->rmass[i] = MY_PI*atom->radius[i]*atom->radius[i] * dvalue; - else - atom->rmass[i] = 4.0*MY_PI/3.0 * - atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue; - else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) { - double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape; - // enable 2d ellipse (versus 3d ellipsoid) when time integration - // options (fix nve/asphere, fix nh/asphere) are also implemented - // if (discflag) - // atom->rmass[i] = MY_PI*shape[0]*shape[1] * dvalue; - // else - atom->rmass[i] = 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue; - } else if (atom->line_flag && atom->line[i] >= 0) { - double length = avec_line->bonus[atom->line[i]].length; - atom->rmass[i] = length * dvalue; - } else if (atom->tri_flag && atom->tri[i] >= 0) { - double *c1 = avec_tri->bonus[atom->tri[i]].c1; - double *c2 = avec_tri->bonus[atom->tri[i]].c2; - double *c3 = avec_tri->bonus[atom->tri[i]].c3; - double c2mc1[3],c3mc1[3]; - MathExtra::sub3(c2,c1,c2mc1); - MathExtra::sub3(c3,c1,c3mc1); - double norm[3]; - MathExtra::cross3(c2mc1,c3mc1,norm); - double area = 0.5 * MathExtra::len3(norm); - atom->rmass[i] = area * dvalue; - } else atom->rmass[i] = dvalue; - } - - // set dipole moment - - else if (keyword == DIPOLE) { - double **mu = atom->mu; - mu[i][0] = xvalue; - mu[i][1] = yvalue; - mu[i][2] = zvalue; - mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + - mu[i][2]*mu[i][2]); - } - - // set magnetic moments - - else if (keyword == SPIN_ATOM) { - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Incorrect value for atomic spin magnitude: {}", - dvalue); - double **sp = atom->sp; - double inorm = 1.0/sqrt(xvalue*xvalue+yvalue*yvalue+zvalue*zvalue); - sp[i][0] = inorm*xvalue; - sp[i][1] = inorm*yvalue; - sp[i][2] = inorm*zvalue; - sp[i][3] = dvalue; - } - - // set electron radius - - else if (keyword == RADIUS_ELECTRON) { - atom->eradius[i] = dvalue; - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Incorrect value for electron radius: {}", dvalue); - } - - // set electron spin - - else if (keyword == SPIN_ELECTRON) { - if ((dvalue == -1) || (dvalue == 1) || (dvalue == 0) || (dvalue == 2) || (dvalue == 3)) - atom->spin[i] = (int)dvalue; - else - error->one(FLERR, Error::NOLASTLINE, "Incorrect value for electron spin: {}", dvalue); - } - - // set quaternion orientation of ellipsoid or tri or body particle or sphere/bpm - // enforce quat rotation vector in z dir for 2d systems - - else if (keyword == QUAT) { - double *quat = nullptr; - double **quat2 = nullptr; - if (avec_ellipsoid && atom->ellipsoid[i] >= 0) - quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; - else if (avec_tri && atom->tri[i] >= 0) - quat = avec_tri->bonus[atom->tri[i]].quat; - else if (avec_body && atom->body[i] >= 0) - quat = avec_body->bonus[atom->body[i]].quat; - else if (atom->quat_flag) - quat2 = atom->quat; - else - error->one(FLERR, Error::NOLASTLINE, "Cannot set quaternion for atom that has none"); - if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0)) - error->one(FLERR, Error::NOLASTLINE, - "Cannot set quaternion with xy components for 2d system"); - - const double theta2 = MY_PI2 * wvalue/180.0; - const double sintheta2 = sin(theta2); - double temp[4]; - temp[0] = cos(theta2); - temp[1] = xvalue * sintheta2; - temp[2] = yvalue * sintheta2; - temp[3] = zvalue * sintheta2; - MathExtra::qnormalize(temp); - if (atom->quat_flag) { - quat2[i][0] = temp[0]; - quat2[i][1] = temp[1]; - quat2[i][2] = temp[2]; - quat2[i][3] = temp[3]; - } else { - quat[0] = temp[0]; - quat[1] = temp[1]; - quat[2] = temp[2]; - quat[3] = temp[3]; - } - } - - // set theta of line particle - - else if (keyword == THETA) { - if (atom->line[i] < 0) - error->one(FLERR, Error::NOLASTLINE, "Cannot set theta for atom that is not a line"); - avec_line->bonus[atom->line[i]].theta = dvalue; - } - - // set angmom or omega of particle - - else if (keyword == ANGMOM) { - atom->angmom[i][0] = xvalue; - atom->angmom[i][1] = yvalue; - atom->angmom[i][2] = zvalue; - } - - else if (keyword == OMEGA) { - atom->omega[i][0] = xvalue; - atom->omega[i][1] = yvalue; - atom->omega[i][2] = zvalue; - } - - // set temperature of particle - - else if (keyword == TEMPERATURE) { - if (dvalue < 0.0) - error->one(FLERR, Error::NOLASTLINE, "Invalid temperature {} in set command", dvalue); - atom->temperature[i] = dvalue; - } - - // reset any or all of 3 image flags - - else if (keyword == IMAGE) { - int xbox = (atom->image[i] & IMGMASK) - IMGMAX; - int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX; - int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX; - if (varflag1) ximage = static_cast(xvalue); - if (varflag2) yimage = static_cast(yvalue); - if (varflag3) zimage = static_cast(zvalue); - if (ximageflag) xbox = ximage; - if (yimageflag) ybox = yimage; - if (zimageflag) zbox = zimage; - atom->image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) | - (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | - (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); - } - - // set the local dielectric constant - - else if (keyword == EPSILON) { - if (dvalue >= 0.0) { - - // assign the new local dielectric constant - // update both the scaled charge value - - atom->epsilon[i] = dvalue; - atom->q_scaled[i] = atom->q[i] / dvalue; - } - } - - // set value for custom property vector or array - - else if (keyword == IVEC) { - atom->ivector[index_custom][i] = ivalue; - } - - else if (keyword == DVEC) { - 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++; + action->count_select = count_select; + action->count_action = count_action; } +} - // update bonus data numbers +/* ---------------------------------------------------------------------- */ - if (keyword == SHAPE) { - bigint nlocal_bonus = avec_ellipsoid->nlocal_bonus; - MPI_Allreduce(&nlocal_bonus,&atom->nellipsoids,1, - MPI_LMP_BIGINT,MPI_SUM,world); +void Set::varparse(const char *name, int m, Action *action) +{ + int ivar = input->variable->find(name+2); + if (ivar < 0) + error->all(FLERR,"Variable name {} for set command does not exist", name); + if (!input->variable->atomstyle(ivar)) + error->all(FLERR,"Variable {} for set command is invalid style", name); + + action->varflag = 1; + + if (m == 1) { + action->varflag1 = 1; action->ivar1 = ivar; + } else if (m == 2) { + action->varflag2 = 1; action->ivar2 = ivar; + } else if (m == 3) { + action->varflag3 = 1; action->ivar3 = ivar; + } else if (m == 4) { + action->varflag4 = 1; action->ivar4 = ivar; } - if (keyword == LENGTH) { - bigint nlocal_bonus = avec_line->nlocal_bonus; - MPI_Allreduce(&nlocal_bonus,&atom->nlines,1,MPI_LMP_BIGINT,MPI_SUM,world); - } - if (keyword == TRI) { - bigint nlocal_bonus = avec_tri->nlocal_bonus; - MPI_Allreduce(&nlocal_bonus,&atom->ntris,1,MPI_LMP_BIGINT,MPI_SUM,world); - } - - // clear up per-atom memory if allocated - - memory->destroy(vec1); - memory->destroy(vec2); - memory->destroy(vec3); - memory->destroy(vec4); } /* ---------------------------------------------------------------------- @@ -1254,7 +608,7 @@ void Set::set(int keyword) make atom result independent of what proc owns it ------------------------------------------------------------------------- */ -void Set::setrandom(int keyword) +void Set::setrandom(int keyword, Action *action) { int i; @@ -1264,7 +618,10 @@ void Set::setrandom(int keyword) auto avec_body = dynamic_cast(atom->style_match("body")); double **x = atom->x; - int seed = ivalue; + + // seed is always set to ivalue1 in process() methods + + int seed = action->ivalue1; auto ranpark = new RanPark(lmp,1); auto ranmars = new RanMars(lmp,seed + comm->me); @@ -1273,6 +630,9 @@ void Set::setrandom(int keyword) if (keyword == TYPE_FRACTION) { int nlocal = atom->nlocal; + double fraction = action->dvalue1; + int newtype = action->ivalue2; + int count = 0; for (i = 0; i < nlocal; i++) if (select[i]) { @@ -1282,44 +642,45 @@ void Set::setrandom(int keyword) count++; } + count_action = count; + // set exact count of atom types to newtype // for TYPE_RATIO, exact = fraction out of total eligible // for TYPE_SUBSET, exact = nsubset out of total eligible } else if (keyword == TYPE_RATIO || keyword == TYPE_SUBSET) { int nlocal = atom->nlocal; + int newtype = action->ivalue2; - // count = number of eligible atoms I own + // convert specified fraction to nsubset of all selected atoms - count = 0; - for (i = 0; i < nlocal; i++) - if (select[i]) count++; - - // convert specified fraction to nsubset - - bigint bcount = count; + bigint bcount = count_select; bigint allcount; MPI_Allreduce(&bcount,&allcount,1,MPI_LMP_BIGINT,MPI_SUM,world); + bigint nsubset; if (keyword == TYPE_RATIO) { + double fraction = action->dvalue1; nsubset = static_cast (fraction * allcount); } else if (keyword == TYPE_SUBSET) { + nsubset = action->bvalue1; if (nsubset > allcount) - error->all(FLERR, Error::NOLASTLINE, "Set type/subset value exceeds eligible atoms"); + error->all(FLERR,"Set type/subset value exceeds eligible atoms"); } // make selection - int *flag = memory->create(flag,count,"set:flag"); - int *work = memory->create(work,count,"set:work"); + int *flag = memory->create(flag,count_select,"set:flag"); + int *work = memory->create(work,count_select,"set:work"); - ranmars->select_subset(nsubset,count,flag,work); + ranmars->select_subset(nsubset,count_select,flag,work); // change types of selected atoms // flag vector from select_subset() is only for eligible atoms - count = 0; + int count = 0; int eligible = 0; + for (i = 0; i < nlocal; i++) { if (!select[i]) continue; if (flag[eligible]) { @@ -1329,6 +690,8 @@ void Set::setrandom(int keyword) eligible++; } + count_action = count; + // clean up memory->destroy(flag); @@ -1340,7 +703,7 @@ void Set::setrandom(int keyword) } else if (keyword == DIPOLE_RANDOM) { double **mu = atom->mu; int nlocal = atom->nlocal; - + double dmag = action->dvalue1; double msq,scale; if (domain->dimension == 3) { @@ -1351,12 +714,11 @@ void Set::setrandom(int keyword) mu[i][1] = ranpark->uniform() - 0.5; mu[i][2] = ranpark->uniform() - 0.5; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; - scale = dvalue/sqrt(msq); + scale = dmag/sqrt(msq); mu[i][0] *= scale; mu[i][1] *= scale; mu[i][2] *= scale; - mu[i][3] = dvalue; - count++; + mu[i][3] = dmag; } } else { @@ -1367,22 +729,20 @@ void Set::setrandom(int keyword) mu[i][1] = ranpark->uniform() - 0.5; mu[i][2] = 0.0; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; - scale = dvalue/sqrt(msq); + scale = dmag/sqrt(msq); mu[i][0] *= scale; mu[i][1] *= scale; - mu[i][3] = dvalue; - count++; + mu[i][3] = dmag; } } - // set spin moments to random orientations in 3d or 2d // spin length is fixed to unity - } else if (keyword == SPIN_RANDOM) { + } else if (keyword == SPIN_ATOM_RANDOM) { double **sp = atom->sp; int nlocal = atom->nlocal; - + double dlen = action->dvalue1; double sp_sq,scale; if (domain->dimension == 3) { @@ -1397,8 +757,7 @@ void Set::setrandom(int keyword) sp[i][0] *= scale; sp[i][1] *= scale; sp[i][2] *= scale; - sp[i][3] = dvalue; - count++; + sp[i][3] = dlen; } } else { @@ -1412,32 +771,35 @@ void Set::setrandom(int keyword) scale = 1.0/sqrt(sp_sq); sp[i][0] *= scale; sp[i][1] *= scale; - sp[i][3] = dvalue; - count++; + sp[i][3] = dlen; } } // set quaternions to random orientations in 3d and 2d } else if (keyword == QUAT_RANDOM) { + int *ellipsoid = atom->ellipsoid; + int *tri = atom->tri; + int *body = atom->body; + double **quat = atom->quat; int nlocal = atom->nlocal; - double *quat; - double **quat2; + int quat_flag = atom->quat_flag; + double *quat_one; if (domain->dimension == 3) { double s,t1,t2,theta1,theta2; for (i = 0; i < nlocal; i++) if (select[i]) { - if (avec_ellipsoid && atom->ellipsoid[i] >= 0) - quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; - else if (avec_tri && atom->tri[i] >= 0) - quat = avec_tri->bonus[atom->tri[i]].quat; - else if (avec_body && atom->body[i] >= 0) - quat = avec_body->bonus[atom->body[i]].quat; - else if (atom->quat_flag) - quat2 = atom->quat; + if (avec_ellipsoid && ellipsoid[i] >= 0) + quat_one = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (avec_tri && tri[i] >= 0) + quat_one = avec_tri->bonus[tri[i]].quat; + else if (avec_body && body[i] >= 0) + quat_one = avec_body->bonus[body[i]].quat; + else if (quat_flag) + quat_one = quat[i]; else - error->one(FLERR, Error::NOLASTLINE, "Cannot set quaternion for atom that has none"); + error->one(FLERR,"Cannot set quaternion for atom that has none"); ranpark->reset(seed,x[i]); s = ranpark->uniform(); @@ -1445,61 +807,46 @@ void Set::setrandom(int keyword) t2 = sqrt(s); theta1 = 2.0*MY_PI*ranpark->uniform(); theta2 = 2.0*MY_PI*ranpark->uniform(); - if (atom->quat_flag) { - quat2[i][0] = cos(theta2)*t2; - quat2[i][1] = sin(theta1)*t1; - quat2[i][2] = cos(theta1)*t1; - quat2[i][3] = sin(theta2)*t2; - } else { - quat[0] = cos(theta2)*t2; - quat[1] = sin(theta1)*t1; - quat[2] = cos(theta1)*t1; - quat[3] = sin(theta2)*t2; - } - count++; + quat_one[0] = cos(theta2)*t2; + quat_one[1] = sin(theta1)*t1; + quat_one[2] = cos(theta1)*t1; + quat_one[3] = sin(theta2)*t2; } } else { double theta2; for (i = 0; i < nlocal; i++) if (select[i]) { - if (avec_ellipsoid && atom->ellipsoid[i] >= 0) - quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; - else if (avec_body && atom->body[i] >= 0) - quat = avec_body->bonus[atom->body[i]].quat; - else if (atom->quat_flag) - quat2 = atom->quat; + if (avec_ellipsoid && ellipsoid[i] >= 0) + quat_one = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (avec_body && body[i] >= 0) + quat_one = avec_body->bonus[body[i]].quat; + else if (quat_flag) + quat_one = quat[i]; else - error->one(FLERR, Error::NOLASTLINE, "Cannot set quaternion for atom that has none"); + error->one(FLERR,"Cannot set quaternion for atom that has none"); ranpark->reset(seed,x[i]); theta2 = MY_PI*ranpark->uniform(); - if (atom->quat_flag) { - quat2[i][0] = cos(theta2); - quat2[i][1] = 0.0; - quat2[i][2] = 0.0; - quat2[i][3] = sin(theta2); - } else { - quat[0] = cos(theta2); - quat[1] = 0.0; - quat[2] = 0.0; - quat[3] = sin(theta2); - } - count++; + quat_one[0] = cos(theta2); + quat_one[1] = 0.0; + quat_one[2] = 0.0; + quat_one[3] = sin(theta2); } } // set theta to random orientation in 2d } else if (keyword == THETA_RANDOM) { + int *line = atom->line; int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) { if (select[i]) { - if (atom->line[i] < 0) - error->one(FLERR, Error::NOLASTLINE, "Cannot set theta for atom that is not a line"); + if (line[i] < 0) + error->one(FLERR,"Cannot set theta for atom that is not a line"); ranpark->reset(seed,x[i]); - avec_line->bonus[atom->line[i]].theta = MY_2PI*ranpark->uniform(); - count++; + avec_line->bonus[line[i]].theta = MY_2PI*ranpark->uniform(); } } } @@ -1510,14 +857,14 @@ void Set::setrandom(int keyword) /* ---------------------------------------------------------------------- */ -void Set::topology(int keyword) +void Set::topology(int keyword, Action *action) { int m,atom1,atom2,atom3,atom4; // error check if (atom->molecular == Atom::TEMPLATE) - error->all(FLERR, Error::NOLASTLINE, "Cannot set bond topology types for atom style template"); + error->all(FLERR,"Cannot set bond topology types for atom style template"); // border swap to acquire ghost atom info // enforce PBC before in case atoms are outside box @@ -1539,18 +886,23 @@ void Set::topology(int keyword) selection(atom->nlocal + atom->nghost); + int count = 0; + // for BOND, each of 2 atoms must be in group if (keyword == BOND) { + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + tagint **bond_atom = atom->bond_atom; int nlocal = atom->nlocal; + + int itype = action->ivalue1; for (int i = 0; i < nlocal; i++) - for (m = 0; m < atom->num_bond[i]; m++) { - atom1 = atom->map(atom->bond_atom[i][m]); - if (atom1 == -1) - error->one(FLERR, Error::NOLASTLINE, "Bond atom missing in set command" - + utils::errorurl(5)); + for (m = 0; m < num_bond[i]; m++) { + atom1 = atom->map(bond_atom[i][m]); + if (atom1 == -1) error->one(FLERR,"Bond atom missing in set command"); if (select[i] && select[atom1]) { - atom->bond_type[i][m] = ivalue; + bond_type[i][m] = itype; count++; } } @@ -1559,17 +911,23 @@ void Set::topology(int keyword) // for ANGLE, each of 3 atoms must be in group if (keyword == ANGLE) { + int *num_angle = atom->num_angle; + int **angle_type = atom->angle_type; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; int nlocal = atom->nlocal; + + int itype = action->ivalue1; for (int i = 0; i < nlocal; i++) - for (m = 0; m < atom->num_angle[i]; m++) { - atom1 = atom->map(atom->angle_atom1[i][m]); - atom2 = atom->map(atom->angle_atom2[i][m]); - atom3 = atom->map(atom->angle_atom3[i][m]); + for (m = 0; m < num_angle[i]; m++) { + atom1 = atom->map(angle_atom1[i][m]); + atom2 = atom->map(angle_atom2[i][m]); + atom3 = atom->map(angle_atom3[i][m]); if (atom1 == -1 || atom2 == -1 || atom3 == -1) - error->one(FLERR, Error::NOLASTLINE, "Angle atom missing in set command" - + utils::errorurl(5)); + error->one(FLERR,"Angle atom missing in set command"); if (select[atom1] && select[atom2] && select[atom3]) { - atom->angle_type[i][m] = ivalue; + angle_type[i][m] = itype; count++; } } @@ -1578,18 +936,25 @@ void Set::topology(int keyword) // for DIHEDRAL, each of 4 atoms must be in group if (keyword == DIHEDRAL) { + int *num_dihedral = atom->num_dihedral; + int **dihedral_type = atom->dihedral_type; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; int nlocal = atom->nlocal; + + int itype = action->ivalue1; for (int i = 0; i < nlocal; i++) - for (m = 0; m < atom->num_dihedral[i]; m++) { - atom1 = atom->map(atom->dihedral_atom1[i][m]); - atom2 = atom->map(atom->dihedral_atom2[i][m]); - atom3 = atom->map(atom->dihedral_atom3[i][m]); - atom4 = atom->map(atom->dihedral_atom4[i][m]); + for (m = 0; m < num_dihedral[i]; m++) { + atom1 = atom->map(dihedral_atom1[i][m]); + atom2 = atom->map(dihedral_atom2[i][m]); + atom3 = atom->map(dihedral_atom3[i][m]); + atom4 = atom->map(dihedral_atom4[i][m]); if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) - error->one(FLERR, Error::NOLASTLINE, "Dihedral atom missing in set command" - + utils::errorurl(5)); + error->one(FLERR,"Dihedral atom missing in set command"); if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) { - atom->dihedral_type[i][m] = ivalue; + dihedral_type[i][m] = itype; count++; } } @@ -1598,43 +963,1976 @@ void Set::topology(int keyword) // for IMPROPER, each of 4 atoms must be in group if (keyword == IMPROPER) { + int *num_improper = atom->num_improper; + int **improper_type = atom->improper_type; + tagint **improper_atom1 = atom->improper_atom1; + tagint **improper_atom2 = atom->improper_atom2; + tagint **improper_atom3 = atom->improper_atom3; + tagint **improper_atom4 = atom->improper_atom4; int nlocal = atom->nlocal; + + int itype = action->ivalue1; for (int i = 0; i < nlocal; i++) - for (m = 0; m < atom->num_improper[i]; m++) { - atom1 = atom->map(atom->improper_atom1[i][m]); - atom2 = atom->map(atom->improper_atom2[i][m]); - atom3 = atom->map(atom->improper_atom3[i][m]); - atom4 = atom->map(atom->improper_atom4[i][m]); + for (m = 0; m < num_improper[i]; m++) { + atom1 = atom->map(improper_atom1[i][m]); + atom2 = atom->map(improper_atom2[i][m]); + atom3 = atom->map(improper_atom3[i][m]); + atom4 = atom->map(improper_atom4[i][m]); if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) - error->one(FLERR, Error::NOLASTLINE, "Improper atom missing in set command" - + utils::errorurl(5)); + error->one(FLERR,"Improper atom missing in set command"); if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) { - atom->improper_type[i][m] = ivalue; + improper_type[i][m] = itype; count++; } } } + + // set count_action for all topology actions + + count_action = count; +} + +// ---------------------------------------------------------------------- +// pairs of process/invoke methods for each keyword +// process method reads args, stores parameters in Action instance +// invoke method resets atoms properties using Action instance +// separate two operations so can be called by either set or fix set command +// ---------------------------------------------------------------------- + +/* ---------------------------------------------------------------------- */ + +void Set::process_angle(int &iarg, int narg, char **arg, Action *action) +{ + if (atom->avec->angles_allow == 0) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set angle", error); + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::ANGLE,lmp); + action->ivalue1 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue1 <= 0 || action->ivalue1 > atom->nangletypes) + error->all(FLERR,"Invalid angle type in set command"); + iarg += 2; +} + +void Set::invoke_angle(Action *action) +{ + topology(ANGLE,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_angmom(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->angmom_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set angmom", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + iarg += 4; +} + +void Set::invoke_angmom(Action *action) +{ + int nlocal = atom->nlocal; + double **angmom = atom->angmom; + + int varflag = action->varflag; + double xvalue,yvalue,zvalue; + if (!action->varflag1) xvalue = action->dvalue1; + if (!action->varflag2) yvalue = action->dvalue2; + if (!action->varflag3) zvalue = action->dvalue3; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) xvalue = vec1[i]; + if (action->varflag1) yvalue = vec2[i]; + if (action->varflag1) zvalue = vec3[i]; + } + + angmom[i][0] = xvalue; + angmom[i][1] = yvalue; + angmom[i][2] = zvalue; + } } /* ---------------------------------------------------------------------- */ -void Set::varparse(const char *name, int m) +void Set::process_bond(int &iarg, int narg, char **arg, Action *action) { - varflag = 1; - int ivar = input->variable->find(name+2); + if (atom->avec->bonds_allow == 0) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set bond", error); - if (ivar < 0) - error->all(FLERR, Error::NOLASTLINE, "Variable name {} for set command does not exist", name); - if (!input->variable->atomstyle(ivar)) - error->all(FLERR, Error::NOLASTLINE, "Variable {} for set command is invalid style", name); + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::BOND,lmp); + action->ivalue1 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue1 <= 0 || action->ivalue1 > atom->nbondtypes) + error->all(FLERR,"Invalid bond type in set command"); - if (m == 1) { - varflag1 = 1; ivar1 = ivar; - } else if (m == 2) { - varflag2 = 1; ivar2 = ivar; - } else if (m == 3) { - varflag3 = 1; ivar3 = ivar; - } else if (m == 4) { - varflag4 = 1; ivar4 = ivar; + iarg += 2; +} + +void Set::invoke_bond(Action *action) +{ + topology(BOND,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_cc(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->tdpd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "set cc", error); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 < 1) error->all(FLERR,"Invalid cc index in set command"); + + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (action->dvalue1 < 0.0) error->all(FLERR,"Invalid cc value in set command"); + } + + iarg += 3; +} + +void Set::invoke_cc(Action *action) +{ + int nlocal = atom->nlocal; + double **cc = atom->cc; + + int cc_index = action->ivalue1 - 1; + // NOTE: need to check if cc_index exceeds cc array allocation + + int varflag = action->varflag; + double ccvalue; + if (!action->varflag1) ccvalue = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + ccvalue = vec1[i]; + if (ccvalue < 0.0) error->all(FLERR,"Invalid cc value in set command"); + } + + cc[i][cc_index] = ccvalue; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_charge(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->q_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set charge", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_charge(Action *action) +{ + int nlocal = atom->nlocal; + double *q = atom->q; + double *q_scaled = atom->q_scaled; + double *epsilon = atom->epsilon; + + int varflag = action->varflag; + double qvalue; + if (!action->varflag1) qvalue = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) qvalue = vec1[i]; + q[i] = qvalue; + + // ensure scaled charges are consistent with new charge value + + if (epsilon) q_scaled[i] = qvalue / epsilon[i]; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_density(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->rmass_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set density", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) error->all(FLERR,"Invalid density in set command"); + } + + action->ivalue1 = 0; + if (strcmp(arg[iarg],"density/disc") == 0) { + action->ivalue1 = 1; + if (domain->dimension != 2) error->all(FLERR,"Set density/disc requires 2d simulation"); + } + + iarg += 2; +} + +void Set::invoke_density(Action *action) +{ + int nlocal = atom->nlocal; + double *rmass = atom->rmass; + double *radius = atom->radius; + int *ellipsoid = atom->ellipsoid; + int *line = atom->line; + int *tri = atom->tri; + + int radius_flag = atom->radius_flag; + int ellipsoid_flag = atom->ellipsoid_flag; + int line_flag = atom->line_flag; + int tri_flag = atom->tri_flag; + + auto avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); + auto avec_line = dynamic_cast(atom->style_match("line")); + auto avec_tri = dynamic_cast(atom->style_match("tri")); + + int varflag = action->varflag; + double density; + if (!action->varflag1) density = action->dvalue1; + int discflag = action->ivalue1; + + // set rmass via density + // if radius > 0.0, treat as sphere or disc + // if shape > 0.0, treat as ellipsoid (or ellipse, when uncomment below) + // if length > 0.0, treat as line + // if area > 0.0, treat as tri + // else set rmass to density directly + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + density = vec1[i]; + if (density <= 0.0) error->one(FLERR,"Invalid density in set command"); + } + + if (radius_flag && radius[i] > 0.0) + if (discflag) rmass[i] = MY_PI*radius[i]*radius[i] * density; + else rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density; + + else if (ellipsoid_flag && ellipsoid[i] >= 0) { + double *shape = avec_ellipsoid->bonus[ellipsoid[i]].shape; + // could enable 2d ellipse (versus 3d ellipsoid) when time integration + // options (fix nve/asphere, fix nh/asphere) are also implemented + // if (discflag) + // atom->rmass[i] = MY_PI*shape[0]*shape[1] * dvalue; + // else + rmass[i] = 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2] * density; + + } else if (line_flag && line[i] >= 0) { + double length = avec_line->bonus[line[i]].length; + rmass[i] = length * density; + + } else if (tri_flag && tri[i] >= 0) { + double *c1 = avec_tri->bonus[tri[i]].c1; + double *c2 = avec_tri->bonus[tri[i]].c2; + double *c3 = avec_tri->bonus[tri[i]].c3; + double c2mc1[3],c3mc1[3]; + MathExtra::sub3(c2,c1,c2mc1); + MathExtra::sub3(c3,c1,c3mc1); + double norm[3]; + MathExtra::cross3(c2mc1,c3mc1,norm); + double area = 0.5 * MathExtra::len3(norm); + rmass[i] = area * density; + + } else rmass[i] = density; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_diameter(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->radius_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set diameter", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid diameter in set command"); + } + + iarg += 2; +} + +void Set::invoke_diameter(Action *action) +{ + int nlocal = atom->nlocal; + double *radius = atom->radius; + + int varflag = action->varflag; + double diam; + if (!action->varflag1) diam = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + diam = vec1[i]; + if (diam < 0.0) error->one(FLERR,"Invalid diameter in set command"); + } + + radius[i] = 0.5 * diam; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_dihedral(int &iarg, int narg, char **arg, Action *action) +{ + if (atom->avec->dihedrals_allow == 0) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set dihedral", error); + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::DIHEDRAL,lmp); + action->ivalue1 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue1 <= 0 || action->ivalue1 > atom->ndihedraltypes) + error->all(FLERR,"Invalid dihedral type in set command"); + + iarg += 2; +} + +void Set::invoke_dihedral(Action *action) +{ + topology(DIHEDRAL,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_dipole(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->mu_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set dipole", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + iarg += 4; +} + +void Set::invoke_dipole(Action *action) +{ + int nlocal = atom->nlocal; + double **mu = atom->mu; + + int varflag = action->varflag; + double xvalue,yvalue,zvalue; + if (!action->varflag1) xvalue = action->dvalue1; + if (!action->varflag2) yvalue = action->dvalue2; + if (!action->varflag3) zvalue = action->dvalue3; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) xvalue = vec1[i]; + if (action->varflag1) yvalue = vec2[i]; + if (action->varflag1) zvalue = vec3[i]; + } + + mu[i][0] = xvalue; + mu[i][1] = yvalue; + mu[i][2] = zvalue; + mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]); + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_dipole_random(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->mu_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "set dipole/random", error); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + action->dvalue1 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (action->ivalue1 <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + if (action->dvalue1 <= 0.0) + error->all(FLERR,"Invalid dipole length in set command"); + + iarg += 3; +} + +void Set::invoke_dipole_random(Action *action) +{ + setrandom(DIPOLE_RANDOM,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_dpd_theta(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->dpd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set dpd/theta", error); + + if (strcmp(arg[iarg+1],"NULL") == 0) action->dvalue1 = -1.0; + else if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->all(FLERR,"Invalid dpd/theta value in set command"); + } + + iarg += 2; +} + +void Set::invoke_dpd_theta(Action *action) +{ + int nlocal = atom->nlocal; + int *type = atom->type; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + double *dpdTheta = atom->dpdTheta; + + double tfactor = force->mvv2e / (domain->dimension * force->boltz); + double onemass; + double vx,vy,vz; + + int varflag = action->varflag; + double theta; + if (!action->varflag1) theta = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + theta = vec1[i]; + if (theta < 0.0) error->one(FLERR,"Invalid dpd/theta value in set command"); + } + + // if theta is negative, NULL was used, set dpdTheta to KE of particle + + if (theta >= 0.0) dpdTheta[i] = theta; + else { + if (rmass) onemass = rmass[i]; + else onemass = mass[type[i]]; + vx = v[i][0]; + vy = v[i][1]; + vz = v[i][2]; + dpdTheta[i] = tfactor * onemass * (vx*vx + vy*vy + vz*vz); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_edpd_cv(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->edpd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set edpd/cv", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->all(FLERR,"Invalid edpd/cv value in set command"); + } + + iarg += 2; +} + +void Set::invoke_edpd_cv(Action *action) +{ + int nlocal = atom->nlocal; + double *edpd_cv = atom->edpd_cv; + + int varflag = action->varflag; + double cv; + if (!action->varflag1) cv = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + cv = vec1[i]; + if (cv < 0.0) error->one(FLERR,"Invalid edpd/cv value in set command"); + } + + edpd_cv[i] = cv; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_edpd_temp(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->edpd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set edpd/temp", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->all(FLERR,"Invalid edpd/temp value in set command"); + } + iarg += 2; +} + +void Set::invoke_edpd_temp(Action *action) +{ + int nlocal = atom->nlocal; + double *edpd_temp = atom->edpd_temp; + + int varflag = action->varflag; + double temp; + if (!action->varflag1) temp = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + temp = vec1[i]; + if (temp < 0.0) error->one(FLERR,"Invalid edpd/temp value in set command"); + } + + edpd_temp[i] = temp; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_epsilon(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->dielectric_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set epsilon", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) error->all(FLERR,"Invalid epsilon in set command"); + } + + iarg += 2; +} + +void Set::invoke_epsilon(Action *action) +{ + int nlocal = atom->nlocal; + double *epsilon = atom->epsilon; + double *q = atom->q; + double *q_scaled = atom->q_scaled; + + int varflag = action->varflag; + double eps; + if (!action->varflag1) eps = action->dvalue1; + + // assign local dielectric constant + // also update scaled charge value + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + eps = vec1[i]; + if (eps <= 0.0) error->one(FLERR,"Invalid epsilon in set command"); + } + + epsilon[i] = eps; + q_scaled[i] = q[i] / eps; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_image(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set image", error); + + if (strcmp(arg[iarg+1],"NULL") == 0) action->ivalue4 = 0; + else { + action->ivalue4 = 1; + if (utils::strmatch(arg[iarg+1],"^v_")) { + if (!domain->xperiodic) + error->all(FLERR,"Cannot set variable image flag for non-periodic dimension"); + varparse(arg[iarg+1],1,action); + } else { + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 && !domain->xperiodic) + error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension"); + } + } + + if (strcmp(arg[iarg+2],"NULL") == 0) action->ivalue5 = 0; + else { + action->ivalue5 = 1; + if (utils::strmatch(arg[iarg+2],"^v_")) { + if (!domain->yperiodic) + error->all(FLERR,"Cannot set variable image flag for non-periodic dimension"); + varparse(arg[iarg+2],2,action); + } else { + action->ivalue2 = utils::inumeric(FLERR,arg[iarg+2],false,lmp); + if (action->ivalue2 && !domain->yperiodic) + error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension"); + } + } + + if (strcmp(arg[iarg+1],"NULL") == 0) action->ivalue6 = 0; + else { + action->ivalue6 = 1; + if (utils::strmatch(arg[iarg+3],"^v_")) { + if (!domain->zperiodic) + error->all(FLERR,"Cannot set variable image flag for non-periodic dimension"); + varparse(arg[iarg+3],3,action); + } else { + action->ivalue3 = utils::inumeric(FLERR,arg[iarg+3],false,lmp); + if (action->ivalue3 && !domain->zperiodic) + error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension"); + } + } + + iarg += 4; +} + +void Set::invoke_image(Action *action) +{ + int nlocal = atom->nlocal; + imageint *image = atom->image; + int xbox,ybox,zbox; + + int ximageflag = action->ivalue4; + int yimageflag = action->ivalue5; + int zimageflag = action->ivalue6; + + int varflag = action->varflag; + int ximage,yimage,zimage; + if (!action->varflag1) ximage = action->ivalue1; + if (!action->varflag2) yimage = action->ivalue2; + if (!action->varflag3) zimage = action->ivalue3; + + // reset any or all of 3 image flags + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) ximage = static_cast (vec1[i]); + if (action->varflag2) yimage = static_cast (vec2[i]); + if (action->varflag3) zimage = static_cast (vec3[i]); + } + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + if (ximageflag) xbox = ximage; + if (yimageflag) ybox = yimage; + if (zimageflag) zbox = zimage; + image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) | + (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | + (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_improper(int &iarg, int narg, char **arg, Action *action) +{ + if (atom->avec->impropers_allow == 0) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set improper", error); + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::IMPROPER,lmp); + action->ivalue1 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue1 <= 0 || action->ivalue1 > atom->nimpropertypes) + error->all(FLERR,"Invalid value in set command"); + + iarg += 2; +} + +void Set::invoke_improper(Action *action) +{ + topology(IMPROPER,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_length(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->line_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set length", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid length in set command"); + } + + iarg += 2; +} + +void Set::invoke_length(Action *action) +{ + int nlocal = atom->nlocal; + auto avec_line = dynamic_cast(atom->style_match("line")); + + int varflag = action->varflag; + double length; + if (!action->varflag1) length = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + length = vec1[i]; + if (length < 0.0) error->one(FLERR,"Invalid length in set command"); + } + + avec_line->set_length(i,length); + } + + // update global line count + + bigint nlocal_bonus = avec_line->nlocal_bonus; + MPI_Allreduce(&nlocal_bonus,&atom->nlines,1,MPI_LMP_BIGINT,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_mass(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->rmass_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set mass", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) error->one(FLERR,"Invalid mass in set command"); + } + + iarg += 2; +} + +void Set::invoke_mass(Action *action) +{ + int nlocal = atom->nlocal; + double *rmass = atom->rmass; + + int varflag = action->varflag; + double mass_one; + if (!action->varflag1) mass_one = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + mass_one = vec1[i]; + if (mass_one < 0.0) error->one(FLERR,"Invalid mass in set command"); + } + + rmass[i] = mass_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_mol(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->molecule_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set mol", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->tvalue1 = utils::tnumeric(FLERR,arg[iarg+1],false,lmp); + if (action->tvalue1 < 0) error->one(FLERR,"Invalid molecule ID in set command"); + } + + iarg += 2; +} + +void Set::invoke_mol(Action *action) +{ + int nlocal = atom->nlocal; + tagint *molecule = atom->molecule; + + int varflag = action->varflag; + tagint molID; + if (!action->varflag1) molID = action->tvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + molID = vec1[i]; + if (molID < 0) error->one(FLERR,"Invalid molecule ID in set command"); + } + + molecule[i] = molID; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_omega(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->omega_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set omega", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + iarg += 4; +} + +void Set::invoke_omega(Action *action) +{ + int nlocal = atom->nlocal; + double **omega = atom->omega; + + int varflag = action->varflag; + double xvalue,yvalue,zvalue; + if (!action->varflag1) xvalue = action->dvalue1; + if (!action->varflag2) yvalue = action->dvalue2; + if (!action->varflag3) zvalue = action->dvalue3; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) xvalue = vec1[i]; + if (action->varflag2) yvalue = vec2[i]; + if (action->varflag3) zvalue = vec3[i]; + } + + omega[i][0] = xvalue; + omega[i][1] = yvalue; + omega[i][2] = zvalue; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_quat(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "set quat", error); + int dimension = domain->dimension; + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (dimension == 2 && action->dvalue1 != 0.0) + error->one(FLERR,"Cannot set quaternion with xy components for 2d system"); + } + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else { + action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (dimension == 2 && action->dvalue2 != 0.0) + error->one(FLERR,"Cannot set quaternion with xy components for 2d system"); + } + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4,action); + else action->dvalue4 = utils::numeric(FLERR,arg[iarg+4],false,lmp); + + iarg += 5; +} + +void Set::invoke_quat(Action *action) +{ + int nlocal = atom->nlocal; + int *ellipsoid = atom->ellipsoid; + int *tri = atom->tri; + int *body = atom->body; + double **quat = atom->quat; + int quat_flag = atom->quat_flag; + + auto avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); + auto avec_tri = dynamic_cast(atom->style_match("tri")); + auto avec_body = dynamic_cast(atom->style_match("body")); + + int dimension = domain->dimension; + double radians,sintheta; + double *quat_one; + + int varflag = action->varflag; + double xvalue,yvalue,zvalue,theta; + if (!action->varflag1) xvalue = action->dvalue1; + if (!action->varflag2) yvalue = action->dvalue2; + if (!action->varflag3) zvalue = action->dvalue3; + if (!action->varflag4) theta = action->dvalue4; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (avec_ellipsoid && ellipsoid[i] >= 0) + quat_one = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (avec_tri && tri[i] >= 0) + quat_one = avec_tri->bonus[tri[i]].quat; + else if (avec_body && body[i] >= 0) + quat_one = avec_body->bonus[body[i]].quat; + else if (quat_flag) + quat_one = quat[i]; + else + error->one(FLERR,"Cannot set quaternion for atom that has none"); + + if (varflag) { + if (action->varflag1) xvalue = vec1[i]; + if (action->varflag2) yvalue = vec2[i]; + if (action->varflag3) zvalue = vec3[i]; + if (action->varflag4) theta = vec4[i]; + if (dimension == 2 && (xvalue != 0.0 || yvalue != 0.0)) + error->one(FLERR,"Cannot set quaternion with xy components for 2d system"); + } + + radians = MY_PI2 * theta/180.0; + sintheta = sin(radians); + quat_one[0] = cos(radians); + quat_one[1] = xvalue * sintheta; + quat_one[2] = yvalue * sintheta; + quat_one[3] = zvalue * sintheta; + MathExtra::qnormalize(quat_one); + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_quat_random(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set quat/random", error); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 <= 0) error->all(FLERR,"Invalid random number seed in set command"); + + iarg += 2; +} + +void Set::invoke_quat_random(Action *action) +{ + setrandom(QUAT_RANDOM,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_radius_election(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->eradius_flag) + error->all(FLERR, "Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "set radius/electron", error); + + if (utils::strmatch(arg[iarg + 1], "^v_")) + varparse(arg[iarg + 1], 1, action); + else { + action->dvalue1 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (action->dvalue1 < 0.0) + error->one(FLERR, iarg + 1, "Invalid electron radius {} in set radius/electron command", + action->dvalue1); + } + + iarg += 2; +} + +void Set::invoke_radius_election(Action *action) +{ + int nlocal = atom->nlocal; + double *eradius = atom->eradius; + + int varflag = action->varflag; + double radius; + if (!action->varflag1) radius = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + radius = vec1[i]; + if (radius < 0.0) error->one(FLERR,"Invalid electron radius in set command"); + } + + eradius[i] = radius; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_shape(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->ellipsoid_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set shape", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid shape in set command"); + } + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else { + action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (action->dvalue2 < 0.0) error->one(FLERR,"Invalid shape in set command"); + } + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else { + action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + if (action->dvalue3 < 0.0) error->one(FLERR,"Invalid shape in set command"); + } + + iarg += 4; +} + +void Set::invoke_shape(Action *action) +{ + int nlocal = atom->nlocal; + auto avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); + + int varflag = action->varflag; + double xvalue,yvalue,zvalue; + if (!action->varflag1) xvalue = action->dvalue1; + if (!action->varflag2) yvalue = action->dvalue2; + if (!action->varflag3) zvalue = action->dvalue3; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) xvalue = vec1[i]; + if (action->varflag1) yvalue = vec2[i]; + if (action->varflag1) zvalue = vec3[i]; + if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0) + error->one(FLERR,"Invalid shape in set command"); + } + + if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) + if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0) + error->one(FLERR,"Invalid shape in set command"); + + avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue); + } + + // update global ellipsoid count + + bigint nlocal_bonus = avec_ellipsoid->nlocal_bonus; + MPI_Allreduce(&nlocal_bonus,&atom->nellipsoids,1,MPI_LMP_BIGINT,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_smd_contact_radius(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->smd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set smd/contact/radius", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid smd/contact/radius in set command"); + } + + iarg += 2; +} + +void Set::invoke_smd_contact_radius(Action *action) +{ + int nlocal = atom->nlocal; + double *contact_radius = atom->contact_radius; + + int varflag = action->varflag; + double radius; + if (!action->varflag1) radius = action->dvalue1; + + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + radius = vec1[i]; + if (radius < 0.0) error->one(FLERR,"Invalid smd/contact/radius in set command"); + } + + contact_radius[i] = radius; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_smd_mass_density(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->smd_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set smd/mass/density", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) error->one(FLERR,"Invalid smd/mass/density in set command"); + } + + iarg += 2; +} + +void Set::invoke_smd_mass_density(Action *action) +{ + int nlocal = atom->nlocal; + double *rmass = atom->rmass; + double *vfrac = atom->vfrac; + + int varflag = action->varflag; + double density; + if (!action->varflag1) density = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + density = vec1[i]; + if (density < 0.0) error->one(FLERR,"Invalid smd/mass/density in set command"); + } + + rmass[i] = vfrac[i] * density; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_sph_cv(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->cv_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/cv", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_sph_cv(Action *action) +{ + int nlocal = atom->nlocal; + double *cv = atom->cv; + + int varflag = action->varflag; + double sph_cv; + if (!action->varflag1) sph_cv = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) sph_cv = vec1[i]; + cv[i] = sph_cv; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_sph_e(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->esph_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/e", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_sph_e(Action *action) +{ + int nlocal = atom->nlocal; + double *esph = atom->esph; + + int varflag = action->varflag; + double sph_e; + if (!action->varflag1) sph_e = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) sph_e = vec1[i]; + esph[i] = sph_e; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_sph_rho(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->rho_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/rho", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_sph_rho(Action *action) +{ + int nlocal = atom->nlocal; + double *rho = atom->rho; + + int varflag = action->varflag; + double sph_rho; + if (!action->varflag1) sph_rho = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) sph_rho = vec1[i]; + rho[i] = sph_rho; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_spin_atom(int &iarg, int narg, char **arg, Action *action) +{ + if ((strcmp(arg[iarg],"spin") == 0) && (comm->me == 0)) + error->warning(FLERR, "Set attribute spin is deprecated -- use spin/atom instead"); + if (!atom->sp_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "set spin/atom", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) + error->all(FLERR,"Invalid spin magnitude {} in set {} command", action->dvalue1, arg[iarg]); + } + + if (utils::strmatch(arg[iarg+2],"^v_")) varparse(arg[iarg+2],2,action); + else action->dvalue2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (utils::strmatch(arg[iarg+3],"^v_")) varparse(arg[iarg+3],3,action); + else action->dvalue3 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4,action); + else action->dvalue4 = utils::numeric(FLERR,arg[iarg+4],false,lmp); + + iarg += 5; +} + +void Set::invoke_spin_atom(Action *action) +{ + int nlocal = atom->nlocal; + double **sp = atom->sp; + double norm; + + int varflag = action->varflag; + double magnitude,xvalue,yvalue,zvalue; + if (!action->varflag1) magnitude = action->dvalue1; + if (!action->varflag2) xvalue = action->dvalue2; + if (!action->varflag3) yvalue = action->dvalue3; + if (!action->varflag4) zvalue = action->dvalue4; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + if (action->varflag1) magnitude = vec1[i]; + if (magnitude < 0.0) + error->one(FLERR,"Invalid spin magnitude in set command"); + if (action->varflag2) xvalue = vec2[i]; + if (action->varflag3) yvalue = vec3[i]; + if (action->varflag4) zvalue = vec4[i]; + } + + if ((xvalue == 0.0) && (yvalue == 0.0) && (zvalue == 0.0)) + error->all(FLERR,"At least one spin vector component must be non-zero"); + + norm = 1.0/sqrt(xvalue*xvalue+yvalue*yvalue+zvalue*zvalue); + sp[i][0] = norm*xvalue; + sp[i][1] = norm*yvalue; + sp[i][2] = norm*zvalue; + sp[i][3] = magnitude; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_spin_atom_random(int &iarg, int narg, char **arg, Action *action) +{ + if ((strcmp(arg[iarg], "spin/random") == 0) && (comm->me == 0)) + error->warning(FLERR, + "Set attribute spin/random is deprecated -- use spin/atom/random instead"); + if (!atom->sp_flag) + error->all(FLERR, "Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "set spin/atom/random", error); + + action->ivalue1 = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + action->dvalue1 = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + + if (action->ivalue1 <= 0) + error->all(FLERR, iarg + 1, "Invalid random number seed {} in set {} command", action->ivalue1, + arg[iarg]); + if (action->dvalue1 <= 0.0) + error->all(FLERR, iarg + 2, "Invalid spin magnitude {} in set {} command", action->dvalue1, + arg[iarg]); + iarg += 3; +} + +void Set::invoke_spin_atom_random(Action *action) +{ + setrandom(SPIN_ATOM_RANDOM, action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_spin_electron(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->spin_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set spin/electron", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 < -1 || action->ivalue1 > 3) + error->one(FLERR,"Invalid electron spin {} in set command", action->ivalue1); + } + + iarg += 2; +} + +void Set::invoke_spin_electron(Action *action) +{ + int nlocal = atom->nlocal; + int *spin = atom->spin; + + int varflag = action->varflag; + int ispin; + if (!action->varflag1) ispin = action->ivalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + ispin = static_cast (vec1[i]); + if (ispin < -1 || ispin > 3) + error->one(FLERR,"Invalid electron spin in set command"); + } + + atom->spin[i] = ispin; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_temperature(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->temperature_flag) + error->all(FLERR,"Cannot set this attribute for this atom style"); + if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid temperature in set command"); + } + + iarg += 2; +} + +void Set::invoke_temperature(Action *action) +{ + int nlocal = atom->nlocal; + double *temperature = atom->temperature; + + int varflag = action->varflag; + double temp; + if (!action->varflag1) temp = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + temp = vec1[i]; + if (temp < 0.0) error->one(FLERR,"Invalid temperature in set command"); + } + + temperature[i] = temp; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_theta(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->line_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set theta", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = DEG2RAD * utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_theta(Action *action) +{ + int nlocal = atom->nlocal; + int *line = atom->line; + + auto avec_line = dynamic_cast(atom->style_match("line")); + + int varflag = action->varflag; + double theta; + if (!action->varflag1) theta = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (line[i] < 0) error->one(FLERR,"Cannot set theta for atom which is not a line"); + if (varflag) theta = vec1[i]; + avec_line->bonus[atom->line[i]].theta = theta; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_theta_random(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->line_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set theta/random", error); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 <= 0) error->all(FLERR,"Invalid random number seed in set command"); + + iarg += 2; +} + +void Set::invoke_theta_random(Action *action) +{ + setrandom(THETA_RANDOM,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_tri(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->tri_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set tri", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 < 0.0) error->one(FLERR,"Invalid tri size in set command"); + } + + iarg += 2; +} + +void Set::invoke_tri(Action *action) +{ + int nlocal = atom->nlocal; + int *tri = atom->tri; + + auto avec_tri = dynamic_cast(atom->style_match("tri")); + + int varflag = action->varflag; + double trisize; + if (!action->varflag1) trisize = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; +#if 0 + // AK: this seems wrong. Isn't the set command supposed *make* this a triangle? + if (tri[i] < 0) error->one(FLERR,"Cannot set tri for atom which is not a triangle"); +#endif + + if (varflag) { + trisize = vec1[i]; + if (trisize < 0.0) error->one(FLERR,"Invalid tri size in set command"); + } + + avec_tri->set_equilateral(i,trisize); + } + + // update bonus tri count + + bigint nlocal_bonus = avec_tri->nlocal_bonus; + MPI_Allreduce(&nlocal_bonus,&atom->ntris,1,MPI_LMP_BIGINT,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_type(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set type", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::ATOM,lmp); + action->ivalue1 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue1 <= 0 || action->ivalue1 > atom->ntypes) + error->one(FLERR,"Invalid atom type in set command"); + } + + iarg += 2; +} + +void Set::invoke_type(Action *action) +{ + int nlocal = atom->nlocal; + int *type = atom->type; + + int varflag = action->varflag; + int itype; + if (!action->varflag1) itype = action->ivalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (action->varflag) { + itype = static_cast (vec1[i]); + if (itype <= 0 || itype > atom->ntypes) + error->one(FLERR,"Invalid atom type in set command"); + } + + type[i] = itype; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_type_fraction(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/fraction", error); + + // random seed must be ivalue1 for use in setrandom() + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::ATOM,lmp); + action->ivalue2 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue2 <= 0 || action->ivalue2 > atom->ntypes) + error->one(FLERR,"Invalid atom type in set command"); + + action->dvalue1 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (action->dvalue1 < 0.0 || action->dvalue1 > 1.0) + error->all(FLERR,"Invalid fraction in set command"); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+3],false,lmp); + if (action->ivalue1 <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + + iarg += 4; +} + +void Set::invoke_type_fraction(Action *action) +{ + setrandom(TYPE_FRACTION,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_type_ratio(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/ratio", error); + + // random seed must be ivalue1 for use in setrandom() + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::ATOM,lmp); + action->ivalue2 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue2 <= 0 || action->ivalue2 > atom->ntypes) + error->all(FLERR,"Invalid atom type in set command"); + + action->dvalue1 = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (action->dvalue1 < 0.0 || action->dvalue1 > 1.0) + error->all(FLERR,"Invalid fraction in set command"); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+3],false,lmp); + if (action->ivalue1 <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + + iarg += 4; +} + +void Set::invoke_type_ratio(Action *action) +{ + setrandom(TYPE_RATIO,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_type_subset(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "set type/subset", error); + + // random seed must be ivalue1 for use in setrandom() + + char *typestr = utils::expand_type(FLERR,arg[iarg+1],Atom::ATOM,lmp); + action->ivalue2 = utils::inumeric(FLERR,typestr?typestr:arg[iarg+1],false,lmp); + delete[] typestr; + if (action->ivalue2 <= 0 || action->ivalue2 > atom->ntypes) + error->all(FLERR,"Invalid atom type in set command"); + + action->bvalue1 = utils::bnumeric(FLERR,arg[iarg+2],false,lmp); + if (action->bvalue1 < 0) + error->all(FLERR,"Invalid subset size in set command"); + + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+3],false,lmp); + if (action->ivalue1 <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + + iarg += 4; +} + +void Set::invoke_type_subset(Action *action) +{ + setrandom(TYPE_SUBSET,action); +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_volume(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->vfrac_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set volume", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (action->dvalue1 <= 0.0) error->all(FLERR,"Invalid volume in set command"); + } + + iarg += 2; +} + +void Set::invoke_volume(Action *action) +{ + int nlocal = atom->nlocal; + double *vfrac = atom->vfrac; + + int varflag = action->varflag; + double vol; + if (!action->varflag1) vol = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + + if (varflag) { + vol = vec1[i]; + if (vol < 0.0) error->one(FLERR,"Invalid volume in set command"); + } + + vfrac[i] = vol; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_vx(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vx", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_vx(Action *action) +{ + int nlocal = atom->nlocal; + double **v = atom->v; + + int varflag = action->varflag; + double vx; + if (!action->varflag1) vx = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) vx = vec1[i]; + v[i][0] = vx; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_vy(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vy", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_vy(Action *action) +{ + int nlocal = atom->nlocal; + double **v = atom->v; + + int varflag = action->varflag; + double vy; + if (!action->varflag1) vy = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) vy = vec1[i]; + v[i][1] = vy; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_vz(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set vz", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_vz(Action *action) +{ + int nlocal = atom->nlocal; + double **v = atom->v; + + int varflag = action->varflag; + double vz; + if (!action->varflag1) vz = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) vz = vec1[i]; + v[i][2] = vz; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_x(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set x", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_x(Action *action) +{ + int nlocal = atom->nlocal; + double **x = atom->x; + + int varflag = action->varflag; + double coord; + if (!action->varflag1) coord = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) coord = vec1[i]; + x[i][0] = coord; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_y(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set y", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_y(Action *action) +{ + int nlocal = atom->nlocal; + double **x = atom->x; + + int varflag = action->varflag; + double coord; + if (!action->varflag1) coord = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) coord = vec1[i]; + x[i][1] = coord; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_z(int &iarg, int narg, char **arg, Action *action) +{ + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set z", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + iarg += 2; +} + +void Set::invoke_z(Action *action) +{ + int nlocal = atom->nlocal; + double **x = atom->x; + + int varflag = action->varflag; + double coord; + if (!action->varflag1) coord = action->dvalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) coord = vec1[i]; + x[i][2] = coord; + } +} + +/* ---------------------------------------------------------------------- */ + +void Set::process_custom(int &iarg, int narg, char **arg, Action *action) +{ + int flag,cols; + ArgInfo argi(arg[iarg],ArgInfo::DNAME|ArgInfo::INAME); + const char *pname = argi.get_name(); + + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set", error); + int index_custom = atom->find_custom(argi.get_name(),flag,cols); + if (index_custom < 0) + error->all(FLERR,"Set keyword or custom property {} does not exist",pname); + action->ivalue2 = index_custom; + + switch (argi.get_type()) { + + case ArgInfo::INAME: + if (flag != 0) error->all(FLERR,"Set command custom property {} is not integer",pname); + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + + if (argi.get_dim() == 0) { + if (cols > 0) + error->all(FLERR,"Set command custom integer property {} is not a vector",pname); + action->keyword = IVEC; + } else if (argi.get_dim() == 1) { + if (cols == 0) + error->all(FLERR,"Set command custom integer property {} is not an array",pname); + int icol_custom = argi.get_index1(); + if (icol_custom <= 0 || icol_custom > cols) + error->all(FLERR,"Set command per-atom custom integer array {} is accessed " + "out-of-range",pname); + action->ivalue3 = icol_custom; + action->keyword = IARRAY; + } else error->all(FLERR,"Illegal set command"); + break; + + case ArgInfo::DNAME: + if (flag != 1) error->all(FLERR,"Custom property {} is not floating-point",argi.get_name()); + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else action->dvalue1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + + if (argi.get_dim() == 0) { + if (cols > 0) + error->all(FLERR,"Set command custom double property {} is not a vector",pname); + action->keyword = DVEC; + } else if (argi.get_dim() == 1) { + if (cols == 0) + error->all(FLERR,"Set command custom double property {} is not an array",pname); + int icol_custom = argi.get_index1(); + if (icol_custom <= 0 || icol_custom > cols) + error->all(FLERR,"Set command per-atom custom double array {} is " + "accessed out-of-range",pname); + action->ivalue3 = icol_custom; + action->keyword = DARRAY; + } else error->all(FLERR,"Illegal set command"); + break; + + default: + error->all(FLERR,"Illegal set command"); + break; + } + + iarg += 2; +} + +void Set::invoke_custom(Action *action) +{ + int nlocal = atom->nlocal; + int ivalue; + double dvalue; + + int varflag = action->varflag; + int index_custom = action->ivalue2; + + if (action->keyword == IVEC) { + if (!varflag) ivalue = action->ivalue1; + int *ivector = atom->ivector[index_custom]; + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) ivalue = static_cast (vec1[i]); + ivector[i] = ivalue; + } + + } else if (action->keyword == DVEC) { + if (!varflag) dvalue = action->dvalue1; + double *dvector = atom->dvector[index_custom]; + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) dvalue = vec1[i]; + dvector[i] = dvalue; + } + + } else if (action->keyword == IARRAY) { + if (!varflag) ivalue = action->ivalue1; + int **iarray = atom->iarray[index_custom]; + int icol_custom = action->ivalue3 - 1; + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) ivalue = static_cast (vec1[i]); + iarray[i][icol_custom] = ivalue; + } + + } else if (action->keyword == DARRAY) { + if (!varflag) dvalue = action->dvalue1; + double **darray = atom->darray[index_custom]; + int icol_custom = action->ivalue3 - 1; + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) dvalue = vec1[i]; + darray[i][icol_custom] = dvalue; + } } } diff --git a/src/set.h b/src/set.h index e2be5d5c1e..65489d9932 100644 --- a/src/set.h +++ b/src/set.h @@ -26,29 +26,168 @@ namespace LAMMPS_NS { class Set : public Command { public: - Set(class LAMMPS *lmp) : Command(lmp) {}; + Set(class LAMMPS *lmp); + ~Set() override; + void command(int, char **) override; - private: - char *id; - int *select; - int style, ivalue, newtype, count, index_custom, icol_custom; - int ximage, yimage, zimage, ximageflag, yimageflag, zimageflag; - int cc_index; - bigint nsubset; - double dvalue, xvalue, yvalue, zvalue, wvalue, fraction; - - int varflag, varflag1, varflag2, varflag3, varflag4; - int ivar1, ivar2, ivar3, ivar4; - double *vec1, *vec2, *vec3, *vec4; - - int discflag; - + void process_args(int, int, char **); void selection(int); - void set(int); - void setrandom(int); - void topology(int); - void varparse(const char *, int); + void invoke_actions(); + + private: + int caller; // SETCOMMAND or FIXSET + + // params for atom selection + + int style; + char *id; + int nlo, nhi; + bigint nlobig, nhibig; + int groupbit; + class Region *region; + + // one Action = one keyword/value pair + + struct Action { + int keyword; + int argindex; + int count_select, count_action; + int varflag; + int varflag1, varflag2, varflag3, varflag4; + int ivar1, ivar2, ivar3, ivar4; + int ivalue1, ivalue2, ivalue3, ivalue4, ivalue5, ivalue6; + tagint tvalue1; + bigint bvalue1; + double dvalue1, dvalue2, dvalue3, dvalue4; + }; + + int naction, maxaction; + Action *actions; + + typedef void (Set::*FnPtrPack)(Action *); + FnPtrPack *invoke_choice; // list of ptrs to invoke functions + + double *vec1, *vec2, *vec3, *vec4; // storage for evaluated peratom variables + int varflag; // 1 if any action uses a per-atom variable + int varflag1, varflag2, varflag3, varflag4; // 1 if any action uses this variable + int maxvariable; // size of peratom variable vectors + + int *select; // flag for selected atoms + int maxselect; // size of select vector + + int count_select; // count of selected atoms on this proc + int count_action; // count of actions on this proc, only if different than selected + + // private functions + + void varparse(const char *, int, Action *); + void setrandom(int, Action *); + void topology(int, Action *); + + // customize by adding a process method + + void process_angle(int &, int, char **, Action *); + void process_angmom(int &, int, char **, Action *); + void process_bond(int &, int, char **, Action *); + void process_cc(int &, int, char **, Action *); + void process_charge(int &, int, char **, Action *); + void process_density(int &, int, char **, Action *); + void process_diameter(int &, int, char **, Action *); + void process_dihedral(int &, int, char **, Action *); + void process_dipole(int &, int, char **, Action *); + void process_dipole_random(int &, int, char **, Action *); + void process_dpd_theta(int &, int, char **, Action *); + void process_edpd_cv(int &, int, char **, Action *); + void process_edpd_temp(int &, int, char **, Action *); + void process_epsilon(int &, int, char **, Action *); + void process_image(int &, int, char **, Action *); + void process_improper(int &, int, char **, Action *); + void process_length(int &, int, char **, Action *); + void process_mass(int &, int, char **, Action *); + void process_mol(int &, int, char **, Action *); + void process_omega(int &, int, char **, Action *); + void process_quat(int &, int, char **, Action *); + void process_quat_random(int &, int, char **, Action *); + void process_radius_election(int &, int, char **, Action *); + void process_shape(int &, int, char **, Action *); + void process_smd_contact_radius(int &, int, char **, Action *); + void process_smd_mass_density(int &, int, char **, Action *); + void process_sph_cv(int &, int, char **, Action *); + void process_sph_e(int &, int, char **, Action *); + void process_sph_rho(int &, int, char **, Action *); + void process_spin_atom(int &, int, char **, Action *); + void process_spin_atom_random(int &, int, char **, Action *); + void process_spin_electron(int &, int, char **, Action *); + void process_temperature(int &, int, char **, Action *); + void process_theta(int &, int, char **, Action *); + void process_theta_random(int &, int, char **, Action *); + void process_tri(int &, int, char **, Action *); + void process_type(int &, int, char **, Action *); + void process_type_fraction(int &, int, char **, Action *); + void process_type_ratio(int &, int, char **, Action *); + void process_type_subset(int &, int, char **, Action *); + void process_volume(int &, int, char **, Action *); + void process_vx(int &, int, char **, Action *); + void process_vy(int &, int, char **, Action *); + void process_vz(int &, int, char **, Action *); + void process_x(int &, int, char **, Action *); + void process_y(int &, int, char **, Action *); + void process_z(int &, int, char **, Action *); + + void process_custom(int &, int, char **, Action *); + + // customize by adding an invoke method + + void invoke_angle(Action *); + void invoke_angmom(Action *); + void invoke_bond(Action *); + void invoke_cc(Action *); + void invoke_charge(Action *); + void invoke_density(Action *); + void invoke_diameter(Action *); + void invoke_dihedral(Action *); + void invoke_dipole(Action *); + void invoke_dipole_random(Action *); + void invoke_dpd_theta(Action *); + void invoke_edpd_cv(Action *); + void invoke_edpd_temp(Action *); + void invoke_epsilon(Action *); + void invoke_image(Action *); + void invoke_improper(Action *); + void invoke_length(Action *); + void invoke_mass(Action *); + void invoke_mol(Action *); + void invoke_omega(Action *); + void invoke_quat(Action *); + void invoke_quat_random(Action *); + void invoke_radius_election(Action *); + void invoke_shape(Action *); + void invoke_smd_contact_radius(Action *); + void invoke_smd_mass_density(Action *); + void invoke_sph_cv(Action *); + void invoke_sph_e(Action *); + void invoke_sph_rho(Action *); + void invoke_spin_atom(Action *); + void invoke_spin_atom_random(Action *); + void invoke_spin_electron(Action *); + void invoke_temperature(Action *); + void invoke_theta(Action *); + void invoke_theta_random(Action *); + void invoke_tri(Action *); + void invoke_type(Action *); + void invoke_type_fraction(Action *); + void invoke_type_ratio(Action *); + void invoke_type_subset(Action *); + void invoke_volume(Action *); + void invoke_vx(Action *); + void invoke_vy(Action *); + void invoke_vz(Action *); + void invoke_x(Action *); + void invoke_y(Action *); + void invoke_z(Action *); + + void invoke_custom(Action *); }; } // namespace LAMMPS_NS diff --git a/unittest/commands/test_set_property.cpp b/unittest/commands/test_set_property.cpp index c306ea4739..26635d225f 100644 --- a/unittest/commands/test_set_property.cpp +++ b/unittest/commands/test_set_property.cpp @@ -90,8 +90,9 @@ TEST_F(SetTest, NoBoxNoAtoms) ASSERT_EQ(compute->vector_atom[0], 0); TEST_FAILURE(".*ERROR: Illegal set command: need at least four.*", command("set type 1 x");); - TEST_FAILURE(".*ERROR: Unknown set command style: xxx.*", command("set xxx 1 x 0.0");); - TEST_FAILURE(".*ERROR: Set keyword or custom property yyy does not exist.*", + TEST_FAILURE(".*ERROR: Unknown set or fix set command style: xxx.*", + command("set xxx 1 x 0.0");); + TEST_FAILURE(".*ERROR: Unrecognized set or fix set command keyword yyy.*", command("set type 1 yyy 0.0");); TEST_FAILURE(".*ERROR: Cannot set attribute spin/atom for atom style atomic.*", @@ -447,9 +448,9 @@ TEST_F(SetTest, EffPackage) EXPECT_EQ(compute->array_atom[6][1], 0.5); EXPECT_EQ(compute->array_atom[7][1], 1.0); - TEST_FAILURE(".*ERROR on proc 0: Incorrect value for electron spin: 0.5.*", + TEST_FAILURE(".*Expected integer parameter instead of '0.5' in input script.*", command("set atom * spin/electron 0.5");); - TEST_FAILURE(".*ERROR on proc 0: Incorrect value for electron radius: -0.5.*", + TEST_FAILURE(".*ERROR on proc 0: Invalid electron radius -0.5 in set.*", command("set atom * radius/electron -0.5");); }