diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst
index fd048d56fc5..9c2f8d606d9 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* = *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*
@@ -375,6 +376,24 @@ 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 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, 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 plane area. If you want to perform
+strain with constant volume, the :doc:`fix deform <fix_deform>` command using
+*volume* keyword is more likely to suit your needs.
+
+.. note::
+   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 strains on the other dimensions as well. It is recommended to
+   perform preliminary NPT equilibration if necessary using 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
 chosen will not move during the simulation.  For example, if the lower
@@ -702,10 +721,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, scaleyz =
+scalexz = scalexy = yes if periodic in second dimension and not coupled to
+barostat, otherwise no.
 
 ----------
 
diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst
index 23741165092..96db81727a6 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
 """"""""
@@ -54,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
@@ -79,6 +81,29 @@ 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 <include>` command.
 
+.. 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
+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 <run>` or :doc:`minimize <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<Howto_type_labels>` 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.  Interactions that have been disabled, e.g., via
+:doc:`fix_shake <fix_shake>`, will not be checked.
+
 ----------
 
 Restrictions
@@ -88,7 +113,7 @@ This command must come after the simulation box is defined by a
 :doc:`read_data <read_data>`, :doc:`read_restart <read_restart>`, or
 :doc:`create_box <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 86660fbc0ed..32e7d435a89 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
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) {
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index 6fe82277a05..6588553b962 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -88,6 +88,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;
@@ -120,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
@@ -335,6 +337,17 @@ 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);
+      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 = 1;
+      iarg += 2;
     } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
       if (iarg+4 > narg)
         utils::missing_cmd_args(FLERR, fmt::format("fix {} fixedpoint", style), error);
@@ -461,6 +474,32 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
+  if (isochoric) {
+    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.");
+    } else if (p_flag[0] + p_flag[1] + p_flag[2] == 0) {
+      error->all(FLERR,"Isochoric simulation requires at least 1 diagonal dimension to be barostated.");
+    }
+    if (dimension == 3) {
+      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 * 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.");
+    }
+  }
+
   if ((tstat_flag && t_period <= 0.0) ||
       (p_flag[0] && p_period[0] <= 0.0) ||
       (p_flag[1] && p_period[1] <= 0.0) ||
@@ -490,6 +529,16 @@ 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/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;
+    }
     no_change_box = 1;
     if (allremap == 0) restart_pbc = 1;
 
@@ -1072,12 +1121,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 old_volume, new_volume;
+
+  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
 
@@ -1148,6 +1201,7 @@ void FixNH::remap()
 
   // scale diagonal components
   // scale tilt factors with cell, if set
+  if (isochoric) isofac = vol_start;
 
   if (p_flag[0]) {
     oldlo = domain->boxlo[0];
@@ -1155,6 +1209,7 @@ void FixNH::remap()
     expfac = exp(dto*omega_dot[0]);
     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]) {
@@ -1163,6 +1218,7 @@ void FixNH::remap()
     expfac = exp(dto*omega_dot[1]);
     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;
   }
 
@@ -1172,10 +1228,44 @@ void FixNH::remap()
     expfac = exp(dto*omega_dot[2]);
     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]) {
+      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]) {
+      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) {
@@ -1271,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;
@@ -1328,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;
@@ -1383,6 +1481,13 @@ void FixNH::restart(char *buf)
       h0_inv[4] = list[n++];
       h0_inv[5] = list[n++];
     }
+    flag = static_cast<int> (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 09134a92caf..5592cde66ac 100644
--- a/src/fix_nh.h
+++ b/src/fix_nh.h
@@ -121,6 +121,9 @@ 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 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
 
diff --git a/src/label_map.cpp b/src/label_map.cpp
index 9eceb69394c..dd0bc94e493 100644
--- a/src/label_map.cpp
+++ b/src/label_map.cpp
@@ -41,7 +41,9 @@ static const char cite_type_label_framework[] =
     " pages =   {3282--3297}\n"
     "}\n\n";
 
-static const std::string empty;
+namespace {
+const std::string empty;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -52,6 +54,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();
 }
 
@@ -72,6 +75,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;
@@ -151,6 +156,33 @@ 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 (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') {
+      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
     error->all(FLERR, "Unknown labelmap keyword {}", tlabel);
 
@@ -509,9 +541,11 @@ int LabelMap::infer_dihedraltype(const std::vector<std::string> &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;
@@ -561,7 +595,8 @@ int LabelMap::infer_impropertype(const std::vector<std::string> &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++) {
@@ -601,8 +636,8 @@ int LabelMap::infer_impropertype(const std::vector<std::string> &mytypes, std::a
 
 int LabelMap::parse_typelabel(int ntypes, const std::string &label, std::vector<std::string> &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;
 }
@@ -781,3 +816,172 @@ void LabelMap::write_map(const std::string &filename)
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   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
+  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;
+      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) {
+          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);
+          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);
+        }
+      }
+    }
+    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");
+  }
+
+  // some angles are not symmetric, like class2
+  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++) {
+        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]);
+        int inferred_type = atom->lmap->infer_angletype(type[atom1], type[atom2], type[atom3]);
+        if (inferred_type != atype) {
+          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);
+          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);
+        }
+      }
+    }
+    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");
+  }
+
+  // some dihedrals are not symmetric, like class2
+  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++) {
+        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]);
+        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 = 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);
+          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);
+        }
+      }
+    }
+    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");
+  }
+
+  // some impropers are not symmetric, like class2
+  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++) {
+        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]);
+        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 = 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);
+          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);
+        }
+      }
+    }
+    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;
+}
diff --git a/src/label_map.h b/src/label_map.h
index 2ac5f8e7a9b..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,6 +297,8 @@ Currently used when combining data from multiple sources with
   std::unordered_map<std::string, int> dtypelabel_map;    //!< Dihedral label -> type mapping
   std::unordered_map<std::string, int> itypelabel_map;    //!< Improper label -> type mapping
 
+  int check_which_labels[4];    //!< Indicate check for bonds, angles, dihedrals, and/or impropers
+
   /*! \struct Lmap2Lmap
    *  \brief Mapping structure between two LabelMaps
    *
diff --git a/src/modify.cpp b/src/modify.cpp
index 36e4355977c..830f66c3a1e 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"
@@ -322,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();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/molecule.cpp b/src/molecule.cpp
index ecbdd4fdd4e..6b9b0a76c85 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -4631,6 +4631,7 @@ void Molecule::skip_lines(int n, char *line, const std::string &section)
 
 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]) {
@@ -4641,6 +4642,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);
@@ -4652,8 +4654,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++) {
@@ -4663,6 +4667,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);
@@ -4675,8 +4680,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++) {
@@ -4687,6 +4694,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);
@@ -4700,8 +4708,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++) {
@@ -4712,6 +4722,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);
@@ -4725,6 +4736,7 @@ void Molecule::check_labels()
           }
         }
       }
+      if (perfect_labels) utils::logmesg(lmp, "All impropers in molecule '{}' have self-consistent type labels\n", id);
     }
   }
 }
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..b35e6f2d682
--- /dev/null
+++ b/unittest/force-styles/tests/fix-timestep-npt_isoch.yaml
@@ -0,0 +1,77 @@
+---
+lammps_version: 11 Feb 2026
+date_generated: Thu Mar  5 08:39:25 2026
+epsilon: 5e-12
+skip_tests: kokkos_omp omp
+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.46317957704565
+global_vector: ! |-
+  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.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
+   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.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
+   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
+...