diff --git a/doc/src/compute_property_atom.rst b/doc/src/compute_property_atom.rst index 5dbf600c36..b03d6eb74e 100644 --- a/doc/src/compute_property_atom.rst +++ b/doc/src/compute_property_atom.rst @@ -128,9 +128,9 @@ Attributes *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom per-atom integer and floating-point vectors or arrays that have been added via the :doc:`fix property/atom ` command. When that command is used specific names are given to each attribute -which are the "name" portion of these attributes. For arrays *i2_name* -and *d2_name*, the column of the array must also be included following -the name in brackets (e.g., d2_xyz[2] or i2_mySpin[3]). +which are the "name" portion of these attributes. For arrays +*i2_name* and *d2_name*, the column of the array must also be included +following the name in brackets (e.g., d2_xyz[2] or i2_mySpin[3]). The additional quantities only accessible via this command, and not directly via the :doc:`dump custom ` command, are as follows. diff --git a/doc/src/dump.rst b/doc/src/dump.rst index 6d13b43200..2d1598e493 100644 --- a/doc/src/dump.rst +++ b/doc/src/dump.rst @@ -805,16 +805,16 @@ computes, fixes, or variables when they are evaluated, so this is a very general means of creating quantities to output to a dump file. The *i_name*, *d_name*, *i2_name*, *d2_name* attributes refer to -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[i], i2_mySpin[i], where :math:`i` is -in the range from 1 to :math:`M`, where :math:`M` is the number of -columns in the custom array). See the discussion above for how :math:`i` -can be specified with a wildcard asterisk to effectively specify -multiple values. +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[i], i2_mySpin[i], where +:math:`i` is in the range from 1 to :math:`M`, where :math:`M` is the +number of columns in the custom array). See the discussion above for +how :math:`i` can be specified with a wildcard asterisk to effectively +specify multiple values. See the :doc:`Modify ` page for information on how to add new compute and fix styles to LAMMPS to calculate per-atom quantities diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 92a78ee3c1..0a040a1655 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -71,6 +71,7 @@ Syntax feature functions = is_available(category,feature), is_active(category,feature), is_defined(category,id) atom value = id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i] atom vector = id, mass, type, mol, radius, q, x, y, z, vx, vy, vz, fx, fy, fz + custom atom property = i_name, d_name, i_name[i], d_name[i], i2_name[i], d2_name[i], i2_name[i][j], d_name[i][j] compute references = c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i] fix references = f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i] variable references = v_name, v_name[i] @@ -514,38 +515,40 @@ is a valid (though strange) variable formula: Specifically, a formula can contain numbers, constants, thermo keywords, math operators, math functions, group functions, region functions, special functions, feature functions, atom values, atom -vectors, compute references, fix references, and references to other +vectors, custom atom properties, compute references, fix references, and references to other variables. -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Number | 0.2, 100, 1.0e20, -15.4, etc | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Constant | PI, version, on, off, true, false, yes, no | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Thermo keywords | vol, pe, ebond, etc | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Math operators | (), -x, x+y, x-y, x\*y, x/y, x\^y, x%y, x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x \|\| y, x \|\^ y, !x | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Math functions | sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), logfreq2(x,y,z), logfreq3(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z) | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Group functions | count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), inertia(ID,dimdim), omega(ID,dim) | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Region functions | count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR) | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Special functions | sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label) | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Feature functions | is_available(category,feature), is_active(category,feature), is_defined(category,id) | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Atom values | id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i] | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Atom vectors | id, mass, type, mol, x, y, z, vx, vy, vz, fx, fy, fz, q | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Compute references | c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i] | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Fix references | f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i] | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| Other variables | v_name, v_name[i] | -+--------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Number | 0.2, 100, 1.0e20, -15.4, etc | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Constant | PI, version, on, off, true, false, yes, no | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Thermo keywords | vol, pe, ebond, etc | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Math operators | (), -x, x+y, x-y, x\*y, x/y, x\^y, x%y, x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x \|\| y, x \|\^ y, !x | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Math functions | sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), logfreq2(x,y,z), logfreq3(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z) | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Group functions | count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), inertia(ID,dimdim), omega(ID,dim) | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Region functions | count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR) | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Special functions | sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label) | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Feature functions | is_available(category,feature), is_active(category,feature), is_defined(category,id) | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Atom values | id[i], mass[i], type[i], mol[i], radius[i], q[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i] | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Atom vectors | id, mass, type, mol, radius, q, x, y, z, vx, vy, vz, fx, fy, fz | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Custom atom properties | i_name, d_name, i_name[i], d_name[i], i2_name[i], d2_name[i], i2_name[i][j], d_name[i][j] | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Compute references | c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i] | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Fix references | f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i] | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| Other variables | v_name, v_name[i] | ++------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ Most of the formula elements produce a scalar value. Some produce a global or per-atom vector of values. Global vectors can be produced @@ -1139,8 +1142,43 @@ defines molecule IDs. Note that many other atom attributes can be used as inputs to a variable by using the :doc:`compute property/atom -` command and then specifying a quantity from -that compute. +` command and then referencing that compute. + +---------- + +Custom atom properties +---------------------- + +Custom atom properties refer to 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 references. References beginning with *i* and *d* refer to +integer and floating point properties respectively. Per-atom vectors +are referenced by *i_name* and *d_name*; per-atom arrays are +referenced by *i2_name* and *d2_name*. + +The various allowed references to integer custom atom properties in +the variable formulas for equal-, vector-, and atom-style variables +are listed in the following table. References to floating point +custom atom properties are the same; just replace the leading "i" with +"d". + ++--------+---------------+------------------------------------------+ +| equal | i_name[I] | element of per-atom vector (I = atom ID) | +| equal | i2_name[I][J] | element of per-atom array (I = atom ID) | ++--------+---------------+------------------------------------------+ +| vector | i_name[I] | element of per-atom vector (I = atom ID) | +| vector | i2_name[I][J] | element of per-atom array (I = atom ID) | ++--------+---------------+------------------------------------------+ +| atom | i_name | per-atom vector | +| atom | i2_name[I] | column of per-atom array | ++--------+--------------+------------------------------------------+ + +The I and J indices in these custom atom property references can be +integers or can be a variable name, specified as v_name, where name is +the name of the variable. The rules for this syntax are the same as +for indices in the "Atom Values and Vectors" discussion above. ---------- diff --git a/src/atom.cpp b/src/atom.cpp index c22556ff9d..3ba5295fa7 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -2681,7 +2681,6 @@ int Atom::add_custom(const char *name, int flag, int cols) ianame[index] = utils::strdup(name); iarray = (int ***) memory->srealloc(iarray,niarray*sizeof(int **),"atom:iarray"); memory->create(iarray[index],nmax,cols,"atom:iarray"); - icols = (int *) memory->srealloc(icols,niarray*sizeof(int),"atom:icols"); icols[index] = cols; @@ -2692,10 +2691,10 @@ int Atom::add_custom(const char *name, int flag, int cols) daname[index] = utils::strdup(name); darray = (double ***) memory->srealloc(darray,ndarray*sizeof(double **),"atom:darray"); memory->create(darray[index],nmax,cols,"atom:darray"); - dcols = (int *) memory->srealloc(dcols,ndarray*sizeof(int),"atom:dcols"); dcols[index] = cols; } + if (index < 0) error->all(FLERR,"Invalid call to Atom::add_custom()"); return index; diff --git a/src/variable.cpp b/src/variable.cpp index 3bb49218fb..c769adf06c 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1345,6 +1345,7 @@ void Variable::copy(int narg, char **from, char **to) special function = sum(x),min(x), ... atom value = x[i], y[i], vx[i], ... atom vector = x, y, vx, ... + custom atom property = i/d_name, i/d_name[i], i/d2_name[i], i/d2_name[i][j] compute = c_ID, c_ID[i], c_ID[i][j] fix = f_ID, f_ID[i], f_ID[i][j] variable = v_name, v_name[i] @@ -1441,7 +1442,9 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // ---------------- // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][], - // v_name, v_name[], exp(), xcm(,), x, x[], PI, vol + // v_name, v_name[], exp(), xcm(,), x, x[], PI, vol, + // i/d_name, i/d_name[], i/d_name[][], + // i/d2_name, i/d2_name[], i/d2_name[][] // ---------------- } else if (isalpha(onechar)) { @@ -2126,8 +2129,118 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) } // ---------------- - // math/group/special/labelmap function or atom value/vector or - // constant or thermo keyword + // custom atom property + // ---------------- + + } else if (utils::strmatch(word,"^[id]2?_")) { + if (domain->box_exist == 0) + print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar); + + int index_custom,type_custom,cols_custom; + if (word[1] == '2') index_custom = atom->find_custom(word+3,type_custom,cols_custom); + else index_custom = atom->find_custom(word+2,type_custom,cols_custom); + + if (index_custom < 0) + print_var_error(FLERR,fmt::format("Invalid custom atom property reference {} in variable formula",word), + ivar); + + // parse zero or one or two trailing brackets + // point i beyond last bracket + // nbracket = # of bracket pairs + // index1,index2 = int inside each bracket pair, possibly an atom ID + + int nbracket; + tagint index1,index2; + if (str[i] != '[') nbracket = 0; + else { + nbracket = 1; + ptr = &str[i]; + index1 = int_between_brackets(ptr,1); + i = ptr-str+1; + if (str[i] == '[') { + nbracket = 2; + ptr = &str[i]; + index2 = int_between_brackets(ptr,1); + i = ptr-str+1; + } + } + + // custom atom property with no bracket + // can only mean use a per-atom vector + + if (nbracket == 0) { + if (cols_custom == 0) { + auto newtree = new Tree(); + treestack[ntreestack++] = newtree; + if (type_custom == 0) { + newtree->type = INTARRAY; + newtree->iarray = atom->ivector[index_custom]; + } else { + newtree->type = ATOMARRAY; + newtree->array = atom->dvector[index_custom]; + } + newtree->nstride = 1; + + } else if (cols_custom) { + print_var_error(FLERR,fmt::format("Invalid custom atom property reference {} in variable formula",word), + ivar); + } + + // custom atom property with one bracket + // can mean either extract a single value from a per-atom vector + // or use a column from a per-atom array + + } else if (nbracket == 1) { + + if (cols_custom == 0) { + if (type_custom == 0) { + custom2global(atom->ivector[index_custom],nullptr,1,index1, + tree,treestack,ntreestack,argstack,nargstack); + } else { + custom2global(nullptr,atom->dvector[index_custom],1,index1, + tree,treestack,ntreestack,argstack,nargstack); + } + + } else if (cols_custom) { + if (index1 <= 0 || index1 > cols_custom) + print_var_error(FLERR,fmt::format("Invalid custom atom property reference {} in variable formula",word), + ivar); + auto newtree = new Tree(); + treestack[ntreestack++] = newtree; + if (type_custom == 0) { + newtree->type = INTARRAY; + newtree->iarray = &atom->iarray[index_custom][0][index1-1]; + } else { + newtree->type = ATOMARRAY; + newtree->array = &atom->darray[index_custom][0][index1-1]; + } + newtree->nstride = cols_custom; + } + + // custom atom property with two brackets + // can only mean extract a single value from per-atom array + + } else if (nbracket == 2) { + if (cols_custom == 0) { + print_var_error(FLERR,fmt::format("Invalid custom atom property reference {} in variable formula",word), + ivar); + + } else if (cols_custom) { + if (index2 <= 0 || index2 > cols_custom) + print_var_error(FLERR,fmt::format("Invalid custom atom property reference {} in variable formula",word), + ivar); + if (type_custom == 0) { + custom2global(&atom->iarray[index_custom][0][index2],nullptr,cols_custom,index1, + tree,treestack,ntreestack,argstack,nargstack); + } else { + custom2global(nullptr,&atom->darray[index_custom][0][index2],cols_custom,index1, + tree,treestack,ntreestack,argstack,nargstack); + } + } + } + + // ---------------- + // math/group/special/labelmap function or atom value/vector or constant or thermo keyword // ---------------- } else { @@ -4646,21 +4759,21 @@ int Variable::feature_function(char *word, char *contents, Tree **tree, Tree **t extract a global value from a per-atom quantity in a formula flag = 0 -> word is an atom vector flag = 1 -> vector is a per-atom compute or fix quantity with nstride - id = global ID of atom, converted to local index + id = global ID of atom, converted here to local index via atom map push result onto tree or arg stack customize by adding an atom vector: - id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q + id,mass,type,mol,radius,q,x,y,z,vx,vy,vz,fx,fy,fz,q ------------------------------------------------------------------------- */ void Variable::peratom2global(int flag, char *word, double *vector, int nstride, tagint id, Tree **tree, Tree **treestack, int &ntreestack, double *argstack, int &nargstack) { - // error check for ID larger than any atom - // int_between_brackets() already checked for ID <= 0 - if (atom->map_style == Atom::MAP_NONE) error->all(FLERR, "Indexed per-atom vector in variable formula without atom map"); + // error check for ID larger than any atom + // int_between_brackets() already checked for ID <= 0 + if (id > atom->map_tag_max) error->all(FLERR,"Variable atom ID is too large"); @@ -4673,17 +4786,27 @@ void Variable::peratom2global(int flag, char *word, double *vector, int nstride, if (index >= 0 && index < atom->nlocal) { if (flag == 0) { - if (strcmp(word,"id") == 0) mine = atom->tag[index]; - else if (strcmp(word,"mass") == 0) { + if (strcmp(word,"id") == 0) { + mine = atom->tag[index]; + } else if (strcmp(word,"mass") == 0) { if (atom->rmass) mine = atom->rmass[index]; else mine = atom->mass[atom->type[index]]; - } - else if (strcmp(word,"type") == 0) mine = atom->type[index]; - else if (strcmp(word,"mol") == 0) { + } else if (strcmp(word,"type") == 0) { + mine = atom->type[index]; + } else if (strcmp(word,"mol") == 0) { if (!atom->molecule_flag) error->one(FLERR,"Variable uses atom property that isn't allocated"); mine = atom->molecule[index]; + } else if (strcmp(word,"radius") == 0) { + if (!atom->radius_flag) + error->one(FLERR,"Variable uses atom property that isn't allocated"); + mine = atom->radius[index]; + } else if (strcmp(word,"q") == 0) { + if (!atom->q_flag) + error->one(FLERR,"Variable uses atom property that isn't allocated"); + mine = atom->q[index]; } + else if (strcmp(word,"x") == 0) mine = atom->x[index][0]; else if (strcmp(word,"y") == 0) mine = atom->x[index][1]; else if (strcmp(word,"z") == 0) mine = atom->x[index][2]; @@ -4693,11 +4816,6 @@ void Variable::peratom2global(int flag, char *word, double *vector, int nstride, else if (strcmp(word,"fx") == 0) mine = atom->f[index][0]; else if (strcmp(word,"fy") == 0) mine = atom->f[index][1]; else if (strcmp(word,"fz") == 0) mine = atom->f[index][2]; - else if (strcmp(word,"q") == 0) { - if (!atom->q_flag) - error->one(FLERR,"Variable uses atom property that isn't allocated"); - mine = atom->q[index]; - } else error->one(FLERR,"Invalid atom vector {} in variable formula", word); } else mine = vector[index*nstride]; @@ -4715,11 +4833,54 @@ void Variable::peratom2global(int flag, char *word, double *vector, int nstride, } else argstack[nargstack++] = value; } +/* ---------------------------------------------------------------------- + extract a global value from a custom atom property in a formula + ivector = ptr to integer per-atom property with nstride + dvector = ptr to floating-point per-atom property with nstride + exactly one if ivector/dvector is non-NULL + id = global ID of atom, converted here to local index via atom map + push result onto tree or arg stack +------------------------------------------------------------------------- */ + +void Variable::custom2global(int *ivector, double *dvector, int nstride, tagint id, Tree **tree, + Tree **treestack, int &ntreestack, double *argstack, int &nargstack) +{ + if (atom->map_style == Atom::MAP_NONE) + error->all(FLERR, "Referenced custom atom property in variable formula without atom map"); + + // error check for ID larger than any atom + // int_between_brackets() already checked for ID <= 0 + + if (id > atom->map_tag_max) + error->all(FLERR,"Variable atom ID is too large"); + + // if ID does not exist, index will be -1 for all procs, + // and mine will be set to 0.0 + + int index = atom->map(id); + + double mine; + if (index >= 0 && index < atom->nlocal) { + if (ivector) mine = ivector[index*nstride]; + else mine = dvector[index*nstride]; + } else mine = 0.0; + + double value; + MPI_Allreduce(&mine,&value,1,MPI_DOUBLE,MPI_SUM,world); + + if (tree) { + auto newtree = new Tree(); + newtree->type = VALUE; + newtree->value = value; + treestack[ntreestack++] = newtree; + } else argstack[nargstack++] = value; +} + /* ---------------------------------------------------------------------- check if word matches an atom vector return 1 if yes, else 0 customize by adding an atom vector: - id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q + id,mass,type,mol,radius,q,x,y,z,vx,vy,vz,fx,fy,fz ------------------------------------------------------------------------- */ int Variable::is_atom_vector(char *word) @@ -4727,8 +4888,9 @@ int Variable::is_atom_vector(char *word) if (strcmp(word,"id") == 0) return 1; if (strcmp(word,"mass") == 0) return 1; if (strcmp(word,"type") == 0) return 1; - if (strcmp(word,"radius") == 0) return 1; if (strcmp(word,"mol") == 0) return 1; + if (strcmp(word,"radius") == 0) return 1; + if (strcmp(word,"q") == 0) return 1; if (strcmp(word,"x") == 0) return 1; if (strcmp(word,"y") == 0) return 1; if (strcmp(word,"z") == 0) return 1; @@ -4738,7 +4900,6 @@ int Variable::is_atom_vector(char *word) if (strcmp(word,"fx") == 0) return 1; if (strcmp(word,"fy") == 0) return 1; if (strcmp(word,"fz") == 0) return 1; - if (strcmp(word,"q") == 0) return 1; return 0; } @@ -4747,7 +4908,7 @@ int Variable::is_atom_vector(char *word) push result onto tree word = atom vector customize by adding an atom vector: - id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q + id,mass,type,mol,radius,q,x,y,z,vx,vy,vz,fx,fy,fz ------------------------------------------------------------------------- */ void Variable::atom_vector(char *word, Tree **tree, Tree **treestack, int &ntreestack) diff --git a/src/variable.h b/src/variable.h index c0aeaebd37..8acfa5bcc7 100644 --- a/src/variable.h +++ b/src/variable.h @@ -144,6 +144,7 @@ class Variable : protected Pointers { int special_function(char *, char *, Tree **, Tree **, int &, double *, int &, int); int feature_function(char *, char *, Tree **, Tree **, int &, double *, int &, int); void peratom2global(int, char *, double *, int, tagint, Tree **, Tree **, int &, double *, int &); + void custom2global(int *, double *, int, tagint, Tree **, Tree **, int &, double *, int &); int is_atom_vector(char *); void atom_vector(char *, Tree **, Tree **, int &); int parse_args(char *, char **);