From 0b1ef95562968061af8b01180e5be9c170d375b6 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 31 Jan 2025 17:11:22 +0200 Subject: [PATCH 01/12] Update methods in region_plane.cpp --- src/region_plane.cpp | 65 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 60 insertions(+), 5 deletions(-) diff --git a/src/region_plane.cpp b/src/region_plane.cpp index 11d96b8266..d65fe2cabe 100644 --- a/src/region_plane.cpp +++ b/src/region_plane.cpp @@ -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); + } } From 1057882126327460f95e29447b11adfd94ce211f Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 31 Jan 2025 17:14:49 +0200 Subject: [PATCH 02/12] Include variables in region_plane.h --- src/region_plane.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/region_plane.h b/src/region_plane.h index 790564e176..be7e68f56a 100644 --- a/src/region_plane.h +++ b/src/region_plane.h @@ -43,6 +43,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(); }; From 72eb284f76ecb1150887c3e31d5910cef1d27911 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 31 Jan 2025 17:22:15 +0200 Subject: [PATCH 03/12] remove whitespace from region_plane.cpp --- src/region_plane.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/region_plane.cpp b/src/region_plane.cpp index d65fe2cabe..597f586d8a 100644 --- a/src/region_plane.cpp +++ b/src/region_plane.cpp @@ -69,7 +69,7 @@ RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : 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); From 418b205362ee3e0aae1941abf53773403a7a3ea6 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 1 Feb 2025 09:30:20 +0200 Subject: [PATCH 04/12] Update region.rst --- doc/src/region.rst | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index 94feee6ad4..2e5bca4c7e 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -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) @@ -210,9 +210,11 @@ equal-style :doc:`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* 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 +*cone* and *plane* all properties can be defined via equal-style variables. +For style *plane*, there is the restriction that all three components +of the direction normal to the plane should be either constant or +defined via 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 From b2b9f2c3e95ae1ee1f57c4da64e56991bdf2351d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 1 Feb 2025 09:31:45 +0200 Subject: [PATCH 05/12] remove whitespace in region.rst --- doc/src/region.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index 2e5bca4c7e..7b95831249 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -213,7 +213,7 @@ 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 styles *cone* and *plane* all properties can be defined via equal-style variables. For style *plane*, there is the restriction that all three components -of the direction normal to the plane should be either constant or +of the direction normal to the plane should be either constant or defined via equal-style variables. If the value is a variable, it should be specified as v_name, where From 150cd216babdc85a656f5279b1b02967dcdb49c6 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 1 Feb 2025 09:46:06 +0200 Subject: [PATCH 06/12] another one whitespace --- doc/src/region.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index 7b95831249..4a6f5ae3b8 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -214,7 +214,7 @@ the cylinder axes can be specified as a equal-style variable. For styles *cone* and *plane* all properties can be defined via equal-style variables. For style *plane*, there is the restriction that all three components of the direction normal to the plane should be either constant or -defined via equal-style variables. +defined via 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 From 262ff223c6f2fe4f4ef8a32bc170987017584b92 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Tue, 4 Feb 2025 13:56:22 +0200 Subject: [PATCH 07/12] make clear that the vector should be all constant or all equal variables --- doc/src/region.rst | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index 4a6f5ae3b8..6a7e22fee7 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -212,9 +212,8 @@ 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 styles *cone* and *plane* all properties can be defined via equal-style variables. -For style *plane*, there is the restriction that all three components -of the direction normal to the plane should be either constant or -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 From 3ed6716b65cbdc3ddc60c1370b1210a163d21205 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 5 Feb 2025 10:35:54 +0200 Subject: [PATCH 08/12] clarify that bad dynamics may occur in the documentation --- doc/src/region.rst | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/doc/src/region.rst b/doc/src/region.rst index 6a7e22fee7..b42fe3d8da 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -226,6 +226,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, it must be ensured + that the variable does not cause any of the region boundaries to move + too far within a single timestep. Otherwise, bad dynamics will occur. + By too far it is meant 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 ` page for a geometric description of triclinic boxes, as defined by LAMMPS, and how to transform these parameters to and from other commonly used From 424c694b6b9d67d9180e30168ce15d4fb02d0203 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 8 Feb 2025 20:56:12 +0200 Subject: [PATCH 09/12] update region.rst --- doc/src/region.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index b42fe3d8da..81107161c4 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -229,8 +229,8 @@ 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, it must be ensured - that the variable does not cause any of the region boundaries to move + 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. By too far it is meant a small fraction of the approximate distance of closest approach between two particles, which for the case of Lennard-Jones From ea973e1d6c90c6f55b0fedd46d99ca8ca4237899 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 8 Feb 2025 20:57:45 +0200 Subject: [PATCH 10/12] Update region.rst --- doc/src/region.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index 81107161c4..6f576ff9b9 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -232,7 +232,7 @@ a time dependent position of the sphere or cylinder region. 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. - By too far it is meant a small fraction of the approximate distance of + "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 From f5099b7c16f813c073eb8dddb885429ce50e0799 Mon Sep 17 00:00:00 2001 From: lasergyro <31989178+lasergyro@users.noreply.github.com> Date: Tue, 18 Feb 2025 17:15:01 +0100 Subject: [PATCH 11/12] Update fix_press_langevin.rst to include full equation 13 Completest equation 13 to match the [paper](https://pubs.aip.org/jcp/article/141/19/194108/152571/Constant-pressure-and-temperature-discrete-time) and the [implementation](https://github.com/lammps/lammps/blob/36e739b734d2c0b12f6ad56da2b3a59a2b35d7b5/src/fix_press_langevin.cpp#L485-L487). --- doc/src/fix_press_langevin.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_press_langevin.rst b/doc/src/fix_press_langevin.rst index 8438d72192..02437bd731 100644 --- a/doc/src/fix_press_langevin.rst +++ b/doc/src/fix_press_langevin.rst @@ -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}} \\ From 9bf19159a4f2b50dfda4531dc3cc24fc24fcf59f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 26 Feb 2025 10:42:52 -0500 Subject: [PATCH 12/12] add block to list of regions with all-variable shape options --- doc/src/region.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/region.rst b/doc/src/region.rst index fbc250ca83..dbe14360ff 100644 --- a/doc/src/region.rst +++ b/doc/src/region.rst @@ -212,7 +212,7 @@ 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 styles -*cone*, *prism*, and *plane* all properties can be defined via +*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.