diff --git a/doc/src/compute_nbond_atom.rst b/doc/src/compute_nbond_atom.rst index f438836534..274d958a10 100644 --- a/doc/src/compute_nbond_atom.rst +++ b/doc/src/compute_nbond_atom.rst @@ -8,10 +8,17 @@ Syntax .. code-block:: LAMMPS - compute ID group-ID nbond/atom + compute ID group-ID nbond/atom keyword value * ID, group-ID are documented in :doc:`compute ` command * nbond/atom = style name of this compute command +* zero or more keyword/value pairs may be appended +* keyword = *bond/type* + + .. parsed-literal:: + + *bond/type* value = *btype* + *btype* = bond type included in count Examples """""""" @@ -19,6 +26,7 @@ Examples .. code-block:: LAMMPS compute 1 all nbond/atom + compute 1 all nbond/atom bond/type 2 Description """"""""""" @@ -31,6 +39,9 @@ the :doc:`Howto broken bonds ` page for more information. The number of bonds will be zero for atoms not in the specified compute group. This compute does not depend on Newton bond settings. +If the keyword *bond/type* is specified, only bonds of *btype* are +counted. + Output info """"""""""" diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index d31c305c20..2edd703299 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -8,12 +8,13 @@ Syntax .. parsed-literal:: - fix ID group-ID rheo/pressure style args + fix ID group-ID rheo/pressure type1 pstyle1 args1 ... typeN pstyleN argsN * ID, group-ID are documented in :doc:`fix ` command * rheo/pressure = style name of this fix command +* one or more types and pressure styles must be appended * types = lists of types (see below) -* style = *linear* or *taitwater* or *cubic* +* pstyle = *linear* or *taitwater* or *cubic* .. parsed-literal:: @@ -32,10 +33,41 @@ Examples Description """"""""""" -This fix... +This fix defines a pressure equation of state for RHEO particles. One can +define different equations of state for different atom types, but an +equation must be specified for every atom type. -Only one instance of fix rheo/pressure can be defined and the fix group must be set to all. +One first defines the atom *types*. A wild-card asterisk can be used in place +of or in conjunction with the *types* argument to set the coefficients for +multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with +no numeric values means all types from 1 to :math:`N`. A leading asterisk +means all types from 1 to n (inclusive). A trailing asterisk means all types +from m to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). +The *types* definition is followed by the pressure style, *pstyle*. Current +options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear +equation of state with a particle pressure :math:`P` calculated as + +.. math:: + + P = c (\rho - \rho_0) + +where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density, +and :math:`\rho` is the current density of a particle. The numerical values of +:math:`c` and :math:`\rho_0` are set in :doc:`fix rheo `. Style *cubic* +is a cubic equation of state which has an extra argument :math:`A_3`, + +.. math:: + + P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . + +Style *taitwater* is Tait's equation of state: + +.. math:: + + P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr]. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -52,7 +84,7 @@ Restrictions This fix must be used with an atom style that includes density such as atom_style rheo or rheo/thermal. This fix must be used in conjuction with :doc:`fix rheo `. The fix group must be -set to all. +set to all. Only one instance of fix rheo/pressure can be defined. This fix is part of the RHEO package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. @@ -60,8 +92,7 @@ LAMMPS was built with that package. See the :doc:`Build package Related commands """""""""""""""" -:doc:`fix rheo/viscosity `, -:doc:`fix rheo/thermal `, +:doc:`fix rheo `, :doc:`pair rheo `, :doc:`compute rheo/property/atom ` diff --git a/doc/src/fix_rheo_thermal.rst b/doc/src/fix_rheo_thermal.rst index b73aeb248e..2ece5521bc 100644 --- a/doc/src/fix_rheo_thermal.rst +++ b/doc/src/fix_rheo_thermal.rst @@ -48,26 +48,38 @@ Examples Description """"""""""" -This fix... +This fix performs time integration of temperature evolution for atom style +rheo/thermal. In addition, it defines multiple thermal properties of +particles and handles melting/solidification, if applicable. For more details +on phase transitions in RHEO, see :doc:`the RHEO howto `. -Each list consists of a series of type -ranges separated by commas. 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 -example, if M = the number of atom types, then an asterisk with no numeric -values means all types from 1 to M. A leading asterisk means all types -from 1 to n (inclusive). A trailing asterisk means all types from n to M -(inclusive). A middle asterisk means all types from m to n (inclusive). -Note that all atom types must be included in exactly one of the N collections. +For each atom type, one can define attributes for the *conductivity*, +*specific/heat*, *latent/heat*, and critical temperature (*Tfreeze*). +The conductivity and specific heat must be defined for all atom types. +The latent heat and critical temperature are optional. However, a +critical temperature must be defined to specify a latent heat. -While the *Tfreeze* keyword is optional, the *conductivity* and -*specific/heat* keywords are mandatory. +For each property, one must first define a list of atom types. A wild-card +asterisk can be used in place of or in conjunction with the *types* argument +to set the coefficients for multiple pairs of atom types. This takes the +form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom +types, then an asterisk with no numeric values means all types from 1 to +:math:`N`. A leading asterisk means all types from 1 to n (inclusive). +A trailing asterisk means all types from m to :math:`N` (inclusive). A +middle asterisk means all types from m to n (inclusive). -Multiple instances of this fix may be defined to apply different -properties to different groups. However, the union of fix groups -across all instances of fix rheo/thermal must cover all atoms. -If there are multiple instances of this fix, any intersections in -the fix groups will lead to incorrect thermal integration. +The *types* definition for each property is followed by the style. Currently, +the only option is *constant*. Style *constant* simply applies a constant value +of respective property to each particle of the assigned type. + +The *react* keyword controls whether bonds are created/deleted when particles +transition between a fluid and solid state. This option only applies to atom +types that have a defined value of *Tfreeze*. When a fluid particle's +temperature drops below *Tfreeze*, bonds of type *btype* are created between +nearby solid particles within a distance of *cut*. The particle's status also +swaps to a solid state. When a solid particle's temperature rises above +*Tfreeze*, all bonds of type *btype* are broken and the particle's tatus swaps +to a fluid state. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -84,7 +96,8 @@ Restrictions This fix must be used with an atom style that includes temperature, heatflow, and conductivity such as atom_tyle rheo/thermal This fix must be used in conjuction with :doc:`fix rheo ` with the -*thermal* setting. +*thermal* setting. The fix group must be set to all. Only one +instance of fix rheo/pressure can be defined. This fix is part of the RHEO package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. @@ -93,8 +106,6 @@ Related commands """""""""""""""" :doc:`fix rheo `, -:doc:`fix rheo/viscosity `, -:doc:`fix rheo/pressure `, :doc:`pair rheo `, :doc:`compute rheo/property/atom ` diff --git a/doc/src/fix_rheo_viscosity.rst b/doc/src/fix_rheo_viscosity.rst index 379b002de1..64c6539acc 100644 --- a/doc/src/fix_rheo_viscosity.rst +++ b/doc/src/fix_rheo_viscosity.rst @@ -8,22 +8,18 @@ Syntax .. parsed-literal:: - fix ID group-ID rheo/viscosity types style args ... + fix ID group-ID rheo/viscosity type1 pstyle1 args1 ... typeN pstyleN argsN * ID, group-ID are documented in :doc:`fix ` command * rheo/viscosity = style name of this fix command +* one or more types and viscosity styles must be appended * types = lists of types (see below) -* style = *constant* or *power* +* vstyle = *constant* .. parsed-literal:: - *constant* args = viscosity (mass/(length*time)) - *power* args = *eta* *gd0* *K* *npow* *tau0* - *eta* = (units) - *gd0* = (units) - *K* = (units) - *npow* = (units) - *tau0* = (units) + *constant* args = *eta* + *eta* = viscosity Examples """""""" @@ -31,20 +27,27 @@ Examples .. code-block:: LAMMPS fix 1 all rheo/viscosity * constant 1.0 - fix 1 all rheo/viscosity 1 constant 1.0 2 power 0.1 1e-2 0.5 0.01 + fix 1 all rheo/viscosity 1 constant 1.0 2 constant 2.0 Description """"""""""" -This fix... +This fix defines a viscosity for RHEO particles. One can define different +viscosities for different atom types, but a viscosity must be specified for +every atom type. -Multiple instances of this fix may be defined to apply different -properties to different groups. However, the union of fix groups -across all instances of fix rheo/viscosity must cover all atoms. -If there are multiple instances of this fix, any intersection -between fix groups will cause the viscosity for the affected atoms -to be calculated multiple times. Any such affected atoms will enabled -up with a viscosity calculated by the latest defined fix. +One first defines the atom *types*. A wild-card asterisk can be used in place +of or in conjunction with the *types* argument to set the coefficients for +multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with +no numeric values means all types from 1 to :math:`N`. A leading asterisk +means all types from 1 to n (inclusive). A trailing asterisk means all types +from m to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). + +The *types* definition is followed by the viscosity style, *vstyle*. Currently, +the only option is *constant*. Style *constant* simply applies a constant value +of the viscosity *eta* to each particle of the assigned type. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -60,7 +63,8 @@ Restrictions This fix must be used with an atom style that includes viscosity such as atom_style rheo or rheo/thermal. This fix must be used in -conjuction with :doc:`fix rheo `. +conjuction with :doc:`fix rheo `. The fix group must be +set to all. Only one instance of fix rheo/viscosity can be defined. This fix is part of the RHEO package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. @@ -69,8 +73,6 @@ Related commands """""""""""""""" :doc:`fix rheo `, -:doc:`fix rheo/pressure `, -:doc:`fix rheo/thermal `, :doc:`pair rheo `, :doc:`compute rheo/property/atom ` diff --git a/src/BPM/compute_nbond_atom.cpp b/src/BPM/compute_nbond_atom.cpp index 4f0fc4c3f0..af1cc2b9bc 100644 --- a/src/BPM/compute_nbond_atom.cpp +++ b/src/BPM/compute_nbond_atom.cpp @@ -25,7 +25,20 @@ using namespace LAMMPS_NS; ComputeNBondAtom::ComputeNBondAtom(LAMMPS *_lmp, int narg, char **arg) : Compute(_lmp, narg, arg), nbond(nullptr) { - if (narg < 3) utils::missing_cmd_args(FLERR, "compute nbond/atom", error); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute nbond/atom", error); + + btype = -1; + + iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg], "bond/type") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute nbond/atom bond/type", error); + btype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else { + error->all(FLERR, "Unknown compute nbond/type command {}", arg[iarg]); + } + } peratom_flag = 1; size_peratom_cols = 0; @@ -78,6 +91,7 @@ void ComputeNBondAtom::compute_peratom() for (i = 0; i < nlocal; i++) { for (j = 0; j < num_bond[i]; j++) { if (bond_type[i][j] <= 0) continue; + if (btype != -1 && bond_type[i][j] != btype) continue; k = atom->map(bond_atom[i][j]); if (k < 0) continue; diff --git a/src/BPM/compute_nbond_atom.h b/src/BPM/compute_nbond_atom.h index e0c2d7ce01..b55ef91e5d 100644 --- a/src/BPM/compute_nbond_atom.h +++ b/src/BPM/compute_nbond_atom.h @@ -35,7 +35,7 @@ class ComputeNBondAtom : public Compute { double memory_usage() override; private: - int nmax; + int nmax, btype; double *nbond; }; diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index 6a71136d9d..d366163459 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "comm.h" +#include "compute_rheo_surface.h" #include "domain.h" #include "error.h" #include "fix_bond_history.h" @@ -42,7 +43,7 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : - BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr) + BondBPM(_lmp), compute_surface(nullptr), k(nullptr), ecrit(nullptr), gamma(nullptr) { partial_flag = 1; comm_reverse = 1; @@ -279,7 +280,7 @@ void BondRHEOShell::compute(int eflag, int vflag) // Communicate changes in nbond if (newton_bond) comm->reverse_comm(this); - for(i = 0; i < nlocal; i++) { + for(int i = 0; i < nlocal; i++) { nbond[i] += dbond[i]; // If it has bonds, no shifting @@ -341,6 +342,18 @@ void BondRHEOShell::init_style() if (comm->ghost_velocity == 0) error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity"); + auto fixes = modify->get_fix_by_style("^rheo$"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell"); + class FixRHEO *fix_rheo = dynamic_cast(fixes[0]); + + if (!fix_rheo->surface_flag) error->all(FLERR, + "Bond rheo/shell requires surface calculation in fix rheo"); + compute_surface = fix_rheo->compute_surface; + + if (fix_rheo->oxidation_fix_defined != 1) + error->all(FLERR, "Need to define fix rheo/oxdiation to use bond rheo/shell"); + // check consistency in values (copy?), swap conditions to rsurf + if (!id_fix_bond_history) { id_fix_bond_history = utils::strdup("HISTORY_RHEO_SHELL"); fix_bond_history = dynamic_cast(modify->replace_fix( diff --git a/src/RHEO/bond_rheo_shell.h b/src/RHEO/bond_rheo_shell.h index 513d481eeb..a6a747a3f1 100644 --- a/src/RHEO/bond_rheo_shell.h +++ b/src/RHEO/bond_rheo_shell.h @@ -49,6 +49,8 @@ class BondRHEOShell : public BondBPM { int index_nb, nmax_store; char *id_fix; + class ComputeRHEOSurface *compute_surface; + void process_ineligibility(int, int); void allocate(); void store_data(); diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 2d6ff1e55a..bdc9951f90 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -29,6 +29,7 @@ #include "domain.h" #include "error.h" #include "fix_rheo.h" +#include "fix_rheo_oxidation.h" #include "fix_rheo_pressure.h" #include "fix_rheo_thermal.h" #include "memory.h" @@ -177,6 +178,11 @@ void ComputeRHEOPropertyAtom::init() fixes = modify->get_fix_by_style("rheo/pressure"); fix_pressure = dynamic_cast(fixes[0]); } + + if (shell_flag) { + fixes = modify->get_fix_by_style("rheo/oxidation"); + fix_oxidation = dynamic_cast(fixes[0]); + } } /* ---------------------------------------------------------------------- */ @@ -381,6 +387,21 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n) /* ---------------------------------------------------------------------- */ +void ComputeRHEOPropertyAtom::pack_nbond_shell(int n) +{ + int *nbond = fix_oxidation->nbond; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nbond[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputeRHEOPropertyAtom::pack_shift_v(int n) { double **vshift = compute_vshift->vshift; @@ -415,7 +436,7 @@ void ComputeRHEOPropertyAtom::pack_gradv(int n) void ComputeRHEOPropertyAtom::pack_atom_style(int n) { - atom->avec->pack_property_atom(avec_index[n],&buf[n],nvalues,groupbit); + atom->avec->pack_property_atom(avec_index[n], &buf[n], nvalues, groupbit); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index fd73b5883f..a66ef4ece6 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -63,6 +63,7 @@ class ComputeRHEOPropertyAtom : public Compute { class FixRHEO *fix_rheo; class FixRHEOPressure *fix_pressure; class FixRHEOThermal *fix_thermal; + class FixRHEOOxidation *fix_oxidation; class ComputeRHEOInterface *compute_interface; class ComputeRHEOKernel *compute_kernel; class ComputeRHEOSurface *compute_surface; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index b0cb1513fc..04b7d33541 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -48,12 +48,14 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : viscosity_fix_defined = 0; pressure_fix_defined = 0; thermal_fix_defined = 0; + oxidation_fix_defined = 0; thermal_flag = 0; rhosum_flag = 0; shift_flag = 0; interface_flag = 0; surface_flag = 0; + oxidation_flag = 0; int i; int n = atom->ntypes; @@ -252,6 +254,7 @@ void FixRHEO::setup(int /*vflag*/) thermal_fix_defined = 0; viscosity_fix_defined = 0; pressure_fix_defined = 0; + oxidation_fix_defined = 0; // Check fixes cover all atoms (may still fail if atoms are created) // FixRHEOPressure currently requires group all diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 33fd0084db..45f74ed4cd 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -51,6 +51,7 @@ class FixRHEO : public Fix { int shift_flag; int interface_flag; int surface_flag; + int oxidation_flag; int viscosity_fix_defined; int pressure_fix_defined; diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index 1628cef13f..f536df7842 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "atom_vec.h" +#include "compute_rheo_surface.h" #include "error.h" #include "fix_rheo.h" #include "force.h" @@ -36,9 +37,9 @@ enum {NONE, CONSTANT}; /* ---------------------------------------------------------------------- */ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) //, fix_bond_history(nullptr) + Fix(lmp, narg, arg), compute_surface(nullptr), fix_rheo(nullptr) { - if (narg != 5) error->all(FLERR,"Illegal fix command"); + if (narg != 6) error->all(FLERR,"Illegal fix command"); cut = utils::numeric(FLERR, arg[3], false, lmp); if (cut <= 0.0) error->all(FLERR, "Illegal bond cutoff {} in fix rheo/oxidation", cut); @@ -46,6 +47,9 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : btype = utils::inumeric(FLERR, arg[4], false, lmp); if (btype < 1 || btype > atom->nbondtypes) error->all(FLERR, "Illegal value {} for bond type in fix rheo/oxidation", btype); + rsurf = utils::numeric(FLERR, arg[5], false, lmp); + if (rsurf <= 0.0) error->all(FLERR, "Illegal surface distance {} in fix rheo/oxidation", cut); + cutsq = cut * cut; } @@ -72,22 +76,21 @@ void FixRHEOOxidation::init() auto fixes = modify->get_fix_by_style("^rheo$"); if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/oxidation"); fix_rheo = dynamic_cast(fixes[0]); - double cut_kernel = fix_rheo->h; - if (cut > cut_kernel) + if (cut > fix_rheo->h) error->all(FLERR, "Bonding length exceeds kernel cutoff"); + if (rsurf >= fix_rheo->h) + error->all(FLERR, "Surface distance must be less than kernel cutoff"); + if (!force->bond) error->all(FLERR, "Must define a bond style with fix rheo/oxidation"); if (!atom->avec->bonds_allow) error->all(FLERR, "Fix rheo/oxidation requires atom bonds"); - //// find instances of bond history to delete data - //histories = modify->get_fix_by_style("BOND_HISTORY"); - //for (auto &ihistory: histories) - // if (strcmp(histories[i]->id, "HISTORY_RHEO_SHELL") == 0) - // fix_bond_history = dynamic_cast(ihistory); -// - //if (!fix_bond_history) - // error->all(FLERR, "Must define bond style rheo/shell to use fix rheo/oxidation"); + int tmp1, tmp2; + index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); + if (index_nb == -1) + error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation"); + nbond = atom->ivector[index_nb]; // need a half neighbor list auto req = neighbor->add_request(this, NeighConst::REQ_DEFAULT); @@ -111,13 +114,14 @@ void FixRHEOOxidation::setup_pre_force(int /*vflag*/) if (!fix_rheo->surface_flag) error->all(FLERR, "fix rheo/oxidation requires surface calculation in fix rheo"); + compute_surface = fix_rheo->compute_surface; pre_force(0); } /* ---------------------------------------------------------------------- */ -void FixRHEOThermal::pre_force(int /*vflag*/) +void FixRHEOOxidation::pre_force(int /*vflag*/) { } @@ -135,9 +139,9 @@ void FixRHEOOxidation::post_integrate() tagint *tag = atom->tag; tagint **bond_atom = atom->bond_atom; - int *status = atom->status; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; + double *rsurface = compute_surface->rsurface; double **x = atom->x; inum = list->inum; @@ -148,7 +152,7 @@ void FixRHEOOxidation::post_integrate() // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (!(status[i] & STATUS_SURFACE)) continue; + if (rsurface[i] > rsurf) continue; tagi = tag[i]; xtmp = x[i][0]; @@ -162,7 +166,7 @@ void FixRHEOOxidation::post_integrate() j = jlist[jj]; j &= NEIGHMASK; - if (!(status[j] & STATUS_SURFACE)) continue; + if (rsurface[j] > rsurf) continue; tagj = tag[j]; delx = xtmp - x[j][0]; @@ -208,10 +212,4 @@ void FixRHEOOxidation::post_integrate() } } } - - //todo: - // allow option to create near-surface bonds - // extract # of bonds in property/atom - // check bond style shell used, get index to bonds, share with compute property atom - // add type option to compute nbond/atom } diff --git a/src/RHEO/fix_rheo_oxidation.h b/src/RHEO/fix_rheo_oxidation.h index ca36a6fdf9..5242d24460 100644 --- a/src/RHEO/fix_rheo_oxidation.h +++ b/src/RHEO/fix_rheo_oxidation.h @@ -36,14 +36,15 @@ class FixRHEOOxidation : public Fix { void setup_pre_force(int) override; void pre_force(int) override; void post_integrate() override; + int *nbond; private: - int btype; - double cut, cutsq; + int btype, index_nb; + double rsurf, cut, cutsq; class NeighList *list; + class ComputeRHEOSurface *compute_surface; class FixRHEO *fix_rheo; - //class FixBondHistory *fix_bond_history; }; } // namespace LAMMPS_NS