tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'github/develop' into fix-dispersion-d3-issues

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-03-01 00:44:14 -05:00
parent ea48161a51 aaa81b2576
commit 38f8299cfb
4 changed files with 86 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/fix_press_langevin.rst
View File

@ -54,7 +54,7 @@ the Langevin equation such as:
f_P = & \frac{N k_B T_{target}}{V} + \frac{1}{V d}\sum_{i=1}^{N} \vec r_i \cdot \vec f_i - P_{target} \\
Q\ddot{L} + \alpha{}\dot{L} = & f_P + \beta(t)\\
L^{n+1} = & L^{n} + bdt\dot{L}^{n} \frac{bdt^{2}}{2Q} \\
L^{n+1} = & L^{n} + bdt\dot{L}^{n} + \frac{bdt^{2}}{2Q} f^{n}_{P} + \frac{bdt}{2Q} \beta^{n+1} \\
\dot{L}^{n+1} = & \alpha\dot{L}^{n} + \frac{dt}{2Q}\left(a f^{n}_{P} + f^{n+1}_{P}\right) + \frac{b}{Q}\beta^{n+1} \\
a = & \frac{1-\frac{\alpha{}dt}{2Q}}{1+\frac{\alpha{}dt}{2Q}} \\
b = & \frac{1}{1+\frac{\alpha{}dt}{2Q}} \\

25
doc/src/region.rst
View File

@ -40,7 +40,7 @@ Syntax
*plane* args = px py pz nx ny nz
px,py,pz = point on the plane (distance units)
nx,ny,nz = direction normal to plane (distance units)
px,py,pz can be a variable (see below)
px,py,pz,nx,ny,nz can be a variable (see below)
*prism* args = xlo xhi ylo yhi zlo zhi xy xz yz
xlo,xhi,ylo,yhi,zlo,zhi = bounds of untilted prism (distance units)
xy = distance to tilt y in x direction (distance units)
@ -211,9 +211,11 @@ equal-style :doc:`variable <variable>`. Likewise, for style *sphere*
and *ellipsoid* the x-, y-, and z- coordinates of the center of the
sphere/ellipsoid can be specified as an equal-style variable. And for
style *cylinder* the two center positions c1 and c2 for the location of
the cylinder axes can be specified as a equal-style variable. For style
*cone* and *prism* all properties can be defined via equal-style variables. For
style *plane* the point can be defined via equal-style variables.
the cylinder axes can be specified as a equal-style variable. For styles
*block*, *cone*, *prism*, and *plane* all properties can be defined via
equal-style variables. For style *plane*, the components of the direction
vector normal to plane should be either all constants or all defined by
equal-style variables.
If the value is a variable, it should be specified as v_name, where
name is the variable name. In this case, the variable will be
@ -226,6 +228,21 @@ keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent radius or have
a time dependent position of the sphere or cylinder region.
.. note::
Whenever a region property, such as a coordinate or an upper/lower
bound, is defined via an equal-style variable, the variable should
not cause any of the region boundaries to move
too far within a single timestep. Otherwise, bad dynamics will occur.
"Too far" means a small fraction of the approximate distance of
closest approach between two particles, which for the case of Lennard-Jones
particles is the distance of the energy minimum while for granular
particles it is their diameter. An example is a rapidly varying direction
vector in region plane since a small change in the normal to plane will
shift the region surface far away from the region point by a large displacement.
Similarly, bad dynamics can also occur for fast changing variables employed
in the move/rotate options.
See the :doc:`Howto tricilinc <Howto_triclinic>` page for a
geometric description of triclinic boxes, as defined by LAMMPS, and
how to transform these parameters to and from other commonly used

65
src/region_plane.cpp
View File

@ -25,9 +25,11 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr)
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr),
nxstr(nullptr), nystr(nullptr), nzstr(nullptr)
{
xvar = yvar = zvar = 0.0;
nxvar = nyvar = nzvar = 0.0;
options(narg - 8, &arg[8]);
@ -61,15 +63,34 @@ RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
zstyle = CONSTANT;
}
int nxstyle = (utils::strmatch(arg[5], "^v_")) ? 1 : 0;
int nystyle = (utils::strmatch(arg[6], "^v_")) ? 1 : 0;
int nzstyle = (utils::strmatch(arg[7], "^v_")) ? 1 : 0;
if (!nxstyle && !nystyle && !nzstyle) nstyle = CONSTANT;
if (nxstyle && nystyle && nzstyle) nstyle = VARIABLE;
if (nstyle == CONSTANT) {
normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} else if (nstyle == VARIABLE) {
normal[0] = 0.0;
normal[1] = 0.0;
normal[2] = 0.0;
nxstr = utils::strdup(arg[5] + 2);
nystr = utils::strdup(arg[6] + 2);
nzstr = utils::strdup(arg[7] + 2);
varshape = 1;
} else {
error->all(FLERR, "The components of the normal vector should be either all variables or all constants");
}
if (varshape) {
variable_check();
RegPlane::shape_update();
}
normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
// enforce unit normal
double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
@ -95,6 +116,9 @@ RegPlane::~RegPlane()
delete[] xstr;
delete[] ystr;
delete[] zstr;
delete[] nxstr;
delete[] nystr;
delete[] nzstr;
delete[] contact;
}
@ -174,6 +198,20 @@ void RegPlane::shape_update()
if (ystyle == VARIABLE) yp = yscale * input->variable->compute_equal(yvar);
if (zstyle == VARIABLE) zp = zscale * input->variable->compute_equal(zvar);
if (nstyle == VARIABLE) {
normal[0] = xscale * input->variable->compute_equal(nxvar);
normal[1] = yscale * input->variable->compute_equal(nyvar);
normal[2] = zscale * input->variable->compute_equal(nzvar);
// enforce unit normal
double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
if (rsq == 0.0) error->all(FLERR, "Illegal region plane normal vector: {} {} {}", normal[0], normal[1], normal[2]);
normal[0] /= sqrt(rsq);
normal[1] /= sqrt(rsq);
normal[2] /= sqrt(rsq);
}
}
/* ----------------------------------------------------------------------
@ -202,4 +240,21 @@ void RegPlane::variable_check()
if (!input->variable->equalstyle(zvar))
error->all(FLERR, "Variable {} for region plane is invalid style", zstr);
}
if (nstyle == VARIABLE) {
nxvar = input->variable->find(nxstr);
if (nxvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nxstr);
if (!input->variable->equalstyle(nxvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nxstr);
nyvar = input->variable->find(nystr);
if (nyvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nystr);
if (!input->variable->equalstyle(nyvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nystr);
nzvar = input->variable->find(nzstr);
if (nzvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nzstr);
if (!input->variable->equalstyle(nzvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nzstr);
}
}

4
src/region_plane.h
View File

@ -45,6 +45,10 @@ class RegPlane : public Region {
int zstyle, zvar;
char *xstr, *ystr, *zstr;
int nstyle;
int nxvar, nyvar, nzvar;
char *nxstr, *nystr, *nzstr;
void variable_check();
};