diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index aaf706b5df..ee03b7e245 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -54,6 +54,7 @@ OPT. * :doc:`oxdna2/fene ` * :doc:`oxrna2/fene ` * :doc:`quartic (o) ` + * :doc:`rheo/shell ` * :doc:`special ` * :doc:`table (o) ` diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 0352ad5374..394a5bee3a 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -122,6 +122,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`reduce ` * :doc:`reduce/chunk ` * :doc:`reduce/region ` + * :doc:`rheo/property/atom ` * :doc:`rigid/local ` * :doc:`saed ` * :doc:`slcsa/atom ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index e89e302673..7053c4809a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -203,6 +203,11 @@ OPT. * :doc:`reaxff/species (k) ` * :doc:`recenter ` * :doc:`restrain ` + * :doc:`rheo ` + * :doc:`rheo/oxidation ` + * :doc:`rheo/pressure ` + * :doc:`rheo/thermal ` + * :doc:`rheo/viscosity ` * :doc:`rhok ` * :doc:`rigid (o) ` * :doc:`rigid/meso ` diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index e7761e7bee..9b56e92819 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -257,6 +257,8 @@ OPT. * :doc:`reaxff (ko) ` * :doc:`rebo (io) ` * :doc:`resquared (go) ` + * :doc:`rheo ` + * :doc:`rheo/solid ` * :doc:`saip/metal (t) ` * :doc:`sdpd/taitwater/isothermal ` * :doc:`smatb ` diff --git a/doc/src/bond_rheo_shell.rst b/doc/src/bond_rheo_shell.rst new file mode 100644 index 0000000000..7f6bab1eab --- /dev/null +++ b/doc/src/bond_rheo_shell.rst @@ -0,0 +1,226 @@ +.. index:: bond_style rheo/shell + +bond_style rheo/shell command +============================= + +Syntax +"""""" + +.. code-block:: LAMMPS + + bond_style rheo/shell keyword value attribute1 attribute2 ... + +* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break* + + .. parsed-literal:: + + *store/local* values = fix_ID N attributes ... + * fix_ID = ID of associated internal fix to store data + * N = prepare data for output every this many timesteps + * attributes = zero or more of the below attributes may be appended + + *id1, id2* = IDs of 2 atoms in the bond + *time* = the timestep the bond broke + *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units) + *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units) + + *overlay/pair* value = *yes* or *no* + bonded particles will still interact with pair forces + + *smooth* value = *yes* or *no* + smooths bond forces near the breaking point + + *normalize* value = *yes* or *no* + normalizes bond forces by the reference length + + *break* value = *yes* or *no* + indicates whether bonds break during a run + +Examples +"""""""" + +.. code-block:: LAMMPS + + bond_style bpm/spring + bond_coeff 1 1.0 0.05 0.1 + + bond_style bpm/spring myfix 1000 time id1 id2 + dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3] + dump_modify 1 write_header no + +Description +""""""""""" + +.. versionadded:: 4May2022 + +The *bpm/spring* bond style computes forces based on +deviations from the initial reference state of the two atoms. The +reference state is stored by each bond when it is first computed in +the setup of a run. Data is then preserved across run commands and is +written to :doc:`binary restart files ` such that restarting +the system will not reset the reference state of a bond. + +This bond style only applies central-body forces which conserve the +translational and rotational degrees of freedom of a bonded set of +particles based on a model described by Clemmer and Robbins +:ref:`(Clemmer) `. The force has a magnitude of + +.. math:: + + F = k (r - r_0) w + +where :math:`k` is a stiffness, :math:`r` is the current distance +and :math:`r_0` is the initial distance between the two particles, and +:math:`w` is an optional smoothing factor discussed below. Bonds will +break at a strain of :math:`\epsilon_c`. This is done by setting +the bond type to 0 such that forces are no longer computed. + +An additional damping force is applied to the bonded +particles. This forces is proportional to the difference in the +normal velocity of particles using a similar construction as +dissipative particle dynamics :ref:`(Groot) `: + +.. math:: + + F_D = - \gamma w (\hat{r} \bullet \vec{v}) + +where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the +radial normal vector, and :math:`\vec{v}` is the velocity difference +between the two particles. + +The smoothing factor :math:`w` can be added or removed by setting the +*smooth* keyword to *yes* or *no*, respectively. It is constructed such +that forces smoothly go to zero, avoiding discontinuities, as bonds +approach the critical strain + +.. math:: + + w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 . + +The following coefficients must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data +` or :doc:`read_restart ` commands: + +* :math:`k` (force/distance units) +* :math:`\epsilon_c` (unit less) +* :math:`\gamma` (force/velocity units) + +If the *normalize* keyword is set to *yes*, the elastic bond force will be +normalized by :math:`r_0` such that :math:`k` must be given in force units. + +By default, pair forces are not calculated between bonded particles. +Pair forces can alternatively be overlaid on top of bond forces by setting +the *overlay/pair* keyword to *yes*. These settings require specific +:doc:`special_bonds ` settings described in the +restrictions. Further details can be found in the :doc:`how to ` +page on BPMs. + +.. versionadded:: 28Mar2023 + +If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break +during a simulation run. This will prevent some unnecessary calculation. +However, if a bond reaches a strain greater than :math:`\epsilon_c`, +it will trigger an error. + +If the *store/local* keyword is used, an internal fix will track bonds that +break during the simulation. Whenever a bond breaks, data is processed +and transferred to an internal fix labeled *fix_ID*. This allows the +local data to be accessed by other LAMMPS commands. Following this optional +keyword, a list of one or more attributes is specified. These include the +IDs of the two atoms in the bond. The other attributes for the two atoms +include the timestep during which the bond broke and the current/initial +center of mass position of the two atoms. + +Data is continuously accumulated over intervals of *N* +timesteps. At the end of each interval, all of the saved accumulated +data is deleted to make room for new data. Individual datum may +therefore persist anywhere between *1* to *N* timesteps depending on +when they are saved. This data can be accessed using the *fix_ID* and a +:doc:`dump local ` command. To ensure all data is output, +the dump frequency should correspond to the same interval of *N* +timesteps. A dump frequency of an integer multiple of *N* can be used +to regularly output a sample of the accumulated data. + +Note that when unbroken bonds are dumped to a file via the +:doc:`dump local ` command, bonds with type 0 (broken bonds) +are not included. +The :doc:`delete_bonds ` command can also be used to +query the status of broken bonds or permanently delete them, e.g.: + +.. code-block:: LAMMPS + + delete_bonds all stats + delete_bonds all bond 0 remove + +---------- + +Restart and other info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This bond style writes the reference state of each bond to +:doc:`binary restart files `. Loading a restart +file will properly restore bonds. However, the reference state is NOT +written to data files. Therefore reading a data file will not +restore bonds and will cause their reference states to be redefined. + +If the *store/local* option is used, an internal fix will calculate +a local vector or local array depending on the number of input values. +The length of the vector or number of rows in the array is the number +of recorded, broken bonds. If a single input is specified, a local +vector is produced. If two or more inputs are specified, a local array +is produced where the number of columns = the number of inputs. The +vector or array can be accessed by any command that uses local values +from a compute as input. See the :doc:`Howto output ` page +for an overview of LAMMPS output options. + +The vector or array will be floating point values that correspond to +the specified attribute. + +The single() function of this bond style returns 0.0 for the energy +of a bonded interaction, since energy is not conserved in these +dissipative potentials. The single() function also calculates an +extra bond quantity, the initial distance :math:`r_0`. This +extra quantity can be accessed by the +:doc:`compute bond/local ` command as *b1*\ . + +Restrictions +"""""""""""" + +This bond style is part of the BPM package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +By default if pair interactions between bonded atoms are to be disabled, +this bond style requires setting + +.. code-block:: LAMMPS + + special_bonds lj 0 1 1 coul 1 1 1 + +and :doc:`newton ` must be set to bond off. If the *overlay/pair* +keyword is set to *yes*, this bond style alternatively requires setting + +.. code-block:: LAMMPS + + special_bonds lj/coul 1 1 1 + +Related commands +"""""""""""""""" + +:doc:`bond_coeff `, :doc:`pair bpm/spring ` + +Default +""""""" + +The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes* + +---------- + +.. _fragment-Clemmer: + +**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022). + +.. _Groot4: + +**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997). diff --git a/doc/src/fix_rheo_oxidation.rst b/doc/src/fix_rheo_oxidation.rst new file mode 100644 index 0000000000..a747ec582f --- /dev/null +++ b/doc/src/fix_rheo_oxidation.rst @@ -0,0 +1,81 @@ +.. index:: fix rheo/oxidation + +fix rheo/oxidation command +========================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID rheo/oxidation cut btype + +* ID, group-ID are documented in :doc:`fix ` command +* rheo/oxidation = style name of this fix command +* cut = maximum bond length (distance units) +* btype = type of bonds created + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all rheo/oxidation 1.5 2 + +Description +""""""""""" + +This fix... + +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. + +While the *Tfreeze* keyword is optional, the *conductivity* and +*specific/heat* keywords are mandatory. + +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. + +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 +"""""""""""" + +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. + +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. + +Related commands +"""""""""""""""" + +:doc:`fix rheo `, +:doc:`fix rheo/viscosity `, +:doc:`fix rheo/pressure `, +:doc:`pair rheo `, +:doc:`compute rheo/property/atom ` + +Default +""""""" + +none diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index b484df7fab..34554497ad 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -50,10 +50,10 @@ BondBPM::BondBPM(LAMMPS *_lmp) : // this is so final order of Modify:fix will conform to input script // BondHistory technically only needs this if updateflag = 1 - id_fix_dummy = utils::strdup("BPM_DUMMY"); + id_fix_dummy = utils::strdup(fmt::format("BPM_DUMMY_{}", instance_total)); modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy)); - id_fix_dummy2 = utils::strdup("BPM_DUMMY2"); + id_fix_dummy2 = utils::strdup(fmt::format("BPM_DUMMY2_{}", instance_total)); modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2)); } diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp new file mode 100644 index 0000000000..6a9f99d8b1 --- /dev/null +++ b/src/RHEO/bond_rheo_shell.cpp @@ -0,0 +1,507 @@ +/* ---------------------------------------------------------------------- + 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 "bond_rheo_shell.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_bond_history.h" +#include "fix_rheo.h" +#include "fix_store_local.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +#define EPSILON 1e-10 + +using namespace LAMMPS_NS; +using namespace RHEO_NS; + +/* ---------------------------------------------------------------------- */ + +BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : + BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr) +{ + partial_flag = 1; + + tform = rmax = -1; + + single_extra = 1; + svector = new double[1]; +} + +/* ---------------------------------------------------------------------- */ + +BondRHEOShell::~BondRHEOShell() +{ + delete[] svector; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(ecrit); + memory->destroy(gamma); + } +} + +/* ---------------------------------------------------------------------- + Store data for a single bond - if bond added after LAMMPS init (e.g. pour) +------------------------------------------------------------------------- */ + +double BondRHEOShell::store_bond(int n, int i, int j) +{ + double **bondstore = fix_bond_history->bondstore; + tagint *tag = atom->tag; + + bondstore[n][0] = 0.0; + bondstore[n][1] = 0.0; + + if (i < atom->nlocal) { + for (int m = 0; m < atom->num_bond[i]; m++) { + if (atom->bond_atom[i][m] == tag[j]) { + fix_bond_history->update_atom_value(i, m, 0, 0.0); + fix_bond_history->update_atom_value(i, m, 1, 0.0); + } + } + } + + if (j < atom->nlocal) { + for (int m = 0; m < atom->num_bond[j]; m++) { + if (atom->bond_atom[j][m] == tag[i]) { + fix_bond_history->update_atom_value(j, m, 0, 0.0); + fix_bond_history->update_atom_value(j, m, 1, 0.0); + } + } + } + + return 0.0; +} + +/* ---------------------------------------------------------------------- + Store data for all bonds called once +------------------------------------------------------------------------- */ + +void BondRHEOShell::store_data() +{ + int i, j, m, type; + int **bond_type = atom->bond_type; + + for (i = 0; i < atom->nlocal; i++) { + for (m = 0; m < atom->num_bond[i]; m++) { + type = bond_type[i][m]; + + //Skip if bond was turned off + if (type < 0) continue; + + // map to find index n + j = atom->map(atom->bond_atom[i][m]); + if (j == -1) error->one(FLERR, "Atom missing in BPM bond"); + + fix_bond_history->update_atom_value(i, m, 0, 0.0); + fix_bond_history->update_atom_value(i, m, 1, 0.0); + } + } + + fix_bond_history->post_neighbor(); +} + +/* ---------------------------------------------------------------------- */ + +void BondRHEOShell::compute(int eflag, int vflag) +{ + + if (!fix_bond_history->stored_flag) { + fix_bond_history->stored_flag = true; + store_data(); + } + + int i1, i2, itmp, n, type; + double delx, dely, delz, delvx, delvy, delvz; + double e, rsq, r, r0, rinv, dr, fbond, dot, t; + double dt = update->dt; + + ev_init(eflag, vflag); + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + tagint *tag = atom->tag; + int *status = atom->status; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + double **bondstore = fix_bond_history->bondstore; + + for (n = 0; n < nbondlist; n++) { + + // skip bond if already broken + if (bondlist[n][2] <= 0) continue; + + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + r0 = bondstore[n][0]; + t = bondstore[n][1]; + + // Ensure pair is always ordered to ensure numerical operations + // are identical to minimize the possibility that a bond straddling + // an mpi grid (newton off) doesn't break on one proc but not the other + if (tag[i2] < tag[i1]) { + itmp = i1; + i1 = i2; + i2 = itmp; + } + + // If bond hasn't been set - set timer to zero + if (t < EPSILON || std::isnan(t)) r0 = store_bond(n, i1, i2); + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + rsq = delx * delx + dely * dely + delz * delz; + r = sqrt(rsq); + + // Bond has not yet formed, check if in range + update timer + if (t < tform) { + + // Check if eligible + if (r > rmax || !(status[i1] & STATUS_SURFACE) || !(status[i2] & STATUS_SURFACE)) { + bondlist[n][2] = 0; + process_ineligibility(i1, i2); + continue; + } + + // Check ellapsed time + bondstore[n][1] += dt; + if (bondstore[n][1] >= tform) { + bondstore[n][0] = r; + r0 = r; + } else { + continue; + } + } + + e = (r - r0) / r0; + if (fabs(e) > ecrit[type]) { + bondlist[n][2] = 0; + process_broken(i1, i2); + continue; + } + + rinv = 1.0 / r; + dr = r - r0; + fbond = 2 * k[type] * (-dr + dr * dr * dr / (r0 * r0 * ecrit[type] * ecrit[type])); + + delvx = v[i1][0] - v[i2][0]; + delvy = v[i1][1] - v[i2][1]; + delvz = v[i1][2] - v[i2][2]; + dot = delx * delvx + dely * delvy + delz * delvz; + fbond -= gamma[type] * dot * rinv; + fbond *= rinv; + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx * fbond; + f[i1][1] += dely * fbond; + f[i1][2] += delz * fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx * fbond; + f[i2][1] -= dely * fbond; + f[i2][2] -= delz * fbond; + } + + if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondRHEOShell::allocate() +{ + allocated = 1; + const int np1 = atom->nbondtypes + 1; + + memory->create(k, np1, "bond:k"); + memory->create(ecrit, np1, "bond:ecrit"); + memory->create(gamma, np1, "bond:gamma"); + + memory->create(setflag, np1, "bond:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondRHEOShell::coeff(int narg, char **arg) +{ + if (narg != 4) error->all(FLERR, "Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + double ecrit_one = utils::numeric(FLERR, arg[2], false, lmp); + double gamma_one = utils::numeric(FLERR, arg[3], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + ecrit[i] = ecrit_one; + gamma[i] = gamma_one; + setflag[i] = 1; + count++; + + if (1.0 + ecrit[i] > max_stretch) max_stretch = 1.0 + ecrit[i]; + } + + if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + check for correct settings and create fix +------------------------------------------------------------------------- */ + +void BondRHEOShell::init_style() +{ + if (comm->ghost_velocity == 0) + error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity"); + + if (!id_fix_bond_history) { + id_fix_bond_history = utils::strdup("HISTORY_RHEO_SHELL"); + fix_bond_history = dynamic_cast(modify->replace_fix( + id_fix_dummy2, fmt::format("{} all BOND_HISTORY 1 2", id_fix_bond_history), 1)); + delete[] id_fix_dummy2; + id_fix_dummy2 = nullptr; + } + + // Reproduce standard functions of BondBPM, removing special restrictions + // Since this bond is intended to be created by fix rheo/oxidation, it + // ignores special statuses + + if (id_fix_store_local) { + auto ifix = modify->get_fix_by_id(id_fix_store_local); + if (!ifix) error->all(FLERR, "Cannot find fix STORE/LOCAL id {}", id_fix_store_local); + if (strcmp(ifix->style, "STORE/LOCAL") != 0) + error->all(FLERR, "Incorrect fix style matched, not STORE/LOCAL: {}", ifix->style); + fix_store_local = dynamic_cast(ifix); + fix_store_local->nvalues = nvalues; + } + + id_fix_update = nullptr; + + if (force->angle || force->dihedral || force->improper) + error->all(FLERR, "Bond style rheo/shell cannot be used with 3,4-body interactions"); + if (atom->molecular == 2) + error->all(FLERR, "Bond style rheo/shell cannot be used with atom style template"); +} + +/* ---------------------------------------------------------------------- */ + +void BondRHEOShell::settings(int narg, char **arg) +{ + BondBPM::settings(narg, arg); + + int iarg; + for (std::size_t i = 0; i < leftover_iarg.size(); i++) { + iarg = leftover_iarg[i]; + if (strcmp(arg[iarg], "t/form") == 0) { + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond rheo/shell command, missing option for t/form"); + tform = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (tform < 0.0) error->all(FLERR, "Illegal bond rheo/shell value for t/form, {}", tform); + i += 1; + } else if (strcmp(arg[iarg], "r/max") == 0) { + if (iarg + 1 > narg) error->all(FLERR, "Illegal bond rheo/shell command, missing option for r/max"); + rmax = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (rmax < 0.0) error->all(FLERR, "Illegal bond rheo/shell value for r/max, {}", rmax); + i += 1; + } else { + error->all(FLERR, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]); + } + } + + if (tform < 0.0) + error->all(FLERR, "Illegal bond rheo/shell command, must specify t/form"); + if (rmax < 0.0) + error->all(FLERR, "Illegal bond rheo/shell command, must specify r/max"); +} + + +/* ---------------------------------------------------------------------- + used to check bond communiction cutoff - not perfect, estimates based on local-local only +------------------------------------------------------------------------- */ + +double BondRHEOShell::equilibrium_distance(int /*i*/) +{ + // Divide out heuristic prefactor added in comm class + return max_stretch * rmax / 1.5; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondRHEOShell::write_restart(FILE *fp) +{ + BondBPM::write_restart(fp); + write_restart_settings(fp); + + fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&ecrit[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&gamma[1], sizeof(double), atom->nbondtypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondRHEOShell::read_restart(FILE *fp) +{ + BondBPM::read_restart(fp); + read_restart_settings(fp); + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &ecrit[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &gamma[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + } + MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&ecrit[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&gamma[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file + ------------------------------------------------------------------------- */ + +void BondRHEOShell::write_restart_settings(FILE *fp) +{ + fwrite(&tform, sizeof(double), 1, fp); + fwrite(&rmax, sizeof(double), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts + ------------------------------------------------------------------------- */ + +void BondRHEOShell::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &rmax, sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&tform, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, world); +} + +/* ---------------------------------------------------------------------- */ + +double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce) +{ + if (type <= 0) return 0.0; + + double r0, t; + for (int n = 0; n < atom->num_bond[i]; n++) { + if (atom->bond_atom[i][n] == atom->tag[j]) { + r0 = fix_bond_history->get_atom_value(i, n, 0); + t = fix_bond_history->get_atom_value(i, n, 1); + } + } + + svector[1] = t; + if (t < tform) return 0.0; + + double r = sqrt(rsq); + double rinv = 1.0 / r; + double dr = r0 - r; + fforce = 2 * k[type] * (dr + dr * dr * dr / (r0 * r0 * ecrit[type] * ecrit[type])); + + double **x = atom->x; + double **v = atom->v; + double delx = x[i][0] - x[j][0]; + double dely = x[i][1] - x[j][1]; + double delz = x[i][2] - x[j][2]; + double delvx = v[i][0] - v[j][0]; + double delvy = v[i][1] - v[j][1]; + double delvz = v[i][2] - v[j][2]; + double dot = delx * delvx + dely * delvy + delz * delvz; + fforce -= gamma[type] * dot * rinv; + fforce *= rinv; + + // set single_extra quantities + + svector[0] = r0; + + return 0.0; +} + +/* ---------------------------------------------------------------------- + Similar to BondBPM->process_broken(), but don't send to FixStoreLocal + ------------------------------------------------------------------------- */ + +void BondRHEOShell::process_ineligibility(int i, int j) +{ + // Manually search and remove from atom arrays + int m, n; + int nlocal = atom->nlocal; + + tagint *tag = atom->tag; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + + if (i < nlocal) { + for (m = 0; m < num_bond[i]; m++) { + if (bond_atom[i][m] == tag[j]) { + bond_type[i][m] = 0; + n = num_bond[i]; + bond_type[i][m] = bond_type[i][n - 1]; + bond_atom[i][m] = bond_atom[i][n - 1]; + fix_bond_history->shift_history(i, m, n - 1); + fix_bond_history->delete_history(i, n - 1); + num_bond[i]--; + break; + } + } + } + + if (j < nlocal) { + for (m = 0; m < num_bond[j]; m++) { + if (bond_atom[j][m] == tag[i]) { + bond_type[j][m] = 0; + n = num_bond[j]; + bond_type[j][m] = bond_type[j][n - 1]; + bond_atom[j][m] = bond_atom[j][n - 1]; + fix_bond_history->shift_history(j, m, n - 1); + fix_bond_history->delete_history(j, n - 1); + num_bond[j]--; + break; + } + } + } +} diff --git a/src/RHEO/bond_rheo_shell.h b/src/RHEO/bond_rheo_shell.h new file mode 100644 index 0000000000..562be6d9a6 --- /dev/null +++ b/src/RHEO/bond_rheo_shell.h @@ -0,0 +1,55 @@ +/* -*- 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 BOND_CLASS +// clang-format off +BondStyle(rheo/shell,BondRHEOShell); +// clang-format on +#else + +#ifndef LMP_BOND_RHEO_SHELL_H +#define LMP_BOND_RHEO_SHELL_H + +#include "bond_bpm.h" + +namespace LAMMPS_NS { + +class BondRHEOShell : public BondBPM { + public: + BondRHEOShell(class LAMMPS *); + ~BondRHEOShell() override; + void compute(int, int) override; + void coeff(int, char **) override; + void init_style() override; + void settings(int, char **) override; + double equilibrium_distance(int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_restart_settings(FILE *) override; + void read_restart_settings(FILE *) override; + double single(int, double, int, int, double &) override; + + protected: + double *k, *ecrit, *gamma; + double tform, rmax; + + void process_ineligibility(int, int); + void allocate(); + void store_data(); + double store_bond(int, int, int); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 3cb2fcf058..31930d655d 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -329,7 +329,8 @@ double ComputeRHEOInterface::correct_rho(int i, int j) { // i is wall, j is fluid //In future may depend on atom type j's pressure equation - return atom->rho[i]; + int itype = atom->type[i]; + return MAX(rho0[itype], atom->rho[i]); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp new file mode 100644 index 0000000000..bd7babedbb --- /dev/null +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -0,0 +1,190 @@ +/* ---------------------------------------------------------------------- + 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. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL) +----------------------------------------------------------------------- */ + +#include "fix_rheo_oxidation.h" + +#include "atom.h" +#include "atom_vec.h" +#include "error.h" +#include "fix_rheo.h" +#include "force.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" + +using namespace LAMMPS_NS; +using namespace RHEO_NS; +using namespace FixConst; +enum {NONE, CONSTANT}; + +/* ---------------------------------------------------------------------- */ + +FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) //, fix_bond_history(nullptr) +{ + if (narg != 5) 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); + + 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); + + cutsq = cut * cut; +} + +/* ---------------------------------------------------------------------- */ + +FixRHEOOxidation::~FixRHEOOxidation() +{ +} + +/* ---------------------------------------------------------------------- */ + +int FixRHEOOxidation::setmask() +{ + int mask = 0; + mask |= POST_INTEGRATE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +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"); + class FixRHEO *fix_rheo = dynamic_cast(fixes[0]); + double cut_kernel = fix_rheo->h; + + if (cut > cut_kernel) + error->all(FLERR, "Bonding length exceeds 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"); + + // need a half neighbor list + auto req = neighbor->add_request(this, NeighConst::REQ_DEFAULT); + req->set_cutoff(cut); +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOOxidation::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOOxidation::post_integrate() +{ + int i, j, n, ii, jj, inum, jnum, bflag; + int *ilist, *jlist, *numneigh, **firstneigh; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; + tagint tagi, tagj; + + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + 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 **x = atom->x; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(status[i] & STATUS_SURFACE)) continue; + + tagi = tag[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(status[j] & STATUS_SURFACE)) continue; + + tagj = tag[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + if (rsq > cutsq) continue; + + // Check if already have an oxide bond + bflag = 0; + for (n = 0; n < num_bond[i]; n++) { + if (bond_type[i][n] == btype && bond_atom[i][n] == tagj) { + bflag = 1; + break; + } + } + if (bflag) continue; + + for (n = 0; n < num_bond[j]; n++) { + if (bond_type[j][n] == btype && bond_atom[j][n] == tagi) { + bflag = 1; + break; + } + } + if (bflag) continue; + + // Add bonds to owned atoms + // If newton bond, add to both, otherwise add to whichever has a smaller tag + if (i < nlocal && (!newton_bond || tagi < tagj)) { + if (num_bond[i] == atom->bond_per_atom) + error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/oxidation for atom {}", tagi); + bond_type[i][num_bond[i]] = btype; + bond_atom[i][num_bond[i]] = tagj; + num_bond[i]++; + } + + if (j < nlocal && (!newton_bond || tagj < tagi)) { + if (num_bond[j] == atom->bond_per_atom) + error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/oxidation for atom {}", tagj); + bond_type[j][num_bond[j]] = btype; + bond_atom[j][num_bond[j]] = tagi; + num_bond[j]++; + } + } + } +} diff --git a/src/RHEO/fix_rheo_oxidation.h b/src/RHEO/fix_rheo_oxidation.h new file mode 100644 index 0000000000..4c81605611 --- /dev/null +++ b/src/RHEO/fix_rheo_oxidation.h @@ -0,0 +1,49 @@ +/* -*- 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(rheo/oxidation,FixRHEOOxidation) +// clang-format on +#else + +#ifndef LMP_FIX_RHEO_OXIDATION_H +#define LMP_FIX_RHEO_OXIDATION_H + +#include "fix.h" + +#include + +namespace LAMMPS_NS { + +class FixRHEOOxidation : public Fix { + public: + FixRHEOOxidation(class LAMMPS *, int, char **); + ~FixRHEOOxidation() override; + int setmask() override; + void init() override; + void init_list(int, class NeighList *) override; + void post_integrate() override; + + private: + int btype; + double cut, cutsq; + class NeighList *list; + + //class FixBondHistory *fix_bond_history; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 0640dd6827..c6d10bcc79 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -271,7 +271,12 @@ void FixRHEOThermal::init() req->set_cutoff(cut_kernel); // find instances of bond history to delete data + // skip history for shell, only exception histories = modify->get_fix_by_style("BOND_HISTORY"); + if (n_histories > 0) + for (int i = 0; i < histories.size(); i++) + if (strcmp(histories[i]->id, "HISTORY_RHEO_SHELL") == 0) + histories.erase(histories.begin() + i); n_histories = histories.size(); } } @@ -644,6 +649,8 @@ double FixRHEOThermal::calc_cv(int i, int itype) if (cv_style[itype] == CONSTANT) { return cv[itype]; } + + return 0.0; } /* ---------------------------------------------------------------------- */ @@ -653,6 +660,8 @@ double FixRHEOThermal::calc_Tc(int i, int itype) if (Tc_style[itype] == CONSTANT) { return Tc[itype]; } + + return 0.0; } /* ---------------------------------------------------------------------- */ @@ -662,6 +671,8 @@ double FixRHEOThermal::calc_L(int i, int itype) if (L_style[itype] == CONSTANT) { return L[itype]; } + + return 0.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/bond_hybrid.cpp b/src/bond_hybrid.cpp index 4e477ab3a6..5f84db1886 100644 --- a/src/bond_hybrid.cpp +++ b/src/bond_hybrid.cpp @@ -259,8 +259,13 @@ void BondHybrid::flags() if (styles[m]) comm_forward = MAX(comm_forward, styles[m]->comm_forward); if (styles[m]) comm_reverse = MAX(comm_reverse, styles[m]->comm_reverse); if (styles[m]) comm_reverse_off = MAX(comm_reverse_off, styles[m]->comm_reverse_off); + if (styles[m]) partial_flag = MAX(partial_flag, styles[m]->partial_flag); } + for (m = 0; m < nstyles; m++) + if (styles[m]->partial_flag != partial_flag) + error->all(FLERR, "Cannot hybridize bond styles with different topology settings"); + init_svector(); }