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

Merge branch 'develop' into general-triclinic

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-03-17 15:56:38 -04:00
parent 01d412bdbf 207a14f351
commit 1bdcf7e42e
15 changed files with 684 additions and 279 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Intro_citing.rst
View File

@ -47,6 +47,8 @@ In addition there are DOIs generated for individual stable releases:
- 3 March 2020 version: `DOI:10.5281/zenodo.3726417 <https://dx.doi.org/10.5281/zenodo.3726417>`_ - 3 March 2020 version: `DOI:10.5281/zenodo.3726417 <https://dx.doi.org/10.5281/zenodo.3726417>`_
- 29 October 2020 version: `DOI:10.5281/zenodo.4157471 <https://dx.doi.org/10.5281/zenodo.4157471>`_ - 29 October 2020 version: `DOI:10.5281/zenodo.4157471 <https://dx.doi.org/10.5281/zenodo.4157471>`_
- 29 September 2021 version: `DOI:10.5281/zenodo.6386596 <https://dx.doi.org/10.5281/zenodo.6386596>`_ - 29 September 2021 version: `DOI:10.5281/zenodo.6386596 <https://dx.doi.org/10.5281/zenodo.6386596>`_
- 23 June 2022 version: `DOI:10.5281/zenodo.10806836 <https://doi.org/10.5281/zenodo.10806836>`_
- 2 August 2023 version: `DOI:10.5281/zenodo.10806852 <https://doi.org/10.5281/zenodo.10806852>`_
Home page Home page
^^^^^^^^^ ^^^^^^^^^

2
doc/src/fix_deform_pressure.rst
View File

@ -91,7 +91,7 @@ corresponding component of the pressure tensor. This option attempts to
maintain a specified target pressure using a linear controller where the maintain a specified target pressure using a linear controller where the
box length :math:`L` evolves according to the equation box length :math:`L` evolves according to the equation
.. parsed-literal:: .. math::
\frac{d L(t)}{dt} = L(t) k (P_t - P) \frac{d L(t)}{dt} = L(t) k (P_t - P)

90
doc/src/fix_indent.rst
View File

@ -8,33 +8,44 @@ Syntax
.. code-block:: LAMMPS .. code-block:: LAMMPS
fix ID group-ID indent K keyword values ... fix ID group-ID indent K gstyle args keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command * ID, group-ID are documented in :doc:`fix <fix>` command
* indent = style name of this fix command * indent = style name of this fix command
* K = force constant for indenter surface (force/distance\^2 units) * K = force constant for indenter surface (force/distance\^2 units)
* one or more keyword/value pairs may be appended * gstyle = *sphere* or *cylinder* or *cone* or *plane*
* keyword = *sphere* or *cylinder* or *plane* or *side* or *units*
.. parsed-literal:: .. parsed-literal::
*sphere* args = x y z R *sphere* args = x y z R
x,y,z = position of center of indenter (distance units) x, y, z = position of center of indenter (distance units)
R = sphere radius of indenter (distance units) R = sphere radius of indenter (distance units)
any of x,y,z,R can be a variable (see below) any of x, y, z, R can be a variable (see below)
*cylinder* args = dim c1 c2 R *cylinder* args = dim c1 c2 R
dim = *x* or *y* or *z* = axis of cylinder dim = *x* or *y* or *z* = axis of cylinder
c1,c2 = coords of cylinder axis in other 2 dimensions (distance units) c1, c2 = coords of cylinder axis in other 2 dimensions (distance units)
R = cylinder radius of indenter (distance units) R = cylinder radius of indenter (distance units)
any of c1,c2,R can be a variable (see below) any of c1,c2,R can be a variable (see below)
*cone* args = dim c1 c2 radlo radhi lo hi
dim = *x* or *y* or *z* = axis of cone
c1, c2 = coords of cone axis in other 2 dimensions (distance units)
radlo,radhi = cone radii at lo and hi end (distance units)
lo,hi = bounds of cone in dim (distance units)
any of c1, c2, radlo, radhi, lo, hi can be a variable (see below)
*plane* args = dim pos side *plane* args = dim pos side
dim = *x* or *y* or *z* = plane perpendicular to this dimension dim = *x* or *y* or *z* = plane perpendicular to this dimension
pos = position of plane in dimension x, y, or z (distance units) pos = position of plane in dimension x, y, or z (distance units)
pos can be a variable (see below) pos can be a variable (see below)
side = *lo* or *hi* side = *lo* or *hi*
* zero or more keyword/value pairs may be appended
* keyword = *side* or *units*
.. parsed-literal::
*side* value = *in* or *out* *side* value = *in* or *out*
*in* = the indenter acts on particles inside the sphere or cylinder *in* = the indenter acts on particles inside the sphere or cylinder or cone
*out* = the indenter acts on particles outside the sphere or cylinder *out* = the indenter acts on particles outside the sphere or cylinder or cone
*units* value = *lattice* or *box* *units* value = *lattice* or *box*
lattice = the geometry is defined in lattice units lattice = the geometry is defined in lattice units
box = the geometry is defined in simulation box units box = the geometry is defined in simulation box units
@ -53,12 +64,12 @@ Description
Insert an indenter within a simulation box. The indenter repels all Insert an indenter within a simulation box. The indenter repels all
atoms in the group that touch it, so it can be used to push into a atoms in the group that touch it, so it can be used to push into a
material or as an obstacle in a flow. Or it can be used as a material or as an obstacle in a flow. Alternatively, it can be used as a
constraining wall around a simulation; see the discussion of the constraining wall around a simulation; see the discussion of the
*side* keyword below. *side* keyword below.
The indenter can either be spherical or cylindrical or planar. You The *gstyle* geometry of the indenter can either be a sphere, a
must set one of those 3 keywords. cylinder, a cone, or a plane.
A spherical indenter exerts a force of magnitude A spherical indenter exerts a force of magnitude
@ -75,15 +86,20 @@ A cylindrical indenter exerts the same force, except that *r* is the
distance from the atom to the center axis of the cylinder. The distance from the atom to the center axis of the cylinder. The
cylinder extends infinitely along its axis. cylinder extends infinitely along its axis.
Spherical and cylindrical indenters account for periodic boundaries in A conical indenter is similar to a cylindrical indenter except that it
two ways. First, the center point of a spherical indenter (x,y,z) or has a finite length (between *lo* and *hi*), and that two different
axis of a cylindrical indenter (c1,c2) is remapped back into the radii (one at each end, *radlo* and *radhi*) can be defined.
simulation box, if the box is periodic in a particular dimension.
This occurs every timestep if the indenter geometry is specified with Spherical, cylindrical, and conical indenters account for periodic
a variable (see below), e.g. it is moving over time. Second, the boundaries in two ways. First, the center point of a spherical
calculation of distance to the indenter center or axis accounts for indenter (x,y,z) or axis of a cylindrical/conical indenter (c1,c2) is
periodic boundaries. Both of these mean that an indenter can remapped back into the simulation box, if the box is periodic in a
effectively move through and straddle one or more periodic boundaries. particular dimension. This occurs every timestep if the indenter
geometry is specified with a variable (see below), e.g. it is moving
over time. Second, the calculation of distance to the indenter center
or axis accounts for periodic boundaries. Both of these mean that an
indenter can effectively move through and straddle one or more
periodic boundaries.
A planar indenter is really an axis-aligned infinite-extent wall A planar indenter is really an axis-aligned infinite-extent wall
exerting the same force on atoms in the system, where *R* is the exerting the same force on atoms in the system, where *R* is the
@ -97,9 +113,13 @@ is specified as *hi*\ .
Any of the 4 quantities defining a spherical indenter's geometry can Any of the 4 quantities defining a spherical indenter's geometry can
be specified as an equal-style :doc:`variable <variable>`, namely *x*, be specified as an equal-style :doc:`variable <variable>`, namely *x*,
*y*, *z*, or *R*\ . Similarly, for a cylindrical indenter, any of *c1*, *y*, *z*, or *R*\ . For a cylindrical indenter, any of the 3
*c2*, or *R*, can be a variable. For a planar indenter, *pos* can be quantities *c1*, *c2*, or *R*, can be a variable. For a conical
a variable. If the value is a variable, it should be specified as indenter, any of the 6 quantities *c1*, *c2*, *radlo*, *radhi*, *lo*,
or *hi* can be a variable. For a planar indenter, the single value
*pos* can be a variable.
If any of these values is a variable, it should be specified as
v_name, where name is the variable name. In this case, the variable v_name, where name is the variable name. In this case, the variable
will be evaluated each timestep, and its value used to define the will be evaluated each timestep, and its value used to define the
indenter geometry. indenter geometry.
@ -110,7 +130,8 @@ command keywords for the simulation box parameters and timestep and
elapsed time. Thus it is easy to specify indenter properties that elapsed time. Thus it is easy to specify indenter properties that
change as a function of time or span consecutive runs in a continuous change as a function of time or span consecutive runs in a continuous
fashion. For the latter, see the *start* and *stop* keywords of the fashion. For the latter, see the *start* and *stop* keywords of the
:doc:`run <run>` command and the *elaplong* keyword of :doc:`thermo_style custom <thermo_style>` for details. :doc:`run <run>` command and the *elaplong* keyword of
:doc:`thermo_style custom <thermo_style>` for details.
For example, if a spherical indenter's x-position is specified as v_x, For example, if a spherical indenter's x-position is specified as v_x,
then this variable definition will keep it's center at a relative then this variable definition will keep it's center at a relative
@ -141,12 +162,13 @@ rate.
If the *side* keyword is specified as *out*, which is the default, If the *side* keyword is specified as *out*, which is the default,
then particles outside the indenter are pushed away from its outer then particles outside the indenter are pushed away from its outer
surface, as described above. This only applies to spherical or surface, as described above. This only applies to spherical,
cylindrical indenters. If the *side* keyword is specified as *in*, cylindrical, and conical indenters. If the *side* keyword is
the action of the indenter is reversed. Particles inside the indenter specified as *in*, the action of the indenter is reversed. Particles
are pushed away from its inner surface. In other words, the indenter inside the indenter are pushed away from its inner surface. In other
is now a containing wall that traps the particles inside it. If the words, the indenter is now a containing wall that traps the particles
radius shrinks over time, it will squeeze the particles. inside it. If the radius shrinks over time, it will squeeze the
particles.
The *units* keyword determines the meaning of the distance units used The *units* keyword determines the meaning of the distance units used
to define the indenter geometry. A *box* value selects standard to define the indenter geometry. A *box* value selects standard
@ -166,10 +188,10 @@ lattice spacings in a variable formula.
The force constant *K* is not affected by the *units* keyword. It is The force constant *K* is not affected by the *units* keyword. It is
always in force/distance\^2 units where force and distance are defined always in force/distance\^2 units where force and distance are defined
by the :doc:`units <units>` command. If you wish K to be scaled by the by the :doc:`units <units>` command. If you wish K to be scaled by
lattice spacing, you can define K with a variable whose formula the lattice spacing, you can define K with a variable whose formula
contains *xlat*, *ylat*, *zlat* keywords of the contains *xlat*, *ylat*, *zlat* keywords of the :doc:`thermo_style
:doc:`thermo_style <thermo_style>` command, e.g. <thermo_style>` command, e.g.
.. code-block:: LAMMPS .. code-block:: LAMMPS

8
src/ELECTRODE/fix_electrode_conp.cpp
View File

@ -915,20 +915,20 @@ void FixElectrodeConp::update_charges()
a = ele_ele_interaction(q_local); a = ele_ele_interaction(q_local);
r = add_nlocalele(b, a); r = add_nlocalele(b, a);
} else { } else {
r = add_nlocalele(r, scale_vector(alpha, y)); r = add_nlocalele(r, scale_vector(alpha, std::move(y)));
} }
auto p = constraint_projection(r); auto p = constraint_projection(r);
double dot_new = dot_nlocalele(r, p); double dot_new = dot_nlocalele(r, p);
d = add_nlocalele(p, scale_vector(dot_new / dot_old, d)); d = add_nlocalele(std::move(p), scale_vector(dot_new / dot_old, d));
delta = dot_nlocalele(r, d); delta = dot_nlocalele(r, d);
dot_old = dot_new; dot_old = dot_new;
} }
recompute_potential(b, q_local); recompute_potential(std::move(b), q_local);
if (delta > cg_threshold && comm->me == 0) error->warning(FLERR, "CG threshold not reached"); if (delta > cg_threshold && comm->me == 0) error->warning(FLERR, "CG threshold not reached");
} else { } else {
error->all(FLERR, "This algorithm is not implemented, yet"); error->all(FLERR, "This algorithm is not implemented, yet");
} }
set_charges(q_local); set_charges(std::move(q_local));
update_time += MPI_Wtime() - start; update_time += MPI_Wtime() - start;
} }

4
src/EXTRA-FIX/fix_deform_pressure.cpp
View File

@ -767,8 +767,7 @@ void FixDeformPressure::apply_box()
if (fabs(v_rate) > max_h_rate) if (fabs(v_rate) > max_h_rate)
v_rate = max_h_rate * v_rate / fabs(v_rate); v_rate = max_h_rate * v_rate / fabs(v_rate);
set_extra[6].cumulative_strain += update->dt * v_rate; scale = (1.0 + update->dt * v_rate);
scale = (1.0 + set_extra[6].cumulative_strain);
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale; shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift; set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
@ -843,7 +842,6 @@ void FixDeformPressure::restart(char *buf)
set_extra[i].saved = set_extra_restart[i].saved; set_extra[i].saved = set_extra_restart[i].saved;
set_extra[i].prior_rate = set_extra_restart[i].prior_rate; set_extra[i].prior_rate = set_extra_restart[i].prior_rate;
set_extra[i].prior_pressure = set_extra_restart[i].prior_pressure; set_extra[i].prior_pressure = set_extra_restart[i].prior_pressure;
set_extra[i].cumulative_strain = set_extra_restart[i].cumulative_strain;
} }
} }

1
src/EXTRA-FIX/fix_deform_pressure.h
View File

@ -51,7 +51,6 @@ class FixDeformPressure : public FixDeform {
struct SetExtra { struct SetExtra {
double ptarget, pgain; double ptarget, pgain;
double prior_pressure, prior_rate; double prior_pressure, prior_rate;
double cumulative_strain;
int saved; int saved;
char *pstr; char *pstr;
int pvar, pvar_flag; int pvar, pvar_flag;

6
src/EXTRA-FIX/fix_nonaffine_displacement.cpp
View File

@ -65,7 +65,8 @@ static const char cite_nonaffine_d2min[] =
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr), D2min(nullptr) Fix(lmp, narg, arg), id_fix(nullptr), fix(nullptr), D2min(nullptr), X(nullptr), Y(nullptr),
F(nullptr), norm(nullptr)
{ {
if (narg < 4) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error); if (narg < 4) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error);
@ -85,7 +86,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char *
} else if (strcmp(arg[iarg + 1], "radius") == 0) { } else if (strcmp(arg[iarg + 1], "radius") == 0) {
cut_style = RADIUS; cut_style = RADIUS;
} else if (strcmp(arg[iarg + 1], "custom") == 0) { } else if (strcmp(arg[iarg + 1], "custom") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement custom", error); if (iarg + 2 > narg)
utils::missing_cmd_args(FLERR,"fix nonaffine/displacement custom", error);
if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD)) if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD))
error->all(FLERR, "Fix nonaffine/displacement with custom cutoff requires neighbor style 'bin' or 'nsq'"); error->all(FLERR, "Fix nonaffine/displacement with custom cutoff requires neighbor style 'bin' or 'nsq'");
cut_style = CUSTOM; cut_style = CUSTOM;

1
src/EXTRA-FIX/fix_nonaffine_displacement.h
View File

@ -57,7 +57,6 @@ class FixNonaffineDisplacement : public Fix {
class NeighList *list; // half neighbor list class NeighList *list; // half neighbor list
void integrate_velocity(); void integrate_velocity();
void calculate_D2Min(); void calculate_D2Min();
void save_reference_state(); void save_reference_state();

5
src/GRANULAR/granular_model.cpp
View File

@ -32,6 +32,7 @@
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
#include <utility>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace Granular_NS; using namespace Granular_NS;
@ -333,11 +334,11 @@ void GranularModel::read_restart(FILE *fp)
utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error); utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&num_char, 1, MPI_INT, 0, world); MPI_Bcast(&num_char, 1, MPI_INT, 0, world);
std::string model_name (num_char, ' '); std::string model_name(num_char, ' ');
if (comm->me == 0) if (comm->me == 0)
utils::sfread(FLERR, const_cast<char*>(model_name.data()), sizeof(char),num_char, fp, nullptr, error); utils::sfread(FLERR, const_cast<char*>(model_name.data()), sizeof(char),num_char, fp, nullptr, error);
MPI_Bcast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, 0, world); MPI_Bcast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, 0, world);
construct_sub_model(model_name, (SubModelType) i); construct_sub_model(std::move(model_name), (SubModelType) i);
if (comm->me == 0) if (comm->me == 0)
utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error); utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error);

2
src/KOKKOS/pair_kokkos.h
View File

@ -1020,7 +1020,7 @@ template<class PairStyle, class Specialisation = void>
EV_FLOAT pair_compute (PairStyle* fpair, NeighListKokkos<typename PairStyle::device_type>* list) { EV_FLOAT pair_compute (PairStyle* fpair, NeighListKokkos<typename PairStyle::device_type>* list) {
EV_FLOAT ev; EV_FLOAT ev;
if (fpair->neighflag == FULL) { if (fpair->neighflag == FULL) {
if (utils::strmatch(fpair->lmp->force->pair_style,"^hybrid/overlay")) { if (utils::strmatch(fpair->lmp->force->pair_style,"^hybrid")) {
fpair->fuse_force_clear_flag = 0; fpair->fuse_force_clear_flag = 0;
ev = pair_compute_neighlist<PairStyle,FULL,0,Specialisation> (fpair,list); ev = pair_compute_neighlist<PairStyle,FULL,0,Specialisation> (fpair,list);
} else { } else {

4
src/REAXFF/fix_qeq_reaxff.cpp
View File

@ -141,7 +141,7 @@ FixQEqReaxFF::FixQEqReaxFF(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays // perform initial allocation of atom-based arrays
// register with Atom class // register with Atom class
reaxff = dynamic_cast<PairReaxFF *>(force->pair_match("^reax..",0)); reaxff = dynamic_cast<PairReaxFF *>(force->pair_match("^reaxff",0));
s_hist = t_hist = nullptr; s_hist = t_hist = nullptr;
atom->add_callback(Atom::GROW); atom->add_callback(Atom::GROW);
@ -217,6 +217,8 @@ void FixQEqReaxFF::pertype_parameters(char *arg)
if (chi == nullptr || eta == nullptr || gamma == nullptr) if (chi == nullptr || eta == nullptr || gamma == nullptr)
error->all(FLERR, "Fix qeq/reaxff could not extract params from pair reaxff"); error->all(FLERR, "Fix qeq/reaxff could not extract params from pair reaxff");
return; return;
} else if (utils::strmatch(arg,"^reax/c")) {
error->all(FLERR, "Fix qeq/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'");
} }
reaxflag = 0; reaxflag = 0;

806
src/fix_indent.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -23,6 +22,7 @@
#include "error.h" #include "error.h"
#include "input.h" #include "input.h"
#include "lattice.h" #include "lattice.h"
#include "math_extra.h"
#include "modify.h" #include "modify.h"
#include "respa.h" #include "respa.h"
#include "update.h" #include "update.h"
@ -34,14 +34,14 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{NONE,SPHERE,CYLINDER,PLANE}; enum { NONE, SPHERE, CYLINDER, PLANE, CONE };
enum{INSIDE,OUTSIDE}; enum { INSIDE, OUTSIDE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr),
xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr) rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr)
{ {
if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error); if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error);
@ -55,22 +55,20 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
respa_level_support = 1; respa_level_support = 1;
ilevel_respa = 0; ilevel_respa = 0;
k = utils::numeric(FLERR,arg[3],false,lmp); k = utils::numeric(FLERR, arg[3], false, lmp);
k3 = k/3.0; if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k);
k3 = k / 3.0;
// read options from end of input line // read geometry of indenter and optional args
options(narg-4,&arg[4]); int iarg = geometry(narg - 4, &arg[4]) + 4;
options(narg - iarg, &arg[iarg]);
// setup scaling // setup scaling
double xscale,yscale,zscale; const double xscale{scaleflag ? domain->lattice->xlattice : 1.0};
if (scaleflag) { const double yscale{scaleflag ? domain->lattice->ylattice : 1.0};
xscale = domain->lattice->xlattice; const double zscale{scaleflag ? domain->lattice->zlattice : 1.0};
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
// apply scaling factors to geometry // apply scaling factors to geometry
@ -79,14 +77,43 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
if (!ystr) yvalue *= yscale; if (!ystr) yvalue *= yscale;
if (!zstr) zvalue *= zscale; if (!zstr) zvalue *= zscale;
if (!rstr) rvalue *= xscale; if (!rstr) rvalue *= xscale;
} else if (istyle == CONE) {
if (!xstr) xvalue *= xscale;
if (!ystr) yvalue *= yscale;
if (!zstr) zvalue *= zscale;
double scaling_factor = 1.0;
switch (cdim) {
case 0:
scaling_factor = xscale;
break;
case 1:
scaling_factor = yscale;
break;
case 2:
scaling_factor = zscale;
break;
}
if (!rlostr) rlovalue *= scaling_factor;
if (!rhistr) rhivalue *= scaling_factor;
if (!lostr) lovalue *= scaling_factor;
if (!histr) hivalue *= scaling_factor;
} else if (istyle == PLANE) { } else if (istyle == PLANE) {
if (cdim == 0 && !pstr) pvalue *= xscale; if (cdim == 0 && !pstr)
else if (cdim == 1 && !pstr) pvalue *= yscale; pvalue *= xscale;
else if (cdim == 2 && !pstr) pvalue *= zscale; else if (cdim == 1 && !pstr)
} else error->all(FLERR,"Unknown fix indent keyword: {}", istyle); pvalue *= yscale;
else if (cdim == 2 && !pstr)
pvalue *= zscale;
} else
error->all(FLERR, "Unknown fix indent keyword: {}", istyle);
varflag = 0; varflag = 0;
if (xstr || ystr || zstr || rstr || pstr) varflag = 1; if (xstr || ystr || zstr || rstr || pstr || rlostr || rhistr || lostr || histr) varflag = 1;
indenter_flag = 0; indenter_flag = 0;
indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
@ -96,11 +123,15 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
FixIndent::~FixIndent() FixIndent::~FixIndent()
{ {
delete [] xstr; delete[] xstr;
delete [] ystr; delete[] ystr;
delete [] zstr; delete[] zstr;
delete [] rstr; delete[] rstr;
delete [] pstr; delete[] pstr;
delete[] rlostr;
delete[] rhistr;
delete[] lostr;
delete[] histr;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -120,43 +151,62 @@ void FixIndent::init()
{ {
if (xstr) { if (xstr) {
xvar = input->variable->find(xstr); xvar = input->variable->find(xstr);
if (xvar < 0) if (xvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", xstr);
error->all(FLERR,"Variable {} for fix indent does not exist", xstr);
if (!input->variable->equalstyle(xvar)) if (!input->variable->equalstyle(xvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", xstr); error->all(FLERR, "Variable {} for fix indent is invalid style", xstr);
} }
if (ystr) { if (ystr) {
yvar = input->variable->find(ystr); yvar = input->variable->find(ystr);
if (yvar < 0) if (yvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", ystr);
error->all(FLERR,"Variable {} for fix indent does not exist", ystr);
if (!input->variable->equalstyle(yvar)) if (!input->variable->equalstyle(yvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", ystr); error->all(FLERR, "Variable {} for fix indent is invalid style", ystr);
} }
if (zstr) { if (zstr) {
zvar = input->variable->find(zstr); zvar = input->variable->find(zstr);
if (zvar < 0) if (zvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", zstr);
error->all(FLERR,"Variable {} for fix indent does not exist", zstr);
if (!input->variable->equalstyle(zvar)) if (!input->variable->equalstyle(zvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", zstr); error->all(FLERR, "Variable {} for fix indent is invalid style", zstr);
} }
if (rstr) { if (rstr) {
rvar = input->variable->find(rstr); rvar = input->variable->find(rstr);
if (rvar < 0) if (rvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rstr);
error->all(FLERR,"Variable {} for fix indent does not exist", rstr);
if (!input->variable->equalstyle(rvar)) if (!input->variable->equalstyle(rvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", rstr); error->all(FLERR, "Variable {} for fix indent is invalid style", rstr);
} }
if (pstr) { if (pstr) {
pvar = input->variable->find(pstr); pvar = input->variable->find(pstr);
if (pvar < 0) if (pvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", pstr);
error->all(FLERR,"Variable {} for fix indent does not exist", pstr);
if (!input->variable->equalstyle(pvar)) if (!input->variable->equalstyle(pvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", pstr); error->all(FLERR, "Variable {} for fix indent is invalid style", pstr);
}
if (rlostr) {
rlovar = input->variable->find(rlostr);
if (rlovar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rlostr);
if (!input->variable->equalstyle(rlovar))
error->all(FLERR, "Variable {} for fix indent is invalid style", rlostr);
}
if (rhistr) {
rhivar = input->variable->find(rhistr);
if (rhivar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rhistr);
if (!input->variable->equalstyle(rhivar))
error->all(FLERR, "Variable {} for fix indent is invalid style", rhistr);
}
if (lostr) {
lovar = input->variable->find(lostr);
if (lovar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", lostr);
if (!input->variable->equalstyle(lovar))
error->all(FLERR, "Variable {} for fix indent is invalid style", lostr);
}
if (histr) {
hivar = input->variable->find(histr);
if (hivar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", histr);
if (!input->variable->equalstyle(hivar))
error->all(FLERR, "Variable {} for fix indent is invalid style", histr);
} }
if (utils::strmatch(update->integrate_style,"^respa")) { if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels-1; ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
} }
} }
@ -164,11 +214,11 @@ void FixIndent::init()
void FixIndent::setup(int vflag) void FixIndent::setup(int vflag)
{ {
if (utils::strmatch(update->integrate_style,"^verlet")) if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag); post_force(vflag);
else { else {
(dynamic_cast<Respa *>(update->integrate))->copy_flevel_f(ilevel_respa); (dynamic_cast<Respa *>(update->integrate))->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0); post_force_respa(vflag, ilevel_respa, 0);
(dynamic_cast<Respa *>(update->integrate))->copy_f_flevel(ilevel_respa); (dynamic_cast<Respa *>(update->integrate))->copy_f_flevel(ilevel_respa);
} }
} }
@ -192,55 +242,53 @@ void FixIndent::post_force(int /*vflag*/)
indenter_flag = 0; indenter_flag = 0;
indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
// spherical indenter // ctr = current indenter centerz
if (istyle == SPHERE) { double ctr[3] = {xvalue, yvalue, zvalue};
// ctr = current indenter center
// remap into periodic box
double ctr[3];
if (xstr) ctr[0] = input->variable->compute_equal(xvar); if (xstr) ctr[0] = input->variable->compute_equal(xvar);
else ctr[0] = xvalue;
if (ystr) ctr[1] = input->variable->compute_equal(yvar); if (ystr) ctr[1] = input->variable->compute_equal(yvar);
else ctr[1] = yvalue;
if (zstr) ctr[2] = input->variable->compute_equal(zvar); if (zstr) ctr[2] = input->variable->compute_equal(zvar);
else ctr[2] = zvalue;
domain->remap(ctr);
double radius;
if (rstr) radius = input->variable->compute_equal(rvar);
else radius = rvalue;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double delx,dely,delz,r,dr,fmag,fx,fy,fz; double delx, dely, delz, r, dr, fmag, fx, fy, fz;
// spherical indenter
if (istyle == SPHERE) {
// remap indenter center into periodic box
domain->remap(ctr);
double radius = rstr ? input->variable->compute_equal(rvar) : rvalue;
if (radius < 0.0) error->all(FLERR, "Illegal fix indent sphere radius: {}", radius);
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
delx = x[i][0] - ctr[0]; delx = x[i][0] - ctr[0];
dely = x[i][1] - ctr[1]; dely = x[i][1] - ctr[1];
delz = x[i][2] - ctr[2]; delz = x[i][2] - ctr[2];
domain->minimum_image(delx,dely,delz); domain->minimum_image(delx, dely, delz);
r = sqrt(delx*delx + dely*dely + delz*delz); r = sqrt(delx * delx + dely * dely + delz * delz);
if (side == OUTSIDE) { if (side == OUTSIDE) {
dr = r - radius; dr = r - radius;
fmag = k*dr*dr; fmag = k * dr * dr;
} else { } else {
dr = radius - r; dr = radius - r;
fmag = -k*dr*dr; fmag = -k * dr * dr;
} }
if (dr >= 0.0) continue; if (dr >= 0.0) continue;
fx = delx*fmag/r; fx = delx * fmag / r;
fy = dely*fmag/r; fy = dely * fmag / r;
fz = delz*fmag/r; fz = delz * fmag / r;
f[i][0] += fx; f[i][0] += fx;
f[i][1] += fy; f[i][1] += fy;
f[i][2] += fz; f[i][2] += fz;
indenter[0] -= k3 * dr*dr*dr; indenter[0] -= k3 * dr * dr * dr;
indenter[1] -= fx; indenter[1] -= fx;
indenter[2] -= fy; indenter[2] -= fy;
indenter[3] -= fz; indenter[3] -= fz;
@ -254,101 +302,120 @@ void FixIndent::post_force(int /*vflag*/)
// remap into periodic box // remap into periodic box
// 3rd coord is just near box for remap(), since isn't used // 3rd coord is just near box for remap(), since isn't used
double ctr[3]; ctr[cdim] = domain->boxlo[cdim];
if (cdim == 0) {
ctr[0] = domain->boxlo[0];
if (ystr) ctr[1] = input->variable->compute_equal(yvar);
else ctr[1] = yvalue;
if (zstr) ctr[2] = input->variable->compute_equal(zvar);
else ctr[2] = zvalue;
} else if (cdim == 1) {
if (xstr) ctr[0] = input->variable->compute_equal(xvar);
else ctr[0] = xvalue;
ctr[1] = domain->boxlo[1];
if (zstr) ctr[2] = input->variable->compute_equal(zvar);
else ctr[2] = zvalue;
} else {
if (xstr) ctr[0] = input->variable->compute_equal(xvar);
else ctr[0] = xvalue;
if (ystr) ctr[1] = input->variable->compute_equal(yvar);
else ctr[1] = yvalue;
ctr[2] = domain->boxlo[2];
}
domain->remap(ctr); domain->remap(ctr);
double radius; double radius{rstr ? input->variable->compute_equal(rvar) : rvalue};
if (rstr) radius = input->variable->compute_equal(rvar); if (radius < 0.0) error->all(FLERR, "Illegal fix indent cylinder radius: {}", radius);
else radius = rvalue;
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delx,dely,delz,r,dr,fmag,fx,fy,fz;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (cdim == 0) { double del[3] = {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]};
delx = 0; del[cdim] = 0;
dely = x[i][1] - ctr[1]; domain->minimum_image(del[0], del[1], del[2]);
delz = x[i][2] - ctr[2]; r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]);
} else if (cdim == 1) {
delx = x[i][0] - ctr[0];
dely = 0;
delz = x[i][2] - ctr[2];
} else {
delx = x[i][0] - ctr[0];
dely = x[i][1] - ctr[1];
delz = 0;
}
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
if (side == OUTSIDE) { if (side == OUTSIDE) {
dr = r - radius; dr = r - radius;
fmag = k*dr*dr; fmag = k * dr * dr;
} else { } else {
dr = radius - r; dr = radius - r;
fmag = -k*dr*dr; fmag = -k * dr * dr;
} }
if (dr >= 0.0) continue; if (dr >= 0.0) continue;
fx = delx*fmag/r; fx = del[0] * fmag / r;
fy = dely*fmag/r; fy = del[1] * fmag / r;
fz = delz*fmag/r; fz = del[2] * fmag / r;
f[i][0] += fx; f[i][0] += fx;
f[i][1] += fy; f[i][1] += fy;
f[i][2] += fz; f[i][2] += fz;
indenter[0] -= k3 * dr*dr*dr; indenter[0] -= k3 * dr * dr * dr;
indenter[1] -= fx; indenter[1] -= fx;
indenter[2] -= fy; indenter[2] -= fy;
indenter[3] -= fz; indenter[3] -= fz;
} }
// conical indenter
} else if (istyle == CONE) {
double radiuslo{rlostr ? input->variable->compute_equal(rlovar) : rlovalue};
if (radiuslo < 0.0) error->all(FLERR, "Illegal fix indent cone lower radius: {}", radiuslo);
double radiushi{rhistr ? input->variable->compute_equal(rhivar) : rhivalue};
if (radiushi < 0.0) error->all(FLERR, "Illegal fix indent cone high radius: {}", radiushi);
double initial_lo{lostr ? input->variable->compute_equal(lovar) : lovalue};
double initial_hi{histr ? input->variable->compute_equal(hivar) : hivalue};
ctr[cdim] = 0.5 * (initial_hi + initial_lo);
domain->remap(ctr);
double hi = ctr[cdim] + 0.5 * (initial_hi - initial_lo);
double lo = ctr[cdim] - 0.5 * (initial_hi - initial_lo);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
delx = x[i][0] - ctr[0];
dely = x[i][1] - ctr[1];
delz = x[i][2] - ctr[2];
domain->minimum_image(delx, dely, delz);
double x0[3] = {delx + ctr[0], dely + ctr[1], delz + ctr[2]};
r = sqrt(delx * delx + dely * dely + delz * delz);
// check if particle is inside or outside the cone
bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0);
if (side == INSIDE && point_inside_cone) continue;
if (side == OUTSIDE && !point_inside_cone) continue;
// find the distance between the point and the cone
if (point_inside_cone) {
DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]);
} else {
DistanceExteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]);
}
// compute the force from the center of the cone
// this is different from how it is done in fix wall/region
dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]);
int force_sign = {point_inside_cone ? 1 : -1};
fmag = force_sign * k * dr * dr;
fx = delx * fmag / r;
fy = dely * fmag / r;
fz = delz * fmag / r;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
indenter[0] -= k3 * dr * dr * dr;
indenter[1] -= fx;
indenter[2] -= fy;
indenter[3] -= fz;
}
}
// planar indenter // planar indenter
} else { } else {
// plane = current plane position // plane = current plane position
double plane; double plane{pstr ? input->variable->compute_equal(pvar) : pvalue};
if (pstr) plane = input->variable->compute_equal(pvar);
else plane = pvalue;
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double dr,fatom;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
dr = planeside * (plane - x[i][cdim]); dr = planeside * (plane - x[i][cdim]);
if (dr >= 0.0) continue; if (dr >= 0.0) continue;
fatom = -planeside * k*dr*dr; fmag = -planeside * k * dr * dr;
f[i][cdim] += fatom; f[i][cdim] += fmag;
indenter[0] -= k3 * dr*dr*dr; indenter[0] -= k3 * dr * dr * dr;
indenter[cdim+1] -= fatom; indenter[cdim + 1] -= fmag;
} }
} }
@ -378,7 +445,7 @@ double FixIndent::compute_scalar()
// only sum across procs one time // only sum across procs one time
if (indenter_flag == 0) { if (indenter_flag == 0) {
MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(indenter, indenter_all, 4, MPI_DOUBLE, MPI_SUM, world);
indenter_flag = 1; indenter_flag = 1;
} }
return indenter_all[0]; return indenter_all[0];
@ -393,113 +460,406 @@ double FixIndent::compute_vector(int n)
// only sum across procs one time // only sum across procs one time
if (indenter_flag == 0) { if (indenter_flag == 0) {
MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(indenter, indenter_all, 4, MPI_DOUBLE, MPI_SUM, world);
indenter_flag = 1; indenter_flag = 1;
} }
return indenter_all[n+1]; return indenter_all[n + 1];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
parse optional parameters at end of input line parse input args for geometry of indenter
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixIndent::options(int narg, char **arg) int FixIndent::geometry(int narg, char **arg)
{ {
if (narg < 0) utils::missing_cmd_args(FLERR, "fix indent", error); if (narg < 0) utils::missing_cmd_args(FLERR, "fix indent", error);
istyle = NONE; istyle = NONE;
xstr = ystr = zstr = rstr = pstr = nullptr; xstr = ystr = zstr = rstr = pstr = nullptr;
xvalue = yvalue = zvalue = rvalue = pvalue = 0.0; xvalue = yvalue = zvalue = rvalue = pvalue = 0.0;
// sphere
if (strcmp(arg[0], "sphere") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error);
if (utils::strmatch(arg[1], "^v_")) {
xstr = utils::strdup(arg[1] + 2);
} else
xvalue = utils::numeric(FLERR, arg[1], false, lmp);
if (utils::strmatch(arg[2], "^v_")) {
ystr = utils::strdup(arg[2] + 2);
} else
yvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
zstr = utils::strdup(arg[3] + 2);
} else
zvalue = utils::numeric(FLERR, arg[3], false, lmp);
if (utils::strmatch(arg[4], "^v_")) {
rstr = utils::strdup(arg[4] + 2);
} else
rvalue = utils::numeric(FLERR, arg[4], false, lmp);
istyle = SPHERE;
return 5;
}
// cylinder
if (strcmp(arg[0], "cylinder") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error);
if (strcmp(arg[1], "x") == 0) {
cdim = 0;
if (utils::strmatch(arg[2], "^v_")) {
ystr = utils::strdup(arg[2] + 2);
} else
yvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
zstr = utils::strdup(arg[3] + 2);
} else
zvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[1], "y") == 0) {
cdim = 1;
if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2] + 2);
} else
xvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
zstr = utils::strdup(arg[3] + 2);
} else
zvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[1], "z") == 0) {
cdim = 2;
if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2] + 2);
} else
xvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
ystr = utils::strdup(arg[3] + 2);
} else
yvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else
error->all(FLERR, "Unknown fix indent cylinder argument: {}", arg[1]);
if (utils::strmatch(arg[4], "^v_")) {
rstr = utils::strdup(arg[4] + 2);
} else
rvalue = utils::numeric(FLERR, arg[4], false, lmp);
istyle = CYLINDER;
return 5;
}
// cone
if (strcmp(arg[0], "cone") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error);
if (strcmp(arg[1], "x") == 0) {
cdim = 0;
if (utils::strmatch(arg[2], "^v_")) {
ystr = utils::strdup(arg[2] + 2);
} else
yvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
zstr = utils::strdup(arg[3] + 2);
} else
zvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[1], "y") == 0) {
cdim = 1;
if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2] + 2);
} else
xvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
zstr = utils::strdup(arg[3] + 2);
} else
zvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[1], "z") == 0) {
cdim = 2;
if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2] + 2);
} else
xvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
ystr = utils::strdup(arg[3] + 2);
} else
yvalue = utils::numeric(FLERR, arg[3], false, lmp);
} else
error->all(FLERR, "Unknown fix indent cone argument: {}", arg[1]);
if (utils::strmatch(arg[4], "^v_")) {
rlostr = utils::strdup(arg[4] + 2);
} else
rlovalue = utils::numeric(FLERR, arg[4], false, lmp);
if (utils::strmatch(arg[5], "^v_")) {
rhistr = utils::strdup(arg[5] + 2);
} else
rhivalue = utils::numeric(FLERR, arg[5], false, lmp);
if (utils::strmatch(arg[6], "^v_")) {
lostr = utils::strdup(arg[6] + 2);
} else
lovalue = utils::numeric(FLERR, arg[6], false, lmp);
if (utils::strmatch(arg[7], "^v_")) {
histr = utils::strdup(arg[7] + 2);
} else
hivalue = utils::numeric(FLERR, arg[7], false, lmp);
istyle = CONE;
return 8;
}
// plane
if (strcmp(arg[0], "plane") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error);
if (strcmp(arg[1], "x") == 0)
cdim = 0;
else if (strcmp(arg[1], "y") == 0)
cdim = 1;
else if (strcmp(arg[1], "z") == 0)
cdim = 2;
else
error->all(FLERR, "Unknown fix indent plane argument: {}", arg[1]);
if (utils::strmatch(arg[2], "^v_")) {
pstr = utils::strdup(arg[2] + 2);
} else
pvalue = utils::numeric(FLERR, arg[2], false, lmp);
if (strcmp(arg[3], "lo") == 0)
planeside = -1;
else if (strcmp(arg[3], "hi") == 0)
planeside = 1;
else
error->all(FLERR, "Unknown fix indent plane argument: {}", arg[3]);
istyle = PLANE;
return 4;
}
// invalid istyle arg
error->all(FLERR, "Unknown fix indent argument: {}", arg[0]);
return 0;
}
/* ----------------------------------------------------------------------
parse optional input args
------------------------------------------------------------------------- */
void FixIndent::options(int narg, char **arg)
{
scaleflag = 1; scaleflag = 1;
side = OUTSIDE; side = OUTSIDE;
int iarg = 0; int iarg = 0;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"sphere") == 0) { if (strcmp(arg[iarg], "units") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error);
if (strcmp(arg[iarg + 1], "box") == 0)
if (utils::strmatch(arg[iarg+1],"^v_")) { scaleflag = 0;
xstr = utils::strdup(arg[iarg+1]+2); else if (strcmp(arg[iarg + 1], "lattice") == 0)
} else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); scaleflag = 1;
if (utils::strmatch(arg[iarg+2],"^v_")) { else
ystr = utils::strdup(arg[iarg+2]+2); error->all(FLERR, "Unknown fix indent units argument: {}", arg[iarg + 1]);
} else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (utils::strmatch(arg[iarg+4],"^v_")) {
rstr = utils::strdup(arg[iarg+4]+2);
} else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
istyle = SPHERE;
iarg += 5;
} else if (strcmp(arg[iarg],"cylinder") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error);
if (strcmp(arg[iarg+1],"x") == 0) {
cdim = 0;
if (utils::strmatch(arg[iarg+2],"^v_")) {
ystr = utils::strdup(arg[iarg+2]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"y") == 0) {
cdim = 1;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"z") == 0) {
cdim = 2;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
ystr = utils::strdup(arg[iarg+3]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[iarg+1]);
if (utils::strmatch(arg[iarg+4],"^v_")) {
rstr = utils::strdup(arg[iarg+4]+2);
} else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
istyle = CYLINDER;
iarg += 5;
} else if (strcmp(arg[iarg],"plane") == 0) {
if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error);
if (strcmp(arg[iarg+1],"x") == 0) cdim = 0;
else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1;
else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+1]);
if (utils::strmatch(arg[iarg+2],"^v_")) {
pstr = utils::strdup(arg[iarg+2]+2);
} else pvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"lo") == 0) planeside = -1;
else if (strcmp(arg[iarg+3],"hi") == 0) planeside = 1;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+3]);
istyle = PLANE;
iarg += 4;
} else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error);
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all(FLERR,"Unknown fix indent units argument: {}", arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"side") == 0) { } else if (strcmp(arg[iarg], "side") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent side", error); if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix indent side", error);
if (strcmp(arg[iarg+1],"in") == 0) side = INSIDE; if (strcmp(arg[iarg + 1], "in") == 0)
else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE; side = INSIDE;
else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]); else if (strcmp(arg[iarg + 1], "out") == 0)
side = OUTSIDE;
else
error->all(FLERR, "Unknown fix indent side argument: {}", arg[iarg + 1]);
iarg += 2; iarg += 2;
} else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]);
} else
error->all(FLERR, "Unknown fix indent argument: {}", arg[iarg]);
} }
} }
/* ----------------------------------------------------------------------
determines if a point is inside (true) or outside (false) of a cone
------------------------------------------------------------------------- */
bool FixIndent::PointInsideCone(int dir, double *center, double lo, double hi, double rlo,
double rhi, double *x)
{
if ((x[dir] > hi) || (x[dir] < lo)) return false;
double del[3] = {x[0] - center[0], x[1] - center[1], x[2] - center[2]};
del[dir] = 0.0;
double dist = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]);
double currentradius = rlo + (x[dir] - lo) * (rhi - rlo) / (hi - lo);
if (dist > currentradius) return false;
return true;
}
/* ----------------------------------------------------------------------
distance between an exterior point and a cone
------------------------------------------------------------------------- */
void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi, double rlo,
double rhi, double &x, double &y, double &z)
{
double xp[3], nearest[3], corner1[3], corner2[3];
double point[3] = {x, y, z};
double del[3] = {x - center[0], y - center[1], z - center[2]};
del[dir] = 0.0;
double r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]);
corner1[0] = center[0] + del[0] * rlo / r;
corner1[1] = center[1] + del[1] * rlo / r;
corner1[2] = center[2] + del[2] * rlo / r;
corner1[dir] = lo;
corner2[0] = center[0] + del[0] * rhi / r;
corner2[1] = center[1] + del[1] * rhi / r;
corner2[2] = center[2] + del[2] * rhi / r;
corner2[dir] = hi;
double corner3[3] = {center[0], center[1], center[2]};
corner3[dir] = lo;
double corner4[3] = {center[0], center[1], center[2]};
corner4[dir] = hi;
// initialize distance to a big number
double distsq = 1.0e20;
// check the first triangle
point_on_line_segment(corner1, corner2, point, xp);
distsq = closest(point, xp, nearest, distsq);
// check the second triangle
point_on_line_segment(corner1, corner3, point, xp);
distsq = closest(point, xp, nearest, distsq);
// check the third triangle
point_on_line_segment(corner2, corner4, point, xp);
distsq = closest(point, xp, nearest, distsq);
x -= nearest[0];
y -= nearest[1];
z -= nearest[2];
return;
}
/* ----------------------------------------------------------------------
distance between an interior point and a cone
------------------------------------------------------------------------- */
void FixIndent::DistanceInteriorPoint(int dir, double *center, double lo, double hi, double rlo,
double rhi, double &x, double &y, double &z)
{
double r, dist_disk, dist_surf;
double surflo[3], surfhi[3], xs[3];
double initial_point[3] = {x, y, z};
double point[3] = {0.0, 0.0, 0.0};
// initial check with the two disks
if ((initial_point[dir] - lo) < (hi - initial_point[dir])) {
dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo);
point[dir] = initial_point[dir] - lo;
} else {
dist_disk = (hi - initial_point[dir]) * (hi - initial_point[dir]);
point[dir] = initial_point[dir] - hi;
}
// check with the points in the conical surface
double del[3] = {x - center[0], y - center[1], z - center[2]};
del[dir] = 0.0;
r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]);
surflo[0] = center[0] + del[0] * rlo / r;
surflo[1] = center[1] + del[1] * rlo / r;
surflo[2] = center[2] + del[2] * rlo / r;
surflo[dir] = lo;
surfhi[0] = center[0] + del[0] * rhi / r;
surfhi[1] = center[1] + del[1] * rhi / r;
surfhi[2] = center[2] + del[2] * rhi / r;
surfhi[dir] = hi;
point_on_line_segment(surflo, surfhi, initial_point, xs);
double dx[3] = {initial_point[0] - xs[0], initial_point[1] - xs[1], initial_point[2] - xs[2]};
dist_surf = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
if (dist_surf < dist_disk) {
x = dx[0];
y = dx[1];
z = dx[2];
} else {
x = point[0];
y = point[1];
z = point[2];
}
return;
}
/* ----------------------------------------------------------------------
helper function extracted from region.cpp
------------------------------------------------------------------------- */
void FixIndent::point_on_line_segment(double *a, double *b, double *c, double *d)
{
double ba[3], ca[3];
MathExtra::sub3(b, a, ba);
MathExtra::sub3(c, a, ca);
double t = MathExtra::dot3(ca, ba) / MathExtra::dot3(ba, ba);
if (t <= 0.0) {
d[0] = a[0];
d[1] = a[1];
d[2] = a[2];
} else if (t >= 1.0) {
d[0] = b[0];
d[1] = b[1];
d[2] = b[2];
} else {
d[0] = a[0] + t * ba[0];
d[1] = a[1] + t * ba[1];
d[2] = a[2] + t * ba[2];
}
}
/* ----------------------------------------------------------------------
helper function extracted from region_cone.cpp
------------------------------------------------------------------------- */
double FixIndent::closest(double *x, double *near, double *nearest, double dsq)
{
double dx = x[0] - near[0];
double dy = x[1] - near[1];
double dz = x[2] - near[2];
double rsq = dx * dx + dy * dy + dz * dz;
if (rsq >= dsq) return dsq;
nearest[0] = near[0];
nearest[1] = near[1];
nearest[2] = near[2];
return rsq;
}

17
src/fix_indent.h
View File

@ -49,7 +49,24 @@ class FixIndent : public Fix {
int cdim, varflag; int cdim, varflag;
int ilevel_respa; int ilevel_respa;
char *rlostr, *rhistr, *lostr, *histr;
int rlovar, rhivar, lovar, hivar;
double rlovalue, rhivalue, lovalue, hivalue;
// methods for argument parsing
int geometry(int, char **);
void options(int, char **); void options(int, char **);
// methods for conical indenter
bool PointInsideCone(int, double *, double, double, double, double, double *);
void DistanceExteriorPoint(int, double *, double, double, double, double,
double &, double &, double &);
void DistanceInteriorPoint(int, double *, double, double, double, double,
double &, double &, double &);
void point_on_line_segment(double *, double *, double *, double *);
double closest(double *, double *, double *, double);
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS

3
src/neighbor.cpp
View File

@ -427,6 +427,9 @@ void Neighbor::init()
} }
} }
} else { } else {
if (!force->pair)
error->all(FLERR, "Cannot use collection/interval command without defining a pairstyle");
if (force->pair->finitecutflag) { if (force->pair->finitecutflag) {
finite_cut_flag = 1; finite_cut_flag = 1;
// If cutoffs depend on finite atom sizes, use radii of intervals to find cutoffs // If cutoffs depend on finite atom sizes, use radii of intervals to find cutoffs

2
src/reader_native.cpp
View File

@ -289,7 +289,7 @@ bigint ReaderNative::read_header(double box[3][3], int &boxinfo, int &triclinic,
labelline = line + strlen("ITEM: ATOMS "); labelline = line + strlen("ITEM: ATOMS ");
} }
Tokenizer tokens(labelline); Tokenizer tokens(std::move(labelline));
std::map<std::string, int> labels; std::map<std::string, int> labels;
nwords = 0; nwords = 0;