ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Prototyping oxidation

Browse Source
This commit is contained in:
jtclemm
2024-03-28 11:31:21 -06:00
parent b3de75da97
commit 010a4c076b
14 changed files with 1137 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_bond.rst
View File

@ -54,6 +54,7 @@ OPT.
* :doc:`oxdna2/fene <bond_oxdna>`
* :doc:`oxrna2/fene <bond_oxdna>`
* :doc:`quartic (o) <bond_quartic>`
* :doc:`rheo/shell <bond_rheo_shell>`
* :doc:`special <bond_special>`
* :doc:`table (o) <bond_table>`

1
doc/src/Commands_compute.rst
View File

@ -122,6 +122,7 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`reduce <compute_reduce>`
* :doc:`reduce/chunk <compute_reduce_chunk>`
* :doc:`reduce/region <compute_reduce>`
* :doc:`rheo/property/atom <compute_rheo_property_atom>`
* :doc:`rigid/local <compute_rigid_local>`
* :doc:`saed <compute_saed>`
* :doc:`slcsa/atom <compute_slcsa_atom>`

5
doc/src/Commands_fix.rst
View File

@ -203,6 +203,11 @@ OPT.
* :doc:`reaxff/species (k) <fix_reaxff_species>`
* :doc:`recenter <fix_recenter>`
* :doc:`restrain <fix_restrain>`
* :doc:`rheo <fix_rheo>`
* :doc:`rheo/oxidation <fix_rheo_oxidation>`
* :doc:`rheo/pressure <fix_rheo_pressure>`
* :doc:`rheo/thermal <fix_rheo_thermal>`
* :doc:`rheo/viscosity <fix_rheo_viscosity>`
* :doc:`rhok <fix_rhok>`
* :doc:`rigid (o) <fix_rigid>`
* :doc:`rigid/meso <fix_rigid_meso>`

2
doc/src/Commands_pair.rst
View File

@ -257,6 +257,8 @@ OPT.
* :doc:`reaxff (ko) <pair_reaxff>`
* :doc:`rebo (io) <pair_airebo>`
* :doc:`resquared (go) <pair_resquared>`
* :doc:`rheo <pair_rheo>`
* :doc:`rheo/solid <pair_rheo_solid>`
* :doc:`saip/metal (t) <pair_saip_metal>`
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
* :doc:`smatb <pair_smatb>`

226
doc/src/bond_rheo_shell.rst Normal file
View File

@ -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 <restart>` 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) <fragment-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) <Groot4>`:
.. 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 <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <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 <special_bonds>` settings described in the
restrictions. Further details can be found in the :doc:`how to <Howto_bpm>`
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 <dump>` 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 <dump>` command, bonds with type 0 (broken bonds)
are not included.
The :doc:`delete_bonds <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 <restart>`. 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 <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 <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
<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 <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 <bond_coeff>`, :doc:`pair bpm/spring <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).

81
doc/src/fix_rheo_oxidation.rst Normal file
View File

@ -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 <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 <restart>`. None of the :doc:`fix_modify <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 <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
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 <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 <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
none

4
src/BPM/bond_bpm.cpp
View File

@ -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));
}

507
src/RHEO/bond_rheo_shell.cpp Normal file
View File

@ -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 <cmath>
#include <cstring>
#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<FixBondHistory *>(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<FixStoreLocal *>(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;
}
}
}
}

55
src/RHEO/bond_rheo_shell.h Normal file
View File

@ -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

3
src/RHEO/compute_rheo_interface.cpp
View File

@ -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]);
}
/* ---------------------------------------------------------------------- */

190
src/RHEO/fix_rheo_oxidation.cpp Normal file
View File

@ -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<FixRHEO *>(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<FixBondHistory *>(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]++;
}
}
}
}

49
src/RHEO/fix_rheo_oxidation.h Normal file
View File

@ -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 <vector>
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

11
src/RHEO/fix_rheo_thermal.cpp
View File

@ -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;
}
/* ---------------------------------------------------------------------- */

5
src/bond_hybrid.cpp
View File

@ -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();
}