From 2fa05928a13e382ee7a23de98be66ec966ed9885 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Tue, 1 Jul 2025 15:30:42 +0200 Subject: [PATCH 01/37] Added some variable to fix_nh --- src/fix_nh.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 3f43c4b4a6f..8c7be2bd2a5 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -483,6 +483,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : if (p_flag[3]) box_change |= BOX_CHANGE_YZ; if (p_flag[4]) box_change |= BOX_CHANGE_XZ; if (p_flag[5]) box_change |= BOX_CHANGE_XY; + if (isochoric) { + box_change |= BOX_CHANGE_X; + box_change |= BOX_CHANGE_Y; + box_change |= BOX_CHANGE_Z; + } no_change_box = 1; if (allremap == 0) restart_pbc = 1; @@ -1071,6 +1076,10 @@ void FixNH::remap() int *mask = atom->mask; int nlocal = atom->nlocal; double *h = domain->h; + double volume; + + if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; + else volume = domain->xprd * domain->yprd; // omega is not used, except for book-keeping From 9231adddfe336729fd7840b7c3977a1df2659950 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 11:38:48 +0200 Subject: [PATCH 02/37] Added the isochoric flag to fix_nh.h --- src/fix_nh.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/fix_nh.h b/src/fix_nh.h index 4741f3bac5f..b1bcd2a8012 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -120,6 +120,7 @@ class FixNH : public Fix { int scalexz; // 1 if xz scaled with lz int scalexy; // 1 if xy scaled with ly int flipflag; // 1 if box flips are invoked as needed + int isochoric // 1 if isochoric NPT simulation int pre_exchange_flag; // set if pre_exchange needed for box flips From 2c6ff1062cd8e2910f18a1bcb99a47a919b7ce27 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 13:41:58 +0200 Subject: [PATCH 03/37] Works better with ; --- src/fix_nh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_nh.h b/src/fix_nh.h index b1bcd2a8012..57f59856a58 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -120,7 +120,7 @@ class FixNH : public Fix { int scalexz; // 1 if xz scaled with lz int scalexy; // 1 if xy scaled with ly int flipflag; // 1 if box flips are invoked as needed - int isochoric // 1 if isochoric NPT simulation + int isochoric; // 1 if isochoric NPT simulation int pre_exchange_flag; // set if pre_exchange needed for box flips From 1ef72cd48a6d76b97aa2157b1ec2cc40b6a9451e Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 13:42:25 +0200 Subject: [PATCH 04/37] Added isochoric option in fix_nh.cpp --- src/fix_nh.cpp | 135 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 128 insertions(+), 7 deletions(-) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 8c7be2bd2a5..ef967b4cfbb 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -87,6 +87,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : dipole_flag = 0; dlm_flag = 0; p_temp_flag = 0; + isochoric = 0; tcomputeflag = 0; pcomputeflag = 0; @@ -334,6 +335,10 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : dlm_flag = 1; } else error->all(FLERR, "Invalid fix {} update argument: {}", style, arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"isochoric") == 0) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix {} isochoric", style), error); + isochoric = utils::logical(FLERR,arg[iarg+1],false,lmp); + iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix {} fixedpoint", style), error); @@ -454,6 +459,23 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : } } + if (isochoric) { + if (dimension == 3 && (p_flag[0]+p_flag[1]+p_flag[2] > 2)) { + error->all(FLERR,"Cannot perform isochoric NPT with all dimensions barostated."); + } else if (dimension == 2 && (p_flag[0] + p_flag[1] > 1)) { + error->all(FLERR,"Cannot perform 2d isochoric NPT with all dimensions barostated."); + } else if (p_flag[0] + p_flag[1] + p_flag[2] == 0) { + error->all(FLERR,"Cannot perform isochoric NPT with no barostated dimension."); + } + if (dimension == 3) { + if (domain->xperiodic * domain->yperiodic * domain->zperiodic == 0) + error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); + } else { + if (domain->xperiodic * domain->yperiodic == 0) + error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); + } + } + if ((tstat_flag && t_period <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) || @@ -483,10 +505,12 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : if (p_flag[3]) box_change |= BOX_CHANGE_YZ; if (p_flag[4]) box_change |= BOX_CHANGE_XZ; if (p_flag[5]) box_change |= BOX_CHANGE_XY; + // In the isochoric case, dimensions can change size while not being + // barostated to maintain volume/surface if (isochoric) { - box_change |= BOX_CHANGE_X; - box_change |= BOX_CHANGE_Y; - box_change |= BOX_CHANGE_Z; + if (not p_flag[0]) box_change |= BOX_CHANGE_X; + if (not p_flag[1]) box_change |= BOX_CHANGE_Y; + if (dimension == 3 && not p_flag[2]) box_change |= BOX_CHANGE_Z; } no_change_box = 1; if (allremap == 0) restart_pbc = 1; @@ -1070,16 +1094,16 @@ void FixNH::remap() { int i; double oldlo,oldhi; - double expfac; + double expfac, isofac; double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; double *h = domain->h; - double volume; + double old_volume, new_volume; - if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; - else volume = domain->xprd * domain->yprd; + if (dimension == 3) old_volume = domain->xprd * domain->yprd * domain->zprd; + else old_volume = domain->xprd * domain->yprd; // omega is not used, except for book-keeping @@ -1150,11 +1174,13 @@ void FixNH::remap() // scale diagonal components // scale tilt factors with cell, if set + isofac = 1.; if (p_flag[0]) { oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; expfac = exp(dto*omega_dot[0]); + isofac /= expfac; domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; } @@ -1163,6 +1189,7 @@ void FixNH::remap() oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; expfac = exp(dto*omega_dot[1]); + isofac /= expfac; domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; if (scalexy) h[5] *= expfac; @@ -1172,12 +1199,106 @@ void FixNH::remap() oldlo = domain->boxlo[2]; oldhi = domain->boxhi[2]; expfac = exp(dto*omega_dot[2]); + isofac /= expfac; domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; if (scalexz) h[4] *= expfac; if (scaleyz) h[3] *= expfac; } + if (isochoric) { + if (dimension == 3) new_volume = domain->xprd * domain->yprd * domain->zprd; + else new_volume = domain->xprd * domain->yprd; + int psum = p_flag[0] + p_flag[1] + p_flag[2]; + if (psum == 1) { + if (p_flag[0]) { + if (dimension == 3) { + isofac = sqrt(isofac); + // Scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + // Scale z + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; + if (scalexz) h[4] *= isofac; + if (scaleyz) h[3] *= isofac; + } else { + // Scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + } + } else if (p_flag[1]) { + if (dimension == 3) { + isofac = sqrt(isofac); + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + // Scale z + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; + if (scalexz) h[4] *= isofac; + if (scaleyz) h[3] *= isofac; + } else { + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + } + } else { // p_flag[2] == 1, should only happen in 3d. + isofac = sqrt(isofac); + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + // Scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + } + } else if (psum == 2) { + if (not p_flag[0]) { + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + } else if (not p_flag[1]) { + // scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + } else { + // scale z + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; + if (scalexz) h[4] *= isofac; + if (scaleyz) h[3] *= isofac; + } + } else { + error->all(FLERR,"Fix {} cannot maintain volume with 3 barostated dimensions", style); + } + } + // off-diagonal components, second half if (pstyle == TRICLINIC) { From 5cb0dc88db0ecb11faf196641a7528df736c6377 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 14:08:17 +0200 Subject: [PATCH 05/37] Added documentation for fix_nh isochoric option --- doc/src/fix_nh.rst | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 0a4076364cb..edf8abb3fef 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -71,6 +71,7 @@ Syntax *scaleyz* value = *yes* or *no* = scale yz with lz *scalexz* value = *yes* or *no* = scale xz with lz *flip* value = *yes* or *no* = allow or disallow box flips when it becomes highly skewed + *isochoric* value = *yes* or *no* = perform NPT simulation at constant volume *fixedpoint* values = x y z x,y,z = perform barostat dilation/contraction around this point (distance units) *update* value = *dipole* or *dipole/dlm* @@ -375,6 +376,19 @@ acquire ghost atoms around a processor's irregular-shaped subdomain. For extreme values of tilt, LAMMPS may also lose atoms and generate an error. +The *isochoric* keyword allows to maintain constant volume during NPT +simulations. This means that the dimensions not coupled to a barostat are used +to compensate the strain induced on barostated dimensions. If one dimension is +coupled to a barostat, then the remaining dimensions are scaled by the same +multiplicative factor so that the volume stays constant. If two dimensions are +coupled to barostats, then the remaining dimension is scaled so that the volume +stays constant. It is not possible to use this keyword if all the dimensions +are coupled to barostats or if none if them is coupled to barostats. In the +case of 2d simulations, only x and y dimensions can be used to maintain a +constant surface. If you want to perform strain with constant volume, the +:doc:`fix deform ` command using *volume* keyword is more likely to +suit your needs. + The *fixedpoint* keyword specifies the fixed point for barostat volume changes. By default, it is the center of the box. Whatever point is chosen will not move during the simulation. For example, if the lower From f0d59fab4d18a3818a27f38e788ceaffe3298ba3 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 14:13:20 +0200 Subject: [PATCH 06/37] Added some details in isochoric keyword documentation --- doc/src/fix_nh.rst | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index edf8abb3fef..02249d18e4d 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -387,7 +387,15 @@ are coupled to barostats or if none if them is coupled to barostats. In the case of 2d simulations, only x and y dimensions can be used to maintain a constant surface. If you want to perform strain with constant volume, the :doc:`fix deform ` command using *volume* keyword is more likely to -suit your needs. +suit your needs. At the moment, it is not possible to use a single barostated +dimension with a single non-barostated dimension to maintain constant volume. + +.. note:: + If large strain are caused by the barostat because the initial configuration + is far from pressure equilibrium or equilibrated too fast, the system will + see large strain on the other dimensions as well. It is recommended to + perform preliminary NPT equilibration if necessary with standard NPT + simulations. The *fixedpoint* keyword specifies the fixed point for barostat volume changes. By default, it is the center of the box. Whatever point is From 01d6cca8b7c8a02df0b910e0d107ec73a0c98922 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 6 Aug 2025 14:49:51 +0200 Subject: [PATCH 07/37] Added isochoric default value to doc. --- doc/src/fix_nh.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 02249d18e4d..90bf77184bb 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -713,10 +713,10 @@ Related commands Default """"""" -The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = 1, -ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none, -flip = yes, scaleyz = scalexz = scalexy = yes if periodic in second -dimension and not coupled to barostat, otherwise no. +The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = 1, ploop = +1, nreset = 0, drag = 0.0, dilate = all, couple = none, flip = yes, isochoric = +no, scaleyz = scalexz = scalexy = yes if periodic in second dimension and not +coupled to barostat, otherwise no. ---------- From cd2d543bfd3dafd8ced94c501e9317dd04635bc0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 21 Feb 2026 09:46:08 -0500 Subject: [PATCH 08/37] replace non-portable text logical operator with traditional one --- src/fix_nh.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index bc5850b520d..424f675e65d 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -475,10 +475,10 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Cannot perform isochoric NPT with no barostated dimension."); } if (dimension == 3) { - if (domain->xperiodic * domain->yperiodic * domain->zperiodic == 0) + if (domain->xperiodic * domain->yperiodic * domain->zperiodic == 0) error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); } else { - if (domain->xperiodic * domain->yperiodic == 0) + if (domain->xperiodic * domain->yperiodic == 0) error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); } } @@ -515,9 +515,9 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : // In the isochoric case, dimensions can change size while not being // barostated to maintain volume/surface if (isochoric) { - if (not p_flag[0]) box_change |= BOX_CHANGE_X; - if (not p_flag[1]) box_change |= BOX_CHANGE_Y; - if (dimension == 3 && not p_flag[2]) box_change |= BOX_CHANGE_Z; + if (!p_flag[0]) box_change |= BOX_CHANGE_X; + if (!p_flag[1]) box_change |= BOX_CHANGE_Y; + if (dimension == 3 && !p_flag[2]) box_change |= BOX_CHANGE_Z; } no_change_box = 1; if (allremap == 0) restart_pbc = 1; @@ -1279,13 +1279,13 @@ void FixNH::remap() if (scalexy) h[5] *= isofac; } } else if (psum == 2) { - if (not p_flag[0]) { + if (!p_flag[0]) { // Scale x oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; - } else if (not p_flag[1]) { + } else if (!p_flag[1]) { // scale y oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; From 11c953205ecfe6c66d2af4e63f12b56c7c498972 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 23 Feb 2026 20:56:05 +0100 Subject: [PATCH 09/37] Added p_isoch array for isochoric option --- src/fix_nh.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/fix_nh.h b/src/fix_nh.h index 108d2533ff0..20e22e8f899 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -121,7 +121,8 @@ class FixNH : public Fix { int scalexz; // 1 if xz scaled with lz int scalexy; // 1 if xy scaled with ly int flipflag; // 1 if box flips are invoked as needed - int isochoric; // 1 if isochoric NPT simulation + int isochoric; // 1 if isochoric NPT simulation + int p_isoch[3]; // 1 if dimension is used for isochoric simulation int pre_exchange_flag; // set if pre_exchange needed for box flips From c4d2516e5c34ae520e8d586c8e6a16dee7d0e4e2 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 23 Feb 2026 20:56:27 +0100 Subject: [PATCH 10/37] Refactored isochoric option to allow for dimension choice --- src/fix_nh.cpp | 140 ++++++++++++++++--------------------------------- 1 file changed, 46 insertions(+), 94 deletions(-) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 424f675e65d..f10d03627d2 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -121,6 +121,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < 6; i++) { p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0; p_flag[i] = 0; + if (i < 3) p_isoch[i] = 0; } // process keywords @@ -338,7 +339,14 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : iarg += 2; } else if (strcmp(arg[iarg],"isochoric") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix {} isochoric", style), error); - isochoric = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg+1], "x") == 0) p_isoch[0] = 1; + else if (strcmp(arg[iarg+1], "y") == 0) p_isoch[1] = 1; + else if (strcmp(arg[iarg+1], "z") == 0) p_isoch[2] = 1; + else if (strcmp(arg[iarg+1], "xy") == 0) p_isoch[0] = p_isoch[1] = 1; + else if (strcmp(arg[iarg+1], "yz") == 0) p_isoch[1] = p_isoch[2] = 1; + else if (strcmp(arg[iarg+1], "xz") == 0) p_isoch[0] = p_isoch[2] = 1; + else error->all(FLERR,"Illegal fix {} isochoric option: {}", style, arg[iarg+1]); + isochoric = true; iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) @@ -467,7 +475,14 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : } if (isochoric) { - if (dimension == 3 && (p_flag[0]+p_flag[1]+p_flag[2] > 2)) { + for (int i; i < 3; i++) { + if (p_flag[i]) { + if (p_isoch[i]) { + error->all(FLERR,"Cannot use barostated dimension as isochoric dimension."); + } + } + } + if (dimension == 3 && (p_flag[0] + p_flag[1] + p_flag[2] > 2)) { error->all(FLERR,"Cannot perform isochoric NPT with all dimensions barostated."); } else if (dimension == 2 && (p_flag[0] + p_flag[1] > 1)) { error->all(FLERR,"Cannot perform 2d isochoric NPT with all dimensions barostated."); @@ -480,6 +495,8 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : } else { if (domain->xperiodic * domain->yperiodic == 0) error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); + if (p_isoch[2]) + error->all(FLERR, "2d Isochoric NPT cannot use z dimension."); } } @@ -515,9 +532,9 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : // In the isochoric case, dimensions can change size while not being // barostated to maintain volume/surface if (isochoric) { - if (!p_flag[0]) box_change |= BOX_CHANGE_X; - if (!p_flag[1]) box_change |= BOX_CHANGE_Y; - if (dimension == 3 && !p_flag[2]) box_change |= BOX_CHANGE_Z; + if (p_isoch[0]) box_change |= BOX_CHANGE_X; + if (p_isoch[1]) box_change |= BOX_CHANGE_Y; + if (dimension == 3 && p_isoch[2]) box_change |= BOX_CHANGE_Z; } no_change_box = 1; if (allremap == 0) restart_pbc = 1; @@ -1214,95 +1231,30 @@ void FixNH::remap() } if (isochoric) { - if (dimension == 3) new_volume = domain->xprd * domain->yprd * domain->zprd; - else new_volume = domain->xprd * domain->yprd; - int psum = p_flag[0] + p_flag[1] + p_flag[2]; - if (psum == 1) { - if (p_flag[0]) { - if (dimension == 3) { - isofac = sqrt(isofac); - // Scale y - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; - domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; - if (scalexy) h[5] *= isofac; - // Scale z - oldlo = domain->boxlo[2]; - oldhi = domain->boxhi[2]; - domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; - domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; - if (scalexz) h[4] *= isofac; - if (scaleyz) h[3] *= isofac; - } else { - // Scale y - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; - domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; - if (scalexy) h[5] *= isofac; - } - } else if (p_flag[1]) { - if (dimension == 3) { - isofac = sqrt(isofac); - // Scale x - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; - domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; - // Scale z - oldlo = domain->boxlo[2]; - oldhi = domain->boxhi[2]; - domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; - domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; - if (scalexz) h[4] *= isofac; - if (scaleyz) h[3] *= isofac; - } else { - // Scale x - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; - domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; - } - } else { // p_flag[2] == 1, should only happen in 3d. - isofac = sqrt(isofac); - // Scale x - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; - domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; - // Scale y - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; - domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; - if (scalexy) h[5] *= isofac; - } - } else if (psum == 2) { - if (!p_flag[0]) { - // Scale x - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; - domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; - } else if (!p_flag[1]) { - // scale y - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; - domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; - if (scalexy) h[5] *= isofac; - } else { - // scale z - oldlo = domain->boxlo[2]; - oldhi = domain->boxhi[2]; - domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; - domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; - if (scalexz) h[4] *= isofac; - if (scaleyz) h[3] *= isofac; - } - } else { - error->all(FLERR,"Fix {} cannot maintain volume with 3 barostated dimensions", style); + int iso_sum = p_isoch[0] + p_isoch[1] + p_isoch[2]; + if (iso_sum == 2) isofac = sqrt(isofac); + if (p_isoch[0]) { + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + } + if (p_isoch[1]) { + // Scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + } + if (p_isoch[2]) { + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; + if (scalexz) h[4] *= isofac; + if (scaleyz) h[3] *= isofac; } } From 7c52c973ca488b6b8e2ece8b5427d59c707d87ba Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 23 Feb 2026 21:03:02 +0100 Subject: [PATCH 11/37] Rephrased the documentation of isochoric option according to changes. --- doc/src/fix_nh.rst | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index a67a440af66..64aec19d25e 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -71,7 +71,7 @@ Syntax *scaleyz* value = *yes* or *no* = scale yz with lz *scalexz* value = *yes* or *no* = scale xz with lz *flip* value = *yes* or *no* = allow or disallow box flips when it becomes highly skewed - *isochoric* value = *yes* or *no* = perform NPT simulation at constant volume + *isochoric* = *x* or *y* or *z* or *xy* or *yz* or *xz* *fixedpoint* values = x y z x,y,z = perform barostat dilation/contraction around this point (distance units) *update* value = *dipole* or *dipole/dlm* @@ -377,24 +377,21 @@ For extreme values of tilt, LAMMPS may also lose atoms and generate an error. The *isochoric* keyword allows to maintain constant volume during NPT -simulations. This means that the dimensions not coupled to a barostat are used -to compensate the strain induced on barostated dimensions. If one dimension is -coupled to a barostat, then the remaining dimensions are scaled by the same -multiplicative factor so that the volume stays constant. If two dimensions are -coupled to barostats, then the remaining dimension is scaled so that the volume -stays constant. It is not possible to use this keyword if all the dimensions -are coupled to barostats or if none if them is coupled to barostats. In the -case of 2d simulations, only x and y dimensions can be used to maintain a -constant surface. If you want to perform strain with constant volume, the -:doc:`fix deform ` command using *volume* keyword is more likely to -suit your needs. At the moment, it is not possible to use a single barostated -dimension with a single non-barostated dimension to maintain constant volume. +simulations. The values following the isochoric dimension indicates the +dimensions to use in that regard: "x" indicates the x dimension, "yz" (no +space) indicates the y and z dimensions. The selected dimensions are scaled to +compensate the strain of the barostat and keep the system at a constant volume +(or surface in 2d). It is not possible to use this keyword if all the +dimensions are coupled to barostats. In the case of 2d simulations, only x and +y dimensions can be used to maintain a constant surface. If you want to perform +strain with constant volume, the :doc:`fix deform ` command using +*volume* keyword is more likely to suit your needs. .. note:: If large strain are caused by the barostat because the initial configuration is far from pressure equilibrium or equilibrated too fast, the system will see large strain on the other dimensions as well. It is recommended to - perform preliminary NPT equilibration if necessary with standard NPT + perform preliminary NPT equilibration if necessary using standard NPT simulations. The *fixedpoint* keyword specifies the fixed point for barostat volume From 00eadbec425619506197cfd64956279f9d915c22 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 23 Feb 2026 21:04:49 +0100 Subject: [PATCH 12/37] Typo --- doc/src/fix_nh.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 64aec19d25e..20189b42502 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -388,9 +388,9 @@ strain with constant volume, the :doc:`fix deform ` command using *volume* keyword is more likely to suit your needs. .. note:: - If large strain are caused by the barostat because the initial configuration + If large strains are caused by the barostat because the initial configuration is far from pressure equilibrium or equilibrated too fast, the system will - see large strain on the other dimensions as well. It is recommended to + see large strains on the other dimensions as well. It is recommended to perform preliminary NPT equilibration if necessary using standard NPT simulations. @@ -722,9 +722,9 @@ Default """"""" The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = 1, ploop = -1, nreset = 0, drag = 0.0, dilate = all, couple = none, flip = yes, isochoric = -no, scaleyz = scalexz = scalexy = yes if periodic in second dimension and not -coupled to barostat, otherwise no. +1, nreset = 0, drag = 0.0, dilate = all, couple = none, flip = yes, scaleyz = +scalexz = scalexy = yes if periodic in second dimension and not coupled to +barostat, otherwise no. ---------- From a977271093cc5a39373bad6e3fd317bf1887a03f Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 23 Feb 2026 21:09:47 +0100 Subject: [PATCH 13/37] Corrected periodic boundary check for isochoric dimensions --- src/fix_nh.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index f10d03627d2..078c4f6ea9d 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -487,14 +487,14 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : } else if (dimension == 2 && (p_flag[0] + p_flag[1] > 1)) { error->all(FLERR,"Cannot perform 2d isochoric NPT with all dimensions barostated."); } else if (p_flag[0] + p_flag[1] + p_flag[2] == 0) { - error->all(FLERR,"Cannot perform isochoric NPT with no barostated dimension."); + error->all(FLERR,"Isochoric simulation requires at least 1 diagonal dimension to be barostated."); } if (dimension == 3) { - if (domain->xperiodic * domain->yperiodic * domain->zperiodic == 0) - error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); + if (domain->xperiodic * p_isoch[0] + domain->yperiodic * p_isoch[1] + domain->zperiodic * p_isoch[2] == 0) + error->all(FLERR, "Isochoric NPT requires periodic boundary conditions along isochoric dimensions."); } else { - if (domain->xperiodic * domain->yperiodic == 0) - error->all(FLERR, "Isochoric NPT requires periodic boundary conditions."); + if (domain->xperiodic * p_isoch[0] + domain->yperiodic * p_isoch[0] == 0) + error->all(FLERR, "Isochoric NPT requires periodic boundary conditions along isochoric dimensions."); if (p_isoch[2]) error->all(FLERR, "2d Isochoric NPT cannot use z dimension."); } From ca277fef824fc3ec688685ec8e2299574db5f70f Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Wed, 25 Feb 2026 11:46:37 +0100 Subject: [PATCH 14/37] Rephrased *isochoric* keyword explanation. Clarified the usage of the *isochoric* keyword in NPT simulations by making the description of the scaling process clearer. --- doc/src/fix_nh.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 20189b42502..9c2f8d606d9 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -376,14 +376,14 @@ acquire ghost atoms around a processor's irregular-shaped subdomain. For extreme values of tilt, LAMMPS may also lose atoms and generate an error. -The *isochoric* keyword allows to maintain constant volume during NPT -simulations. The values following the isochoric dimension indicates the +The *isochoric* keyword allows to maintain constant volume when barostating +up to two dimensions with this fix. The values following the isochoric keyword indicates the dimensions to use in that regard: "x" indicates the x dimension, "yz" (no -space) indicates the y and z dimensions. The selected dimensions are scaled to -compensate the strain of the barostat and keep the system at a constant volume -(or surface in 2d). It is not possible to use this keyword if all the +space) indicates the y and z dimensions, etc. The selected dimensions are scaled to +compensate the strain induced by the barostat and keep the system at a constant volume +(or area in 2d). It is not possible to use this keyword if all the dimensions are coupled to barostats. In the case of 2d simulations, only x and -y dimensions can be used to maintain a constant surface. If you want to perform +y dimensions can be used to maintain a constant plane area. If you want to perform strain with constant volume, the :doc:`fix deform ` command using *volume* keyword is more likely to suit your needs. From 356de88b462cd41686ca2cf9f600e8f72fab6437 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 2 Mar 2026 11:57:55 +0100 Subject: [PATCH 15/37] Changed the algorithm to recompute scaling factor with regard to reference volume. --- src/fix_nh.cpp | 41 ++++++++++++++++++++++++++++++++--------- src/fix_nh.h | 1 + 2 files changed, 33 insertions(+), 9 deletions(-) diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 078c4f6ea9d..65eb025efef 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -346,7 +346,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1], "yz") == 0) p_isoch[1] = p_isoch[2] = 1; else if (strcmp(arg[iarg+1], "xz") == 0) p_isoch[0] = p_isoch[2] = 1; else error->all(FLERR,"Illegal fix {} isochoric option: {}", style, arg[iarg+1]); - isochoric = true; + isochoric = 1; iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) @@ -530,8 +530,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : if (p_flag[4]) box_change |= BOX_CHANGE_XZ; if (p_flag[5]) box_change |= BOX_CHANGE_XY; // In the isochoric case, dimensions can change size while not being - // barostated to maintain volume/surface + // barostated to maintain volume/area if (isochoric) { + if (dimension == 3) vol_start = domain->xprd * domain->yprd * domain->zprd; + else vol_start = domain->xprd * domain->yprd; + if (p_isoch[0]) box_change |= BOX_CHANGE_X; if (p_isoch[1]) box_change |= BOX_CHANGE_Y; if (dimension == 3 && p_isoch[2]) box_change |= BOX_CHANGE_Z; @@ -1198,24 +1201,24 @@ void FixNH::remap() // scale diagonal components // scale tilt factors with cell, if set - isofac = 1.; + if (isochoric) isofac = vol_start; if (p_flag[0]) { oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; expfac = exp(dto*omega_dot[0]); - isofac /= expfac; domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; + if (isochoric) isofac /= domain->boxhi[0] - domain->boxlo[0]; } if (p_flag[1]) { oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; expfac = exp(dto*omega_dot[1]); - isofac /= expfac; domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; + if (isochoric) isofac /= domain->boxhi[1] - domain->boxlo[1]; if (scalexy) h[5] *= expfac; } @@ -1223,25 +1226,30 @@ void FixNH::remap() oldlo = domain->boxlo[2]; oldhi = domain->boxhi[2]; expfac = exp(dto*omega_dot[2]); - isofac /= expfac; domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; + if (isochoric) isofac /= domain->boxhi[2] - domain->boxlo[2]; if (scalexz) h[4] *= expfac; if (scaleyz) h[3] *= expfac; } if (isochoric) { + + // We remove remaining dimensions so that only scale factors are left + // in isofac + for (int i = 0; i < 3; i++) { + if (p_isoch[i] || !p_flag[i]) isofac /= (domain->boxhi[i]-domain->boxlo[i]); + } int iso_sum = p_isoch[0] + p_isoch[1] + p_isoch[2]; if (iso_sum == 2) isofac = sqrt(isofac); + if (p_isoch[0]) { - // Scale x oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; } if (p_isoch[1]) { - // Scale y oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; @@ -1353,8 +1361,9 @@ int FixNH::size_restart_global() int nsize = 2; if (tstat_flag) nsize += 1 + 2*mtchain; if (pstat_flag) { - nsize += 16 + 2*mpchain; + nsize += 17 + 2*mpchain; if (deviatoric_flag) nsize += 6; + if (isochoric) nsize += 4; } return nsize; @@ -1410,6 +1419,13 @@ int FixNH::pack_restart_data(double *list) list[n++] = h0_inv[4]; list[n++] = h0_inv[5]; } + list[n++] = isochoric; + if (isochoric) { + list[n++] = p_isoch[0]; + list[n++] = p_isoch[1]; + list[n++] = p_isoch[2]; + list[n++] = vol_start; + } } return n; @@ -1465,6 +1481,13 @@ void FixNH::restart(char *buf) h0_inv[4] = list[n++]; h0_inv[5] = list[n++]; } + flag = static_cast (list[n++]); + if (flag) { + p_isoch[0] = list[n++]; + p_isoch[1] = list[n++]; + p_isoch[2] = list[n++]; + vol_start = list[n++]; + } } } diff --git a/src/fix_nh.h b/src/fix_nh.h index 20e22e8f899..5592cde66ac 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -123,6 +123,7 @@ class FixNH : public Fix { int flipflag; // 1 if box flips are invoked as needed int isochoric; // 1 if isochoric NPT simulation int p_isoch[3]; // 1 if dimension is used for isochoric simulation + double vol_start; // reference volume for isochoric simulation int pre_exchange_flag; // set if pre_exchange needed for box flips From f06ae1d9613c5488537c298d6faa6327053002a0 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Mon, 2 Mar 2026 11:58:41 +0100 Subject: [PATCH 16/37] =?UTF-8?q?Added=20isochoric=20NPT=20option=20to=20t?= =?UTF-8?q?he=20OPENMP=20nos=C3=A9-hoover=20version.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/OPENMP/fix_nh_omp.cpp | 40 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/OPENMP/fix_nh_omp.cpp b/src/OPENMP/fix_nh_omp.cpp index 64b70f131a8..b33c499e1eb 100644 --- a/src/OPENMP/fix_nh_omp.cpp +++ b/src/OPENMP/fix_nh_omp.cpp @@ -45,6 +45,7 @@ using dbl3_t = struct { double x,y,z; }; void FixNHOMP::remap() { double oldlo,oldhi,expfac; + double isofac; double * const * _noalias const x = atom->x; const int * _noalias const mask = atom->mask; @@ -121,6 +122,8 @@ void FixNHOMP::remap() } } + if (isochoric) isofac = vol_start; + // scale diagonal components // scale tilt factors with cell, if set @@ -128,16 +131,20 @@ void FixNHOMP::remap() oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; expfac = exp(dto*omega_dot[0]); + isofac /= expfac; domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; + if (isochoric) isofac /= domain->boxhi[0] - domain->boxlo[0]; } if (p_flag[1]) { oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; expfac = exp(dto*omega_dot[1]); + isofac /= expfac; domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; + if (isochoric) isofac /= domain->boxhi[1] - domain->boxlo[1]; if (scalexy) h[5] *= expfac; } @@ -145,12 +152,45 @@ void FixNHOMP::remap() oldlo = domain->boxlo[2]; oldhi = domain->boxhi[2]; expfac = exp(dto*omega_dot[2]); + isofac /= expfac; domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; + if (isochoric) isofac /= domain->boxhi[2] - domain->boxlo[2]; if (scalexz) h[4] *= expfac; if (scaleyz) h[3] *= expfac; } + if (isochoric) { + for (int i = 0; i < 3; i++) { + if (p_isoch[i] || !p_flag[i]) isofac /= (domain->boxhi[i]-domain->boxlo[i]); + } + int iso_sum = p_isoch[0] + p_isoch[1] + p_isoch[2]; + if (iso_sum == 2) isofac = sqrt(isofac); + if (p_isoch[0]) { + // Scale x + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + domain->boxlo[0] = (oldlo-fixedpoint[0])*isofac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*isofac + fixedpoint[0]; + } + if (p_isoch[1]) { + // Scale y + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + domain->boxlo[1] = (oldlo-fixedpoint[1])*isofac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*isofac + fixedpoint[1]; + if (scalexy) h[5] *= isofac; + } + if (p_isoch[2]) { + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + domain->boxlo[2] = (oldlo-fixedpoint[2])*isofac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*isofac + fixedpoint[2]; + if (scalexz) h[4] *= isofac; + if (scaleyz) h[3] *= isofac; + } + } + // off-diagonal components, second half if (pstyle == TRICLINIC) { From 877f8989144a73b75d0067904754455f426129a0 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Thu, 5 Mar 2026 13:24:40 +0100 Subject: [PATCH 17/37] Added falsy test file for isochoric npt option --- .../tests/fix-timestep-npt_isoch.yaml | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 unittest/force-styles/tests/fix-timestep-npt_isoch.yaml diff --git a/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml b/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml new file mode 100644 index 00000000000..9bb1145f8e5 --- /dev/null +++ b/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml @@ -0,0 +1,78 @@ +--- +lammps_version: 11 Feb 2026 +tags: generated +date_generated: Tue Feb 24 22:21:53 2026 +epsilon: 2e-12 +skip_tests: +prerequisites: ! | + atom full + fix npt +pre_commands: ! "" +post_commands: ! | + fix test solute npt temp 50.0 ${t_target} 1.0 x 1.0 1.0 100.0 isochoric yz dilate solute +input_file: in.fourmol +natoms: 29 +global_scalar: 328.4631795770457 +global_vector: ! |- + 37 2.0163300805844036 6.868643347719453 40.009180447755085 1.3817019243250745 3.174879058675895 39.53076155633592 0.08220041314309319 0 0 0.06739859756637817 0 0 0.07396099294609239 0.0013777525606211414 -0.0002323660741360924 0.12551892778874418 0.005390245194737114 -0.00023178792751798444 33.65766302261048 1.3649414080498525 7.95065114472879 15.93386223945713 1.0015379719701611 155.26851851065686 3.316011832025461e-17 0 0 101.55405864200037 0 0 0.014697578072112753 0.00027378791194084884 -4.6175941937593954e-05 11.740666503448047 0.02165166767693729 4.00364844049439e-05 -0.045336760079463065 +run_pos: ! |2 + 1 -4.2450026720914469e-01 2.3765488067337897e+00 -1.8553951506250499e-01 + 2 2.0620037399528446e-01 2.8301681103454497e+00 -8.4475981511827758e-01 + 3 -8.8638127401289601e-01 1.1776321251958146e+00 -6.2456828151456190e-01 + 4 -1.8405314926523220e+00 1.4194856960730240e+00 -1.2249586360582549e+00 + 5 -1.1049197679538540e+00 8.8510927694432073e-01 3.6157867187936521e-01 + 6 1.6672379177118302e-01 2.4093319920558898e-01 -1.2213919064025323e+00 + 7 2.4386800671995701e-01 -2.4742649792130145e-02 -2.4276585353553299e+00 + 8 1.1431675163651747e+00 -4.7464766396768709e-01 -6.5305684324807434e-01 + 9 1.3709735202206712e+00 -2.5039312060211660e-01 2.4358209406108866e-01 + 10 2.0836025556741902e+00 -1.3951422822128796e+00 -9.6042562259698983e-01 + 11 1.8171879155801172e+00 -1.9179751428337113e+00 -1.8350762321894063e+00 + 12 3.1359492314593886e+00 -4.7799278708981685e-01 -1.5795327784541051e+00 + 13 4.2677801202291565e+00 -8.6758823062012969e-01 -1.5970995171423015e+00 + 14 2.7020060611324688e+00 -4.0475204816924926e-01 -2.5770873677445190e+00 + 15 3.1003099963414122e+00 5.2852411352995610e-01 -1.2118431383428288e+00 + 16 2.7670057300891884e+00 -2.3168661302253746e+00 -2.5231155287456453e-02 + 17 2.2857702249417500e+00 -2.0181056016090082e+00 1.1069767682551275e+00 + 18 2.1384791188033843e+00 3.0177261773770208e+00 -3.5160827596876225e+00 + 19 1.5349125211132961e+00 2.6315969880333707e+00 -4.2472859440220647e+00 + 20 2.7641167828863153e+00 3.6833419064000221e+00 -3.9380850623312638e+00 + 21 4.9064454390208301e+00 -4.0751205255383196e+00 -3.6215576073601046e+00 + 22 4.3687453488627543e+00 -4.2054270536772504e+00 -4.4651491269372565e+00 + 23 5.7374928154769504e+00 -3.5763355905184966e+00 -3.8820297194230728e+00 + 24 2.0684115301174013e+00 3.1518221747664397e+00 3.1554242678474576e+00 + 25 1.2998381073113014e+00 3.2755513587518097e+00 2.5092990173114837e+00 + 26 2.5807438597688113e+00 4.0120175892854135e+00 3.2133398379059099e+00 + 27 -1.9613581876744359e+00 -4.3556300596085160e+00 2.1101467673534788e+00 + 28 -2.7406520384725965e+00 -4.0207251278130975e+00 1.5828689861678511e+00 + 29 -1.3108232656499081e+00 -3.5992986322410760e+00 2.2680459788743503e+00 +run_vel: ! |2 + 1 3.1079869415450291e-03 7.3871555358947247e-03 2.3942177040894141e-03 + 2 1.2292885324635353e-03 5.9496206849300227e-04 -1.3209048610036884e-03 + 3 -3.3735579283857915e-03 -6.7301140990895696e-03 -2.4073605207107162e-03 + 4 3.5000970748313005e-04 -1.9822758230524052e-03 4.9426282386654690e-04 + 5 -3.1524024773206304e-03 -3.9424753875481463e-03 3.2195242469472132e-05 + 6 -1.2497154284030917e-02 1.5421825842421117e-02 1.9392600476126957e-02 + 7 6.1730125377847155e-04 -4.2752514399470223e-03 -2.1780483103692065e-02 + 8 2.0588959379309571e-03 -8.5640007755960284e-04 9.7792059940798783e-03 + 9 6.9835687868646736e-04 2.5942857996567631e-03 7.3668971579547143e-03 + 10 9.2214978485370240e-03 -8.5665005632045033e-03 -4.7888901436116517e-03 + 11 -1.5672555728599635e-03 -1.5574913339948266e-03 -1.0077848765233858e-03 + 12 7.2781311154948233e-04 -3.1346133568288417e-04 -5.3916246566641638e-04 + 13 -1.9351169409221212e-03 1.4104094126296746e-03 -3.3832590292616298e-04 + 14 1.4628968356840087e-03 -1.1268515315162020e-03 -1.4151104844813652e-03 + 15 1.4371750061283967e-04 -2.4290706916696225e-04 1.8666352958446584e-03 + 16 7.4725903664237536e-03 -5.2483339515333599e-03 -1.8325238298268338e-02 + 17 -5.3984346144911495e-03 3.7214246344641011e-03 1.4622593169796826e-02 + 18 -6.0936815808025862e-04 -9.3774557532468582e-04 -3.3558072507805731e-04 + 19 -6.9919768291957119e-04 -3.6060777270430031e-03 4.2833405289822791e-03 + 20 4.7777805013736515e-03 5.1003745845520452e-03 1.8002873923729241e-03 + 21 -9.5568188553430398e-04 1.6594630943762931e-04 -1.8199788009966615e-04 + 22 -3.3137518957653462e-03 -2.8683968287936054e-03 3.6384389958326871e-03 + 23 2.4209481134686401e-04 -4.5457709985051130e-03 2.7663581642115042e-03 + 24 2.5447450568861086e-04 4.8412447786110117e-04 -4.8021914527341357e-04 + 25 4.3722771097312743e-03 -4.5184411669545515e-03 2.5200952006556795e-03 + 26 -1.9250110555001179e-03 -3.0342169883610837e-03 3.5062814567984532e-03 + 27 -2.6510179146429716e-04 3.6306203629019116e-04 -5.6235585400647747e-04 + 28 -2.3068708109787484e-04 -8.5663070212203200e-04 2.1302563179109169e-03 + 29 -2.5054744388303732e-03 -1.6773997805290820e-04 2.8436699761004796e-03 +... From 7de590e81e751bc73e8e1f9dc498b88118c3ed3f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 5 Mar 2026 08:42:05 -0500 Subject: [PATCH 18/37] flush small values for stress or global scalars or vectors to zero to avoid false positives --- unittest/force-styles/test_fix_timestep.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/unittest/force-styles/test_fix_timestep.cpp b/unittest/force-styles/test_fix_timestep.cpp index eb1f1c80362..83d2acace79 100644 --- a/unittest/force-styles/test_fix_timestep.cpp +++ b/unittest/force-styles/test_fix_timestep.cpp @@ -197,6 +197,9 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // run_stress, if enabled if (ifix->thermo_virial) { auto *stress = ifix->virial; + // avoid false positives on tiny stresses. force to zero instead. + for (int i = 0; i < 6; ++i) + if (fabs(stress[i]) < 1.0e-13) stress[i] = 0.0; block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0], stress[1], stress[2], stress[3], stress[4], stress[5]); writer.emit_block("run_stress", block); @@ -205,6 +208,8 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // global scalar if (ifix->scalar_flag) { double value = ifix->compute_scalar(); + // avoid false positives on tiny values. force to zero instead. + if (fabs(value) < 1.0e-13) value = 0.0; writer.emit("global_scalar", value); } @@ -212,8 +217,13 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) if (ifix->vector_flag) { int num = ifix->size_vector; block = std::to_string(num); - for (int i = 0; i < num; ++i) - block += fmt::format(" {}", ifix->compute_vector(i)); + double value; + for (int i = 0; i < num; ++i) { + // avoid false positives on tiny values. force to zero instead. + value = ifix->compute_vector(i); + if (fabs(value) < 1.0e-13) value = 0.0; + block += fmt::format(" {:23.16e}", value); + } writer.emit_block("global_vector", block); } } From 7520387c2be7bda81bc0a55e7e8d867b6a4971a5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 5 Mar 2026 08:42:51 -0500 Subject: [PATCH 19/37] update test reference --- .../tests/fix-timestep-npt_isoch.yaml | 79 +++++++++---------- 1 file changed, 39 insertions(+), 40 deletions(-) diff --git a/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml b/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml index 9bb1145f8e5..b35e6f2d682 100644 --- a/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml +++ b/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml @@ -1,9 +1,8 @@ --- lammps_version: 11 Feb 2026 -tags: generated -date_generated: Tue Feb 24 22:21:53 2026 -epsilon: 2e-12 -skip_tests: +date_generated: Thu Mar 5 08:39:25 2026 +epsilon: 5e-12 +skip_tests: kokkos_omp omp prerequisites: ! | atom full fix npt @@ -12,27 +11,27 @@ post_commands: ! | fix test solute npt temp 50.0 ${t_target} 1.0 x 1.0 1.0 100.0 isochoric yz dilate solute input_file: in.fourmol natoms: 29 -global_scalar: 328.4631795770457 +global_scalar: 328.46317957704565 global_vector: ! |- - 37 2.0163300805844036 6.868643347719453 40.009180447755085 1.3817019243250745 3.174879058675895 39.53076155633592 0.08220041314309319 0 0 0.06739859756637817 0 0 0.07396099294609239 0.0013777525606211414 -0.0002323660741360924 0.12551892778874418 0.005390245194737114 -0.00023178792751798444 33.65766302261048 1.3649414080498525 7.95065114472879 15.93386223945713 1.0015379719701611 155.26851851065686 3.316011832025461e-17 0 0 101.55405864200037 0 0 0.014697578072112753 0.00027378791194084884 -4.6175941937593954e-05 11.740666503448047 0.02165166767693729 4.00364844049439e-05 -0.045336760079463065 + 37 2.0163300805843996e+00 6.8686433477194528e+00 4.0009180447754908e+01 1.3817019243250837e+00 3.1748790586759470e+00 3.9530761556335840e+01 8.2200413143093287e-02 0.0000000000000000e+00 0.0000000000000000e+00 6.7398597566378296e-02 0.0000000000000000e+00 0.0000000000000000e+00 7.3960992946092569e-02 1.3777525606211501e-03 -2.3236607413609238e-04 1.2551892778874446e-01 5.3902451947371406e-03 -2.3178792751798425e-04 3.3657663022610407e+01 1.3649414080498525e+00 7.9506511447287549e+00 1.5933862239457341e+01 1.0015379719701940e+00 1.5526851851065626e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0155405864200073e+02 0.0000000000000000e+00 0.0000000000000000e+00 1.4697578072112789e-02 2.7378791194085052e-04 -4.6175941937593947e-05 1.1740666503448098e+01 2.1651667676937509e-02 4.0036484404943838e-05 -4.5336760079463030e-02 run_pos: ! |2 - 1 -4.2450026720914469e-01 2.3765488067337897e+00 -1.8553951506250499e-01 - 2 2.0620037399528446e-01 2.8301681103454497e+00 -8.4475981511827758e-01 - 3 -8.8638127401289601e-01 1.1776321251958146e+00 -6.2456828151456190e-01 - 4 -1.8405314926523220e+00 1.4194856960730240e+00 -1.2249586360582549e+00 - 5 -1.1049197679538540e+00 8.8510927694432073e-01 3.6157867187936521e-01 - 6 1.6672379177118302e-01 2.4093319920558898e-01 -1.2213919064025323e+00 - 7 2.4386800671995701e-01 -2.4742649792130145e-02 -2.4276585353553299e+00 - 8 1.1431675163651747e+00 -4.7464766396768709e-01 -6.5305684324807434e-01 - 9 1.3709735202206712e+00 -2.5039312060211660e-01 2.4358209406108866e-01 - 10 2.0836025556741902e+00 -1.3951422822128796e+00 -9.6042562259698983e-01 - 11 1.8171879155801172e+00 -1.9179751428337113e+00 -1.8350762321894063e+00 - 12 3.1359492314593886e+00 -4.7799278708981685e-01 -1.5795327784541051e+00 - 13 4.2677801202291565e+00 -8.6758823062012969e-01 -1.5970995171423015e+00 - 14 2.7020060611324688e+00 -4.0475204816924926e-01 -2.5770873677445190e+00 - 15 3.1003099963414122e+00 5.2852411352995610e-01 -1.2118431383428288e+00 - 16 2.7670057300891884e+00 -2.3168661302253746e+00 -2.5231155287456453e-02 - 17 2.2857702249417500e+00 -2.0181056016090082e+00 1.1069767682551275e+00 + 1 -4.2450026720914380e-01 2.3765488067337799e+00 -1.8553951506250765e-01 + 2 2.0620037399528446e-01 2.8301681103454399e+00 -8.4475981511827669e-01 + 3 -8.8638127401289601e-01 1.1776321251958102e+00 -6.2456828151456101e-01 + 4 -1.8405314926523220e+00 1.4194856960730124e+00 -1.2249586360582558e+00 + 5 -1.1049197679538540e+00 8.8510927694431807e-01 3.6157867187936255e-01 + 6 1.6672379177118302e-01 2.4093319920558276e-01 -1.2213919064025323e+00 + 7 2.4386800671995701e-01 -2.4742649792133697e-02 -2.4276585353553282e+00 + 8 1.1431675163651747e+00 -4.7464766396769242e-01 -6.5305684324807523e-01 + 9 1.3709735202206712e+00 -2.5039312060212371e-01 2.4358209406108955e-01 + 10 2.0836025556741902e+00 -1.3951422822128867e+00 -9.6042562259698983e-01 + 11 1.8171879155801172e+00 -1.9179751428337148e+00 -1.8350762321894072e+00 + 12 3.1359492314593886e+00 -4.7799278708982129e-01 -1.5795327784541033e+00 + 13 4.2677801202291565e+00 -8.6758823062013590e-01 -1.5970995171423006e+00 + 14 2.7020060611324688e+00 -4.0475204816925725e-01 -2.5770873677445199e+00 + 15 3.1003099963414122e+00 5.2852411352994988e-01 -1.2118431383428314e+00 + 16 2.7670057300891884e+00 -2.3168661302253772e+00 -2.5231155287457341e-02 + 17 2.2857702249417500e+00 -2.0181056016090118e+00 1.1069767682551248e+00 18 2.1384791188033843e+00 3.0177261773770208e+00 -3.5160827596876225e+00 19 1.5349125211132961e+00 2.6315969880333707e+00 -4.2472859440220647e+00 20 2.7641167828863153e+00 3.6833419064000221e+00 -3.9380850623312638e+00 @@ -46,23 +45,23 @@ run_pos: ! |2 28 -2.7406520384725965e+00 -4.0207251278130975e+00 1.5828689861678511e+00 29 -1.3108232656499081e+00 -3.5992986322410760e+00 2.2680459788743503e+00 run_vel: ! |2 - 1 3.1079869415450291e-03 7.3871555358947247e-03 2.3942177040894141e-03 - 2 1.2292885324635353e-03 5.9496206849300227e-04 -1.3209048610036884e-03 - 3 -3.3735579283857915e-03 -6.7301140990895696e-03 -2.4073605207107162e-03 - 4 3.5000970748313005e-04 -1.9822758230524052e-03 4.9426282386654690e-04 - 5 -3.1524024773206304e-03 -3.9424753875481463e-03 3.2195242469472132e-05 - 6 -1.2497154284030917e-02 1.5421825842421117e-02 1.9392600476126957e-02 - 7 6.1730125377847155e-04 -4.2752514399470223e-03 -2.1780483103692065e-02 - 8 2.0588959379309571e-03 -8.5640007755960284e-04 9.7792059940798783e-03 - 9 6.9835687868646736e-04 2.5942857996567631e-03 7.3668971579547143e-03 - 10 9.2214978485370240e-03 -8.5665005632045033e-03 -4.7888901436116517e-03 - 11 -1.5672555728599635e-03 -1.5574913339948266e-03 -1.0077848765233858e-03 - 12 7.2781311154948233e-04 -3.1346133568288417e-04 -5.3916246566641638e-04 - 13 -1.9351169409221212e-03 1.4104094126296746e-03 -3.3832590292616298e-04 - 14 1.4628968356840087e-03 -1.1268515315162020e-03 -1.4151104844813652e-03 - 15 1.4371750061283967e-04 -2.4290706916696225e-04 1.8666352958446584e-03 - 16 7.4725903664237536e-03 -5.2483339515333599e-03 -1.8325238298268338e-02 - 17 -5.3984346144911495e-03 3.7214246344641011e-03 1.4622593169796826e-02 + 1 3.1079869415452143e-03 7.3871555358953128e-03 2.3942177040896695e-03 + 2 1.2292885324635978e-03 5.9496206849300531e-04 -1.3209048610038083e-03 + 3 -3.3735579283860023e-03 -6.7301140990903008e-03 -2.4073605207110106e-03 + 4 3.5000970748310560e-04 -1.9822758230523024e-03 4.9426282386654658e-04 + 5 -3.1524024773206404e-03 -3.9424753875481446e-03 3.2195242469553752e-05 + 6 -1.2497154284031036e-02 1.5421825842421211e-02 1.9392600476127026e-02 + 7 6.1730125377847686e-04 -4.2752514399470275e-03 -2.1780483103692155e-02 + 8 2.0588959379311735e-03 -8.5640007755977620e-04 9.7792059940798870e-03 + 9 6.9835687868644817e-04 2.5942857996567484e-03 7.3668971579546735e-03 + 10 9.2214978485368627e-03 -8.5665005632043108e-03 -4.7888901436116188e-03 + 11 -1.5672555728599583e-03 -1.5574913339948741e-03 -1.0077848765234028e-03 + 12 7.2781311154947268e-04 -3.1346133568292423e-04 -5.3916246566644110e-04 + 13 -1.9351169409221175e-03 1.4104094126296922e-03 -3.3832590292615772e-04 + 14 1.4628968356840187e-03 -1.1268515315162148e-03 -1.4151104844813249e-03 + 15 1.4371750061283972e-04 -2.4290706916690652e-04 1.8666352958447035e-03 + 16 7.4725903664238256e-03 -5.2483339515334154e-03 -1.8325238298268442e-02 + 17 -5.3984346144911937e-03 3.7214246344641237e-03 1.4622593169796923e-02 18 -6.0936815808025862e-04 -9.3774557532468582e-04 -3.3558072507805731e-04 19 -6.9919768291957119e-04 -3.6060777270430031e-03 4.2833405289822791e-03 20 4.7777805013736515e-03 5.1003745845520452e-03 1.8002873923729241e-03 From b26822a480fb5515da4af6a5877b7d5a1fde0c93 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 13 Mar 2026 23:05:15 -0400 Subject: [PATCH 20/37] add label_map check for consistent type checks verify self-consistency between type labels and constituent atom types. not actually used anywhere yet --- src/label_map.cpp | 105 ++++++++++++++++++++++++++++++++++++++++++++++ src/label_map.h | 2 + 2 files changed, 107 insertions(+) diff --git a/src/label_map.cpp b/src/label_map.cpp index c05aec0d41e..f0545c98f93 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -781,3 +781,108 @@ void LabelMap::write_map(const std::string &filename) fclose(fp); } } + +/* ---------------------------------------------------------------------- + check type label self-consistency +------------------------------------------------------------------------- */ + +void LabelMap::check_labels() +{ + int *type = atom->type; + tagint *tag = atom->tag; + // in rare cases, bonds are not symmetric. only check if newton on for bonds + if (force->newton_bond && check_which_labels[0]) { + for (int i = 0; i < atom->nlocal; i++) { + int atom1 = i; + for (int j = 0; j < atom->num_bond[i]; j++) { + int btype = atom->bond_type[i][j]; + int atom2 = atom->map(atom->bond_atom[i][j]); + if (atom2<0) printf("atom2 oops %d\n",atom2); + int inferred_type = atom->lmap->infer_bondtype(type[atom1], type[atom2]); + if (inferred_type != btype) { + std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); + std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); + std::string blabel = atom->lmap->find_label(btype, Atom::BOND); + if (inferred_type == -btype) + error->warning(FLERR, "Bond betwen atoms {}, {} has constituent atom types ({}, {}) in reverse order compared " + "to its bond type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); + else error->warning(FLERR, "Bond between atoms {}, {} has constituent atom types ({}, {}) that do not match " + "its type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); + } + } + } + } + // some angles are not symmetric, like class2 + if (check_which_labels[1]) { + for (int i = 0; i < atom->nlocal; i++) { + for (int j = 0; j < atom->num_angle[i]; j++) { + int atype = atom->angle_type[i][j]; + int atom1 = atom->map(atom->angle_atom1[i][j]); + int atom2 = atom->map(atom->angle_atom2[i][j]); + int atom3 = atom->map(atom->angle_atom3[i][j]); + int inferred_type = atom->lmap->infer_angletype(type[atom1], type[atom2], type[atom3]); + if (inferred_type != atype) { + std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); + std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); + std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); + std::string alabel = atom->lmap->find_label(atype, Atom::ANGLE); + if (inferred_type == -atype) + error->warning(FLERR, "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) in reverse order compared " + "to its angle type label ({})", tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, atom3_label, alabel); + else error->warning(FLERR, "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) that do not match its " + "type label ({})", tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, atom3_label, alabel); + } + } + } + } + // some dihedrals are not symmetric, like class2 + if (check_which_labels[2]) { + for (int i = 0; i < atom->nlocal; i++) { + for (int j = 0; j < atom->num_dihedral[i]; j++) { + int dtype = atom->dihedral_type[i][j]; + int atom1 = atom->map(atom->dihedral_atom1[i][j]); + int atom2 = atom->map(atom->dihedral_atom2[i][j]); + int atom3 = atom->map(atom->dihedral_atom3[i][j]); + int atom4 = atom->map(atom->dihedral_atom4[i][j]); + int inferred_type = atom->lmap->infer_dihedraltype(type[atom1], type[atom2], type[atom3], type[atom4]); + if (inferred_type != dtype) { + std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); + std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); + std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); + std::string atom4_label = atom->lmap->find_label(type[atom4], Atom::ATOM); + std::string dlabel = atom->lmap->find_label(dtype, Atom::DIHEDRAL); + if (inferred_type == -dtype) + error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) in reverse order compared to its " + "dihedral type label ({})", atom1, atom2, atom3, atom4, atom1_label, atom2_label, atom3_label, atom4_label, dlabel); + else error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) that do not match its " + "dihedral label ({})", atom1, atom2, atom3, atom4, atom1_label, atom2_label, atom3_label, atom4_label, dlabel); + } + } + } + } + // some impropers are not symmetric, like class2 + if (check_which_labels[3]) { + for (int i = 0; i < atom->nlocal; i++) { + for (int j = 0; j < atom->num_improper[i]; j++) { + int itype = atom->improper_type[i][j]; + int atom1 = atom->map(atom->improper_atom1[i][j]); + int atom2 = atom->map(atom->improper_atom2[i][j]); + int atom3 = atom->map(atom->improper_atom3[i][j]); + int atom4 = atom->map(atom->improper_atom4[i][j]); + int inferred_type = atom->lmap->infer_impropertype(type[atom1], type[atom2], type[atom3], type[atom4]); + if (inferred_type != itype) { + std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); + std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); + std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); + std::string atom4_label = atom->lmap->find_label(type[atom4], Atom::ATOM); + std::string ilabel = atom->lmap->find_label(itype, Atom::IMPROPER); + if (inferred_type == -itype) + error->warning(FLERR, "Improper containing atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) in a different order compared to its " + "improper type label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, ilabel); + else error->warning(FLERR, "Improper containing atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) that do not match its " + "improper label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, ilabel); + } + } + } + } +} diff --git a/src/label_map.h b/src/label_map.h index 2ac5f8e7a9b..287bac25a62 100644 --- a/src/label_map.h +++ b/src/label_map.h @@ -355,6 +355,8 @@ Currently used when combining data from multiple sources with * * \param filename file name */ void write_map(const std::string &filename); + + void check_labels(); //!< Check if type labels are self-consistent }; } // namespace LAMMPS_NS From 1899d3fca4717cc01c1f2441226a2b3ee2712785 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 00:04:41 -0400 Subject: [PATCH 21/37] add rest of label_map check --- src/label_map.cpp | 13 +++++++++++++ src/label_map.h | 3 +++ 2 files changed, 16 insertions(+) diff --git a/src/label_map.cpp b/src/label_map.cpp index f0545c98f93..25e84605e7e 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -71,6 +71,8 @@ LabelMap::~LabelMap() void LabelMap::reset_type_labels() { + for (int i = 0; i < 4; i++) check_which_labels[i] = 0; + typelabel_map.clear(); typelabel.resize(natomtypes); delete[] lmap2lmap.atom; @@ -150,6 +152,17 @@ void LabelMap::modify_lmap(int narg, char **arg) if (narg != 2) error->all(FLERR, "Incorrect number of arguments for labelmap write command"); write_map(arg[1]); return; + } else if (tlabel == "check_labels") { + if (strchr(arg[1], 'b')) + check_which_labels[0] = 1; + if (strchr(arg[1], 'a')) + check_which_labels[1] = 1; + if (strchr(arg[1], 'd')) + check_which_labels[2] = 1; + if (strchr(arg[1], 'i')) + check_which_labels[3] = 1; + check_labels(); + return; } else error->all(FLERR, "Unknown labelmap keyword {}", tlabel); diff --git a/src/label_map.h b/src/label_map.h index 287bac25a62..f16bedb6c67 100644 --- a/src/label_map.h +++ b/src/label_map.h @@ -293,6 +293,9 @@ Currently used when combining data from multiple sources with std::unordered_map dtypelabel_map; //!< Dihedral label -> type mapping std::unordered_map itypelabel_map; //!< Improper label -> type mapping + int checkflag; //!< Flag to check for self-consistent type labels + int check_which_labels[4]; //!< Indicate check for bonds, angles, dihedrals, and/or impropers + /*! \struct Lmap2Lmap * \brief Mapping structure between two LabelMaps * From 1df0be5e59c06dbf30435444ab2498a89a1bbc89 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 00:27:44 -0400 Subject: [PATCH 22/37] Add runtime check for type label consistency... ...directly to Modify --- src/label_map.h | 7 ++++--- src/modify.cpp | 5 +++++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/label_map.h b/src/label_map.h index f16bedb6c67..854402d9a84 100644 --- a/src/label_map.h +++ b/src/label_map.h @@ -62,6 +62,8 @@ class LabelMap : protected Pointers { int nimpropertypes); ~LabelMap() override; + int checkflag; //!< Flag to check for self-consistent type labels + /*! Process labelmap command from input script * \verbatim embed:rst @@ -279,6 +281,8 @@ Currently used when combining data from multiple sources with * \param fp File pointer for writing */ void write_restart(FILE *); + void check_labels(); //!< Check if type labels are self-consistent + /*! @} */ protected: @@ -293,7 +297,6 @@ Currently used when combining data from multiple sources with std::unordered_map dtypelabel_map; //!< Dihedral label -> type mapping std::unordered_map itypelabel_map; //!< Improper label -> type mapping - int checkflag; //!< Flag to check for self-consistent type labels int check_which_labels[4]; //!< Indicate check for bonds, angles, dihedrals, and/or impropers /*! \struct Lmap2Lmap @@ -358,8 +361,6 @@ Currently used when combining data from multiple sources with * * \param filename file name */ void write_map(const std::string &filename); - - void check_labels(); //!< Check if type labels are self-consistent }; } // namespace LAMMPS_NS diff --git a/src/modify.cpp b/src/modify.cpp index 36e4355977c..c4e3d6abad0 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -25,6 +25,7 @@ #include "group.h" #include "input.h" #include "memory.h" +#include "label_map.h" #include "region.h" #include "update.h" #include "variable.h" @@ -299,6 +300,10 @@ void Modify::init() if (comm->me == 0 && checkall) error->warning(FLERR, "One or more atoms are time integrated more than once" + utils::errorurl(32)); + + // runtime check for type label self-consistency + + if (atom->labelmapflag && atom->lmap->checkflag) atom->lmap->check_labels(); } /* ---------------------------------------------------------------------- From ccda46b70723b7d86e3bba859dd35cde2219e77c Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 15:44:25 -0400 Subject: [PATCH 23/37] move check to Modify::setup, must be after neighboring --- src/label_map.cpp | 28 +++++++++++++++++++++++++++- src/modify.cpp | 8 ++++---- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 25e84605e7e..9003a75ba7a 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -120,6 +120,8 @@ void LabelMap::modify_lmap(int narg, char **arg) if (lmp->citeme) lmp->citeme->add(cite_type_label_framework); + checkflag = 0; + int ntypes; std::vector *labels; std::unordered_map *labels_map; @@ -161,7 +163,7 @@ void LabelMap::modify_lmap(int narg, char **arg) check_which_labels[2] = 1; if (strchr(arg[1], 'i')) check_which_labels[3] = 1; - check_labels(); + checkflag = 1; return; } else error->all(FLERR, "Unknown labelmap keyword {}", tlabel); @@ -804,6 +806,8 @@ void LabelMap::check_labels() int *type = atom->type; tagint *tag = atom->tag; // in rare cases, bonds are not symmetric. only check if newton on for bonds + bool globally_perfect_labels; + bool perfect_labels = true; if (force->newton_bond && check_which_labels[0]) { for (int i = 0; i < atom->nlocal; i++) { int atom1 = i; @@ -813,6 +817,7 @@ void LabelMap::check_labels() if (atom2<0) printf("atom2 oops %d\n",atom2); int inferred_type = atom->lmap->infer_bondtype(type[atom1], type[atom2]); if (inferred_type != btype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string blabel = atom->lmap->find_label(btype, Atom::BOND); @@ -824,8 +829,13 @@ void LabelMap::check_labels() } } } + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels) + utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } + MPI_Barrier(MPI_COMM_WORLD); // some angles are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[1]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_angle[i]; j++) { @@ -835,6 +845,7 @@ void LabelMap::check_labels() int atom3 = atom->map(atom->angle_atom3[i][j]); int inferred_type = atom->lmap->infer_angletype(type[atom1], type[atom2], type[atom3]); if (inferred_type != atype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -847,8 +858,13 @@ void LabelMap::check_labels() } } } + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels) + utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } + MPI_Barrier(MPI_COMM_WORLD); // some dihedrals are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[2]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_dihedral[i]; j++) { @@ -859,6 +875,7 @@ void LabelMap::check_labels() int atom4 = atom->map(atom->dihedral_atom4[i][j]); int inferred_type = atom->lmap->infer_dihedraltype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != dtype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -872,8 +889,13 @@ void LabelMap::check_labels() } } } + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels) + utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } + MPI_Barrier(MPI_COMM_WORLD); // some impropers are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[3]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_improper[i]; j++) { @@ -884,6 +906,7 @@ void LabelMap::check_labels() int atom4 = atom->map(atom->improper_atom4[i][j]); int inferred_type = atom->lmap->infer_impropertype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != itype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -897,5 +920,8 @@ void LabelMap::check_labels() } } } + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels) + utils::logmesg(lmp, "All impropers in the simulation have self-consistent type labels\n"); } } diff --git a/src/modify.cpp b/src/modify.cpp index c4e3d6abad0..830f66c3a1e 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -300,10 +300,6 @@ void Modify::init() if (comm->me == 0 && checkall) error->warning(FLERR, "One or more atoms are time integrated more than once" + utils::errorurl(32)); - - // runtime check for type label self-consistency - - if (atom->labelmapflag && atom->lmap->checkflag) atom->lmap->check_labels(); } /* ---------------------------------------------------------------------- @@ -327,6 +323,10 @@ void Modify::setup(int vflag) for (int i = 0; i < nfix; i++) fix[i]->setup(vflag); else if (update->whichflag == 2) for (int i = 0; i < nfix; i++) fix[i]->min_setup(vflag); + + // runtime check for type label self-consistency + + if (atom->labelmapflag && atom->lmap->checkflag) atom->lmap->check_labels(); } /* ---------------------------------------------------------------------- From 6e62141b21275ca381c51028d84fcdffe22e02bf Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 18:19:48 -0400 Subject: [PATCH 24/37] copy/paste error --- src/label_map.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 9003a75ba7a..df280b59b60 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -883,9 +883,9 @@ void LabelMap::check_labels() std::string dlabel = atom->lmap->find_label(dtype, Atom::DIHEDRAL); if (inferred_type == -dtype) error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) in reverse order compared to its " - "dihedral type label ({})", atom1, atom2, atom3, atom4, atom1_label, atom2_label, atom3_label, atom4_label, dlabel); + "dihedral type label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, dlabel); else error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) that do not match its " - "dihedral label ({})", atom1, atom2, atom3, atom4, atom1_label, atom2_label, atom3_label, atom4_label, dlabel); + "dihedral label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, dlabel); } } } From f1a74fc441d68f732aef51e92e219a8bad6e7ef6 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 18:36:32 -0400 Subject: [PATCH 25/37] add verbose success to molecule 'check_types' --- src/molecule.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/molecule.cpp b/src/molecule.cpp index 5206f0ee7ce..8f09dacd90d 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -4243,6 +4243,7 @@ void Molecule::skip_lines(int n, char *line, const std::string §ion) void Molecule::check_labels() { + bool perfect_labels = true; if (atom->labelmapflag) { // in rare cases, bonds are not symmetric. only check if newton on for bonds if (force->newton_bond && check_which_labels[0]) { @@ -4253,6 +4254,7 @@ void Molecule::check_labels() int atom2 = bond_atom[i][j]; int inferred_type = atom->lmap->infer_bondtype(type[atom1-1], type[atom2-1]); if (inferred_type != btype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1-1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2-1], Atom::ATOM); std::string blabel = atom->lmap->find_label(btype, Atom::BOND); @@ -4264,8 +4266,10 @@ void Molecule::check_labels() } } } + if (perfect_labels) utils::logmesg(lmp, "All bonds in molecule '{}' have self-consistent type labels\n", id); } // some angles are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[1]) { for (int i = 0; i < natoms; i++) { for (int j = 0; j < num_angle[i]; j++) { @@ -4275,6 +4279,7 @@ void Molecule::check_labels() int atom3 = angle_atom3[i][j]; int inferred_type = atom->lmap->infer_angletype(type[atom1-1], type[atom2-1], type[atom3-1]); if (inferred_type != atype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1-1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2-1], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3-1], Atom::ATOM); @@ -4287,8 +4292,10 @@ void Molecule::check_labels() } } } + if (perfect_labels) utils::logmesg(lmp, "All angles in molecule '{}' have self-consistent type labels\n", id); } // some dihedrals are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[2]) { for (int i = 0; i < natoms; i++) { for (int j = 0; j < num_dihedral[i]; j++) { @@ -4299,6 +4306,7 @@ void Molecule::check_labels() int atom4 = dihedral_atom4[i][j]; int inferred_type = atom->lmap->infer_dihedraltype(type[atom1-1], type[atom2-1], type[atom3-1], type[atom4-1]); if (inferred_type != dtype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1-1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2-1], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3-1], Atom::ATOM); @@ -4312,8 +4320,10 @@ void Molecule::check_labels() } } } + if (perfect_labels) utils::logmesg(lmp, "All dihedrals in molecule '{}' have self-consistent type labels\n", id); } // some impropers are not symmetric, like class2 + perfect_labels = true; if (check_which_labels[3]) { for (int i = 0; i < natoms; i++) { for (int j = 0; j < num_improper[i]; j++) { @@ -4324,6 +4334,7 @@ void Molecule::check_labels() int atom4 = improper_atom4[i][j]; int inferred_type = atom->lmap->infer_impropertype(type[atom1-1], type[atom2-1], type[atom3-1], type[atom4-1]); if (inferred_type != itype) { + perfect_labels = false; std::string atom1_label = atom->lmap->find_label(type[atom1-1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2-1], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3-1], Atom::ATOM); @@ -4337,6 +4348,7 @@ void Molecule::check_labels() } } } + if (perfect_labels) utils::logmesg(lmp, "All impropers in molecule '{}' have self-consistent type labels\n", id); } } } From a44868d483327542a32195424adcf1962e01c2b3 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 14 Mar 2026 18:48:24 -0400 Subject: [PATCH 26/37] Update labelmap.rst --- doc/src/labelmap.rst | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index 23741165092..452fc9f71ee 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -10,7 +10,7 @@ Syntax labelmap option args -* *option* = *atom* or *bond* or *angle* or *dihedral* or *improper* or *clear* or *write* +* *option* = *atom* or *bond* or *angle* or *dihedral* or *improper* or *clear* or *write* or *check_labels* .. parsed-literal:: @@ -18,6 +18,8 @@ Syntax *write* arg = filename *atom* or *bond* or *angle* or *dihedral* or *improper* args = list of one or more numeric-type/type-label pairs + *check_labels* value = string + string = string containing any of the following characters: 'b', 'a', 'd', or 'i' Examples """""""" @@ -79,6 +81,24 @@ label mappings to a file as a sequence of *labelmap* commands, so the file can be copied into a new LAMMPS input file or read in using the :doc:`include ` command. +The *check_labels* keyword provides a warning if the type label of a +bond, angle, dihedral, or improper defined in the simulation is not +consistent with the atom types of its constituent atoms. This check is +performed at runtime, i.e., when the :doc:`run ` command is called. +The *check_labels* value is a single string that should contain one or +more of the characters 'b', 'a', 'd', and 'i', which correspond to +bonds, angles, dihedrals, and impropers, respectively. For example, the +keyword/value pair 'check_labels badi' will check all the type labels of +all higher-order interactions, while 'check_labels adi' will only check +type labels for angles, dihedrals, and impropers. The *check_labels* +keyword requires a specific :doc:`type label` format +to infer the types of higher-order interactions. Bond, angle, dihedral, +and improper type labels must contain their constituent atom types +delimited by hyphens, e.g., 'c2-c2-c2-n' for a dihedral that contains +three atoms of type 'c2' and one atom of 'n'. If the constituent atoms +do not have these atom types in the proper order, a warning will be +generated when using this *check_types* keyword. + ---------- Restrictions From 7cf3b20d3fb34ca4fcc726940cd72d1f2fade1a5 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 15 Mar 2026 00:01:17 -0400 Subject: [PATCH 27/37] parse keyword/arg with syntax checks --- src/label_map.cpp | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index df280b59b60..693fd78b373 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -155,14 +155,28 @@ void LabelMap::modify_lmap(int narg, char **arg) write_map(arg[1]); return; } else if (tlabel == "check_labels") { - if (strchr(arg[1], 'b')) - check_which_labels[0] = 1; - if (strchr(arg[1], 'a')) - check_which_labels[1] = 1; - if (strchr(arg[1], 'd')) - check_which_labels[2] = 1; - if (strchr(arg[1], 'i')) - check_which_labels[3] = 1; + if (narg != 2) error->all(FLERR, "Incorrect number of arguments for labelmap write command"); + int i = 0; + char option; + while ((option = arg[1][i++]) != '\0') { + switch (option) { + case 'b': + check_which_labels[0] = 1; + break; + case 'a': + check_which_labels[1] = 1; + break; + case 'd': + check_which_labels[2] = 1; + break; + case 'i': + check_which_labels[3] = 1; + break; + default: + error->all(FLERR, "Labelmap command: Illegal check_labels option {}", option); + break; + } + } checkflag = 1; return; } else From fe65ba6e256718a78b394a858b4e2ec803026f5d Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 15 Mar 2026 00:14:19 -0400 Subject: [PATCH 28/37] bool->int due to old MPI version issues --- src/label_map.cpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 693fd78b373..2621fe0f272 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -820,8 +820,8 @@ void LabelMap::check_labels() int *type = atom->type; tagint *tag = atom->tag; // in rare cases, bonds are not symmetric. only check if newton on for bonds - bool globally_perfect_labels; - bool perfect_labels = true; + int globally_perfect_labels; + int perfect_labels = 1; if (force->newton_bond && check_which_labels[0]) { for (int i = 0; i < atom->nlocal; i++) { int atom1 = i; @@ -831,7 +831,7 @@ void LabelMap::check_labels() if (atom2<0) printf("atom2 oops %d\n",atom2); int inferred_type = atom->lmap->infer_bondtype(type[atom1], type[atom2]); if (inferred_type != btype) { - perfect_labels = false; + perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string blabel = atom->lmap->find_label(btype, Atom::BOND); @@ -843,13 +843,13 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); if (comm->me == 0 && globally_perfect_labels) utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); // some angles are not symmetric, like class2 - perfect_labels = true; + perfect_labels = 1; if (check_which_labels[1]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_angle[i]; j++) { @@ -859,7 +859,7 @@ void LabelMap::check_labels() int atom3 = atom->map(atom->angle_atom3[i][j]); int inferred_type = atom->lmap->infer_angletype(type[atom1], type[atom2], type[atom3]); if (inferred_type != atype) { - perfect_labels = false; + perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -872,13 +872,13 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); if (comm->me == 0 && globally_perfect_labels) utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); // some dihedrals are not symmetric, like class2 - perfect_labels = true; + perfect_labels = 1; if (check_which_labels[2]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_dihedral[i]; j++) { @@ -889,7 +889,7 @@ void LabelMap::check_labels() int atom4 = atom->map(atom->dihedral_atom4[i][j]); int inferred_type = atom->lmap->infer_dihedraltype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != dtype) { - perfect_labels = false; + perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -903,13 +903,13 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); if (comm->me == 0 && globally_perfect_labels) utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); // some impropers are not symmetric, like class2 - perfect_labels = true; + perfect_labels = 1; if (check_which_labels[3]) { for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_improper[i]; j++) { @@ -920,7 +920,7 @@ void LabelMap::check_labels() int atom4 = atom->map(atom->improper_atom4[i][j]); int inferred_type = atom->lmap->infer_impropertype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != itype) { - perfect_labels = false; + perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); @@ -934,7 +934,7 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_C_BOOL, MPI_LAND, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); if (comm->me == 0 && globally_perfect_labels) utils::logmesg(lmp, "All impropers in the simulation have self-consistent type labels\n"); } From bf1ccfda16e8544d56bb2f9aa88decfd4ce715b7 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 15 Mar 2026 12:14:27 -0400 Subject: [PATCH 29/37] more old MPI issues --- src/label_map.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 2621fe0f272..b92216434c2 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -843,8 +843,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == 1) utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -872,8 +872,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == 1) utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -903,8 +903,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == 1) utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -934,8 +934,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == 1) utils::logmesg(lmp, "All impropers in the simulation have self-consistent type labels\n"); } } From 87bdf45b6547820e948c37a056b0ca6c246ebb13 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 15 Mar 2026 12:30:49 -0400 Subject: [PATCH 30/37] ok, this is a STUBS issue too lazy to update the STUBS library --- src/label_map.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index b92216434c2..3a2690f051f 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -843,8 +843,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels == 1) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -872,8 +872,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels == 1) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -903,8 +903,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels == 1) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } MPI_Barrier(MPI_COMM_WORLD); @@ -934,8 +934,8 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_PROD, 0, MPI_COMM_WORLD); - if (comm->me == 0 && globally_perfect_labels == 1) + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All impropers in the simulation have self-consistent type labels\n"); } } From 6609f240722159754ddfadcf86a579246aa8ed96 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 17 Mar 2026 16:05:03 -0400 Subject: [PATCH 31/37] Fix typo in labelmap.rst regarding check_labels --- doc/src/labelmap.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index 452fc9f71ee..334c6403bb0 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -97,7 +97,7 @@ and improper type labels must contain their constituent atom types delimited by hyphens, e.g., 'c2-c2-c2-n' for a dihedral that contains three atoms of type 'c2' and one atom of 'n'. If the constituent atoms do not have these atom types in the proper order, a warning will be -generated when using this *check_types* keyword. +generated when using this *check_labels* keyword. ---------- From 6cff77b18256658cc689ddd212f4abd7c11084a8 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 18 Mar 2026 12:36:34 -0400 Subject: [PATCH 32/37] thanks, Copilot! --- doc/src/labelmap.rst | 34 +++++++++++++++++++--------------- src/label_map.cpp | 24 ++++++++++++------------ 2 files changed, 31 insertions(+), 27 deletions(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index 452fc9f71ee..10a60fb7e7e 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -81,23 +81,27 @@ label mappings to a file as a sequence of *labelmap* commands, so the file can be copied into a new LAMMPS input file or read in using the :doc:`include ` command. +.. versionchanged:: TBD + The *check_labels* keyword provides a warning if the type label of a bond, angle, dihedral, or improper defined in the simulation is not -consistent with the atom types of its constituent atoms. This check is -performed at runtime, i.e., when the :doc:`run ` command is called. -The *check_labels* value is a single string that should contain one or -more of the characters 'b', 'a', 'd', and 'i', which correspond to -bonds, angles, dihedrals, and impropers, respectively. For example, the -keyword/value pair 'check_labels badi' will check all the type labels of -all higher-order interactions, while 'check_labels adi' will only check -type labels for angles, dihedrals, and impropers. The *check_labels* -keyword requires a specific :doc:`type label` format -to infer the types of higher-order interactions. Bond, angle, dihedral, -and improper type labels must contain their constituent atom types -delimited by hyphens, e.g., 'c2-c2-c2-n' for a dihedral that contains -three atoms of type 'c2' and one atom of 'n'. If the constituent atoms -do not have these atom types in the proper order, a warning will be -generated when using this *check_types* keyword. +consistent with the atom types of its constituent atoms. This +consistency check is performed only once, when the simulation is +initialized by the first :doc:`run ` or :doc:`minimize ` +command after this *labelmap* command. The *check_labels* value is a +single string that should contain one or more of the characters 'b', +'a', 'd', and 'i', which correspond to bonds, angles, dihedrals, and +impropers, respectively. For example, the keyword/value pair +'check_labels badi' will check all the type labels of all higher-order +interactions, while 'check_labels adi' will only check type labels for +angles, dihedrals, and impropers. The *check_labels* keyword requires a +specific :doc:`type label` format to infer the types +of higher-order interactions. Bond, angle, dihedral, and improper type +labels must contain their constituent atom types delimited by hyphens, +e.g., 'c2-c2-c2-n' for a dihedral that contains three atoms of type 'c2' +and one atom of 'n'. If the constituent atoms do not have these atom +types in the proper order, a warning will be generated when using this +*check_types* keyword. ---------- diff --git a/src/label_map.cpp b/src/label_map.cpp index 3a2690f051f..9b78c142af2 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -51,6 +51,7 @@ LabelMap::LabelMap(LAMMPS *_lmp, int _natomtypes, int _nbondtypes, int _nanglety { lmap2lmap.atom = lmap2lmap.bond = lmap2lmap.angle = lmap2lmap.dihedral = lmap2lmap.improper = nullptr; + checkflag = 0; reset_type_labels(); } @@ -120,8 +121,6 @@ void LabelMap::modify_lmap(int narg, char **arg) if (lmp->citeme) lmp->citeme->add(cite_type_label_framework); - checkflag = 0; - int ntypes; std::vector *labels; std::unordered_map *labels_map; @@ -155,7 +154,8 @@ void LabelMap::modify_lmap(int narg, char **arg) write_map(arg[1]); return; } else if (tlabel == "check_labels") { - if (narg != 2) error->all(FLERR, "Incorrect number of arguments for labelmap write command"); + if (narg != 2) error->all(FLERR, "Incorrect number of arguments for labelmap check_labels command"); + for (int j = 0; j < 4; j++) check_which_labels[j] = 0; int i = 0; char option; while ((option = arg[1][i++]) != '\0') { @@ -828,7 +828,6 @@ void LabelMap::check_labels() for (int j = 0; j < atom->num_bond[i]; j++) { int btype = atom->bond_type[i][j]; int atom2 = atom->map(atom->bond_atom[i][j]); - if (atom2<0) printf("atom2 oops %d\n",atom2); int inferred_type = atom->lmap->infer_bondtype(type[atom1], type[atom2]); if (inferred_type != btype) { perfect_labels = 0; @@ -836,18 +835,18 @@ void LabelMap::check_labels() std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string blabel = atom->lmap->find_label(btype, Atom::BOND); if (inferred_type == -btype) - error->warning(FLERR, "Bond betwen atoms {}, {} has constituent atom types ({}, {}) in reverse order compared " + error->warning(FLERR, "Bond between atoms {}, {} has constituent atom types ({}, {}) in reverse order compared " "to its bond type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); else error->warning(FLERR, "Bond between atoms {}, {} has constituent atom types ({}, {}) that do not match " "its type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, world); if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(world); // some angles are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[1]) { @@ -872,11 +871,11 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, world); if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(world); // some dihedrals are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[2]) { @@ -903,11 +902,11 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, world); if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(world); // some impropers are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[3]) { @@ -934,8 +933,9 @@ void LabelMap::check_labels() } } } - MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&perfect_labels, &globally_perfect_labels, 1, MPI_INT, MPI_SUM, 0, world); if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All impropers in the simulation have self-consistent type labels\n"); } + checkflag = 0; } From cd432b9f4793241f65edbefd7c7eb9109a2dd5eb Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 18 Mar 2026 15:56:14 -0400 Subject: [PATCH 33/37] skip interactions that have been turned off thanks, Axel! --- doc/src/labelmap.rst | 5 ++++- src/label_map.cpp | 4 ++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index d40c8623cbe..03b781c3c87 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -101,7 +101,10 @@ labels must contain their constituent atom types delimited by hyphens, e.g., 'c2-c2-c2-n' for a dihedral that contains three atoms of type 'c2' and one atom of 'n'. If the constituent atoms do not have these atom types in the proper order, a warning will be generated when using this -*check_labels* keyword. +*check_types* keyword. Interactions that have been disabled, e.g., via +:doc:`fix_shake `, will not be checked. + +---------- Restrictions """""""""""" diff --git a/src/label_map.cpp b/src/label_map.cpp index 9b78c142af2..99a6e47ec4b 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -827,6 +827,7 @@ void LabelMap::check_labels() int atom1 = i; for (int j = 0; j < atom->num_bond[i]; j++) { int btype = atom->bond_type[i][j]; + if (btype < 1) continue; int atom2 = atom->map(atom->bond_atom[i][j]); int inferred_type = atom->lmap->infer_bondtype(type[atom1], type[atom2]); if (inferred_type != btype) { @@ -853,6 +854,7 @@ void LabelMap::check_labels() for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_angle[i]; j++) { int atype = atom->angle_type[i][j]; + if (atype < 1) continue; int atom1 = atom->map(atom->angle_atom1[i][j]); int atom2 = atom->map(atom->angle_atom2[i][j]); int atom3 = atom->map(atom->angle_atom3[i][j]); @@ -882,6 +884,7 @@ void LabelMap::check_labels() for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_dihedral[i]; j++) { int dtype = atom->dihedral_type[i][j]; + if (dtype < 1) continue; int atom1 = atom->map(atom->dihedral_atom1[i][j]); int atom2 = atom->map(atom->dihedral_atom2[i][j]); int atom3 = atom->map(atom->dihedral_atom3[i][j]); @@ -913,6 +916,7 @@ void LabelMap::check_labels() for (int i = 0; i < atom->nlocal; i++) { for (int j = 0; j < atom->num_improper[i]; j++) { int itype = atom->improper_type[i][j]; + if (itype < 1) continue; int atom1 = atom->map(atom->improper_atom1[i][j]); int atom2 = atom->map(atom->improper_atom2[i][j]); int atom3 = atom->map(atom->improper_atom3[i][j]); From c12e30561f32d4c697e884386afc9a795154d661 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 27 Mar 2026 16:33:04 -0400 Subject: [PATCH 34/37] apply clang-format --- src/label_map.cpp | 92 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 67 insertions(+), 25 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 99a6e47ec4b..846daf03c65 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -40,7 +40,9 @@ static const char cite_type_label_framework[] = " pages = {3282--3297}\n" "}\n\n"; -static const std::string empty; +namespace { +const std::string empty; +} /* ---------------------------------------------------------------------- */ @@ -154,7 +156,8 @@ void LabelMap::modify_lmap(int narg, char **arg) write_map(arg[1]); return; } else if (tlabel == "check_labels") { - if (narg != 2) error->all(FLERR, "Incorrect number of arguments for labelmap check_labels command"); + if (narg != 2) + error->all(FLERR, "Incorrect number of arguments for labelmap check_labels command"); for (int j = 0; j < 4; j++) check_which_labels[j] = 0; int i = 0; char option; @@ -537,9 +540,11 @@ int LabelMap::infer_dihedraltype(const std::vector &mytypes) status = parse_typelabel(4, dtypelabel[i], dtypes); if (status != -1) { if (mytypes[0] == dtypes[0] && mytypes[1] == dtypes[1] && mytypes[2] == dtypes[2] && - mytypes[3] == dtypes[3]) return i + 1; + mytypes[3] == dtypes[3]) + return i + 1; if (mytypes[3] == dtypes[0] && mytypes[2] == dtypes[1] && mytypes[1] == dtypes[2] && - mytypes[0] == dtypes[3]) out = -(i + 1); + mytypes[0] == dtypes[3]) + out = -(i + 1); } } return out; @@ -589,7 +594,8 @@ int LabelMap::infer_impropertype(const std::vector &mytypes, std::a status = parse_typelabel(4, itypelabel[i], itypes); if (status != -1) { if (mytypes[0] == itypes[0] && mytypes[1] == itypes[1] && mytypes[2] == itypes[2] && - mytypes[3] == itypes[3]) return i + 1; + mytypes[3] == itypes[3]) + return i + 1; navail_types = 4; avail_types = mytypes; for (int j = 0; j < 4; j++) { @@ -629,8 +635,8 @@ int LabelMap::infer_impropertype(const std::vector &mytypes, std::a int LabelMap::parse_typelabel(int ntypes, const std::string &label, std::vector &types) { - auto out = Tokenizer(label,"-").as_vector(); - if ((int)out.size() != ntypes) return -1; + auto out = Tokenizer(label, "-").as_vector(); + if ((int) out.size() != ntypes) return -1; types = std::move(out); return 1; } @@ -836,10 +842,17 @@ void LabelMap::check_labels() std::string atom2_label = atom->lmap->find_label(type[atom2], Atom::ATOM); std::string blabel = atom->lmap->find_label(btype, Atom::BOND); if (inferred_type == -btype) - error->warning(FLERR, "Bond between atoms {}, {} has constituent atom types ({}, {}) in reverse order compared " - "to its bond type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); - else error->warning(FLERR, "Bond between atoms {}, {} has constituent atom types ({}, {}) that do not match " - "its type label ({})", tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); + error->warning(FLERR, + "Bond between atoms {}, {} has constituent atom types ({}, {}) in " + "reverse order compared " + "to its bond type label ({})", + tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); + else + error->warning( + FLERR, + "Bond between atoms {}, {} has constituent atom types ({}, {}) that do not match " + "its type label ({})", + tag[atom1], tag[atom2], atom1_label, atom2_label, blabel); } } } @@ -866,10 +879,19 @@ void LabelMap::check_labels() std::string atom3_label = atom->lmap->find_label(type[atom3], Atom::ATOM); std::string alabel = atom->lmap->find_label(atype, Atom::ANGLE); if (inferred_type == -atype) - error->warning(FLERR, "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) in reverse order compared " - "to its angle type label ({})", tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, atom3_label, alabel); - else error->warning(FLERR, "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) that do not match its " - "type label ({})", tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, atom3_label, alabel); + error->warning(FLERR, + "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) " + "in reverse order compared " + "to its angle type label ({})", + tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, + atom3_label, alabel); + else + error->warning(FLERR, + "Angle between atoms {}, {}, {} has constituent atom types ({}, {}, {}) " + "that do not match its " + "type label ({})", + tag[atom1], tag[atom2], tag[atom3], atom1_label, atom2_label, + atom3_label, alabel); } } } @@ -889,7 +911,8 @@ void LabelMap::check_labels() int atom2 = atom->map(atom->dihedral_atom2[i][j]); int atom3 = atom->map(atom->dihedral_atom3[i][j]); int atom4 = atom->map(atom->dihedral_atom4[i][j]); - int inferred_type = atom->lmap->infer_dihedraltype(type[atom1], type[atom2], type[atom3], type[atom4]); + int inferred_type = + atom->lmap->infer_dihedraltype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != dtype) { perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); @@ -898,10 +921,19 @@ void LabelMap::check_labels() std::string atom4_label = atom->lmap->find_label(type[atom4], Atom::ATOM); std::string dlabel = atom->lmap->find_label(dtype, Atom::DIHEDRAL); if (inferred_type == -dtype) - error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) in reverse order compared to its " - "dihedral type label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, dlabel); - else error->warning(FLERR, "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) that do not match its " - "dihedral label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, dlabel); + error->warning(FLERR, + "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, " + "{}, {}, {}) in reverse order compared to its " + "dihedral type label ({})", + tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, + atom3_label, atom4_label, dlabel); + else + error->warning(FLERR, + "Dihedral between atoms {}, {}, {}, {} has constituent atom types ({}, " + "{}, {}, {}) that do not match its " + "dihedral label ({})", + tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, + atom3_label, atom4_label, dlabel); } } } @@ -921,7 +953,8 @@ void LabelMap::check_labels() int atom2 = atom->map(atom->improper_atom2[i][j]); int atom3 = atom->map(atom->improper_atom3[i][j]); int atom4 = atom->map(atom->improper_atom4[i][j]); - int inferred_type = atom->lmap->infer_impropertype(type[atom1], type[atom2], type[atom3], type[atom4]); + int inferred_type = + atom->lmap->infer_impropertype(type[atom1], type[atom2], type[atom3], type[atom4]); if (inferred_type != itype) { perfect_labels = 0; std::string atom1_label = atom->lmap->find_label(type[atom1], Atom::ATOM); @@ -930,10 +963,19 @@ void LabelMap::check_labels() std::string atom4_label = atom->lmap->find_label(type[atom4], Atom::ATOM); std::string ilabel = atom->lmap->find_label(itype, Atom::IMPROPER); if (inferred_type == -itype) - error->warning(FLERR, "Improper containing atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) in a different order compared to its " - "improper type label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, ilabel); - else error->warning(FLERR, "Improper containing atoms {}, {}, {}, {} has constituent atom types ({}, {}, {}, {}) that do not match its " - "improper label ({})", tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, atom3_label, atom4_label, ilabel); + error->warning(FLERR, + "Improper containing atoms {}, {}, {}, {} has constituent atom types " + "({}, {}, {}, {}) in a different order compared to its " + "improper type label ({})", + tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, + atom3_label, atom4_label, ilabel); + else + error->warning(FLERR, + "Improper containing atoms {}, {}, {}, {} has constituent atom types " + "({}, {}, {}, {}) that do not match its " + "improper label ({})", + tag[atom1], tag[atom2], tag[atom3], tag[atom4], atom1_label, atom2_label, + atom3_label, atom4_label, ilabel); } } } From bc504e1051f9be808a8225e34eb4c8120aa24ce4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 27 Mar 2026 16:33:15 -0400 Subject: [PATCH 35/37] remove redundant barriers --- src/label_map.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/label_map.cpp b/src/label_map.cpp index 846daf03c65..498880159bb 100644 --- a/src/label_map.cpp +++ b/src/label_map.cpp @@ -860,7 +860,7 @@ void LabelMap::check_labels() if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All bonds in the simulation have self-consistent type labels\n"); } - MPI_Barrier(world); + // some angles are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[1]) { @@ -899,7 +899,7 @@ void LabelMap::check_labels() if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All angles in the simulation have self-consistent type labels\n"); } - MPI_Barrier(world); + // some dihedrals are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[2]) { @@ -941,7 +941,7 @@ void LabelMap::check_labels() if (comm->me == 0 && globally_perfect_labels == comm->nprocs) utils::logmesg(lmp, "All dihedrals in the simulation have self-consistent type labels\n"); } - MPI_Barrier(world); + // some impropers are not symmetric, like class2 perfect_labels = 1; if (check_which_labels[3]) { From 10d7a2ae89da563557c4334bf05bcbb5521ddf49 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 27 Mar 2026 16:36:04 -0400 Subject: [PATCH 36/37] we use versionadded for added features --- doc/src/labelmap.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index 03b781c3c87..cd08036a98b 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -81,7 +81,7 @@ label mappings to a file as a sequence of *labelmap* commands, so the file can be copied into a new LAMMPS input file or read in using the :doc:`include ` command. -.. versionchanged:: TBD +.. versionadded:: TBD The *check_labels* keyword provides a warning if the type label of a bond, angle, dihedral, or improper defined in the simulation is not From d43206291eb72e221c183ac435fc2df3d1a954b2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 27 Mar 2026 16:40:14 -0400 Subject: [PATCH 37/37] small doc and spelling tweaks --- doc/src/labelmap.rst | 14 +++++++------- doc/utils/sphinx-config/false_positives.txt | 2 ++ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst index cd08036a98b..96db81727a6 100644 --- a/doc/src/labelmap.rst +++ b/doc/src/labelmap.rst @@ -56,12 +56,12 @@ or '[' and ']', parenthesis '(' and ')', dash '-', underscore '_', plus '+' and equals '=' signs and more. They must not contain blanks or any other whitespace. Note that type labels must be put in single or double quotation marks if they contain the '#' character or if they contain a -double (") or single quotation mark ('). If the label contains both -a single and a double quotation mark, then triple quotation (""") must -be used. When enclosing a type label with quotation marks, the -LAMMPS input parser may require adding leading or trailing blanks -around the type label so it can identify the enclosing quotation -marks. Those blanks will be removed when defining the label. +double (") or single quotation mark ('). If the label contains both a +single and a double quotation mark, then triple quotation (""") must be +used. When enclosing a type label with quotation marks, the LAMMPS +input parser may require adding leading or trailing blanks around the +type label so it can identify the enclosing quotation marks. Those +blanks will be removed when defining the label. A *labelmap* command can only modify the label map for one type-kind (atom types, bond types, etc). Any number of numeric-type/type-label @@ -113,7 +113,7 @@ This command must come after the simulation box is defined by a :doc:`read_data `, :doc:`read_restart `, or :doc:`create_box ` command. -Label maps are currently not supported when using the KOKKOS package. +Label maps are considered experimental when using the KOKKOS package. Related commands """""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 41238197591..d66b3ba4c3d 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -37,6 +37,7 @@ Addington addstep addtorque adf +adi Adhikari adiam adim @@ -238,6 +239,7 @@ Babaei Baconnier backcolor Baczewski +badi Baer Bagchi Bagi