diff --git a/doc/src/Modify_pair.rst b/doc/src/Modify_pair.rst
index 7263b8fd48..6913204504 100644
--- a/doc/src/Modify_pair.rst
+++ b/doc/src/Modify_pair.rst
@@ -12,24 +12,24 @@ includes some optional methods to enable its use with rRESPA.
 
 Here is a brief description of the class methods in pair.h:
 
-+---------------------------------+-------------------------------------------------------------------+
-| compute                         | workhorse routine that computes pairwise interactions             |
-+---------------------------------+-------------------------------------------------------------------+
-| settings                        | reads the input script line with arguments you define             |
-+---------------------------------+-------------------------------------------------------------------+
-| coeff                           | set coefficients for one i,j type pair                            |
-+---------------------------------+-------------------------------------------------------------------+
-| init_one                        | perform initialization for one i,j type pair                      |
-+---------------------------------+-------------------------------------------------------------------+
-| init_style                      | initialization specific to this pair style                        |
-+---------------------------------+-------------------------------------------------------------------+
-| write & read_restart            | write/read i,j pair coeffs to restart files                       |
-+---------------------------------+-------------------------------------------------------------------+
-| write & read_restart_settings   | write/read global settings to restart files                       |
-+---------------------------------+-------------------------------------------------------------------+
-| single                          | force and energy of a single pairwise interaction between 2 atoms |
-+---------------------------------+-------------------------------------------------------------------+
-| compute_inner/middle/outer      | versions of compute used by rRESPA                                |
-+---------------------------------+-------------------------------------------------------------------+
++---------------------------------+---------------------------------------------------------------------+
+| compute                         | workhorse routine that computes pairwise interactions               |
++---------------------------------+---------------------------------------------------------------------+
+| settings                        | reads the input script line with arguments you define               |
++---------------------------------+---------------------------------------------------------------------+
+| coeff                           | set coefficients for one i,j type pair                              |
++---------------------------------+---------------------------------------------------------------------+
+| init_one                        | perform initialization for one i,j type pair                        |
++---------------------------------+---------------------------------------------------------------------+
+| init_style                      | initialization specific to this pair style                          |
++---------------------------------+---------------------------------------------------------------------+
+| write & read_restart            | write/read i,j pair coeffs to restart files                         |
++---------------------------------+---------------------------------------------------------------------+
+| write & read_restart_settings   | write/read global settings to restart files                         |
++---------------------------------+---------------------------------------------------------------------+
+| single                          | force/r and energy of a single pairwise interaction between 2 atoms |
++---------------------------------+---------------------------------------------------------------------+
+| compute_inner/middle/outer      | versions of compute used by rRESPA                                  |
++---------------------------------+---------------------------------------------------------------------+
 
 The inner/middle/outer routines are optional.
diff --git a/doc/src/pair_nm.rst b/doc/src/pair_nm.rst
index 09bb571da5..2256c7f220 100644
--- a/doc/src/pair_nm.rst
+++ b/doc/src/pair_nm.rst
@@ -195,5 +195,4 @@ none
 
 .. _Dietz:
 
-**(Dietz)** J.D. Dietz, R.S. Hoy, "Facile equilibration of well-entangled
-semi-flexible bead-spring polymer melts" arXiv:2109.11001
+**(Dietz)** Dietz and Hoy, J. Chem Phys, 156, 014103 (2022).
diff --git a/doc/src/pair_python.rst b/doc/src/pair_python.rst
index 3d087565be..35e07dbd11 100644
--- a/doc/src/pair_python.rst
+++ b/doc/src/pair_python.rst
@@ -126,11 +126,11 @@ and *compute_energy*, which both take 3 numerical arguments:
 * itype = the (numerical) type of the first atom
 * jtype = the (numerical) type of the second atom
 
-This functions need to compute the force and the energy, respectively,
-and use the result as return value. The functions need to use the
-*pmap* dictionary to convert the LAMMPS atom type number to the symbolic
-value of the internal potential parameter data structure. Following
-the *LJCutMelt* example, here are the two functions:
+This functions need to compute the (scaled) force and the energy,
+respectively, and use the result as return value. The functions need
+to use the *pmap* dictionary to convert the LAMMPS atom type number
+to the symbolic value of the internal potential parameter data structure.
+Following the *LJCutMelt* example, here are the two functions:
 
 .. code-block:: python
 
@@ -154,10 +154,10 @@ the *LJCutMelt* example, here are the two functions:
 
    for consistency with the C++ pair styles in LAMMPS, the
    *compute_force* function follows the conventions of the Pair::single()
-   methods and does not return the full force, but the force scaled by
-   the distance between the two atoms, so this value only needs to be
-   multiplied by delta x, delta y, and delta z to conveniently obtain the
-   three components of the force vector between these two atoms.
+   methods and does not return the pairwise force directly, but the force
+   divided by the distance between the two atoms, so this value only needs
+   to be  multiplied by delta x, delta y, and delta z to conveniently obtain
+   the three components of the force vector between these two atoms.
 
 ----------
 
diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp
index f9280f4a82..147a63512a 100644
--- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp
+++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp
@@ -47,7 +47,7 @@ void BondFENENM::compute(int eflag, int vflag)
 {
   int i1, i2, n, type;
   double delx, dely, delz, ebond, fbond;
-  double rsq, r0sq, rlogarg, sr2, sr6;
+  double rsq, r0sq, rlogarg, sr6;
   double r;
 
   ebond = sr6 = 0.0;
diff --git a/src/EXTRA-MOLECULE/dihedral_table_cut.cpp b/src/EXTRA-MOLECULE/dihedral_table_cut.cpp
index b132104bdd..522eff8626 100644
--- a/src/EXTRA-MOLECULE/dihedral_table_cut.cpp
+++ b/src/EXTRA-MOLECULE/dihedral_table_cut.cpp
@@ -82,20 +82,14 @@ enum { //GSL status return codes.
   GSL_EBADLEN  = 19
 };
 
-
-
 // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
 //           with n control points at xa[], ya[], with parameters y2a[].
 //           The xa[] must be monotonically increasing and their
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splintD(double const *xa,
-                          double const *ya,
-                          double const *y2a,
-                          int n,
-                          double period,
-                          double x)
+static double cyc_splintD(double const *xa, double const *ya, double const *y2a,
+                          int n, double period, double x)
 {
   int klo = -1;
   int khi = n; // (not n-1)
@@ -490,8 +484,7 @@ void DihedralTableCut::coeff(int narg, char **arg)
   if (tb->ninput < 2)
     error->all(FLERR,"Invalid dihedral table length: {}",arg[5]);
   else if ((tb->ninput == 2) && (tabstyle == SPLINE))
-    error->all(FLERR,"Invalid dihedral spline table length: {} "
-                                 "(Try linear)",arg[5]);
+    error->all(FLERR,"Invalid dihedral spline table length: {} (Try linear)",arg[5]);
 
   // check for monotonicity
   for (int i=0; i < tb->ninput-1; i++) {
@@ -509,12 +502,10 @@ void DihedralTableCut::coeff(int narg, char **arg)
   double phihi = tb->phifile[tb->ninput-1];
   if (tb->use_degrees) {
     if ((phihi - philo) >= 360)
-      error->all(FLERR,"Dihedral table angle range must be < 360 "
-                                   "degrees ({})",arg[5]);
+      error->all(FLERR,"Dihedral table angle range must be < 360 degrees ({})",arg[5]);
   } else {
     if ((phihi - philo) >= MY_2PI)
-      error->all(FLERR,"Dihedral table angle range must be < 2*PI "
-                                   "radians ({})",arg[5]);
+      error->all(FLERR,"Dihedral table angle range must be < 2*PI radians ({})",arg[5]);
   }
 
   // convert phi from degrees to radians
@@ -532,9 +523,9 @@ void DihedralTableCut::coeff(int narg, char **arg)
   // We also want the angles to be sorted in increasing order.
   // This messy code fixes these problems with the user's data:
   {
-    double *phifile_tmp = new double [tb->ninput];  //temporary arrays
-    double *ffile_tmp = new double [tb->ninput];  //used for sorting
-    double *efile_tmp = new double [tb->ninput];
+    double *phifile_tmp = new double[tb->ninput];  //temporary arrays
+    double *ffile_tmp = new double[tb->ninput];  //used for sorting
+    double *efile_tmp = new double[tb->ninput];
 
     // After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
     // If so, find the discontinuity:
diff --git a/src/GPU/fix_nve_gpu.cpp b/src/GPU/fix_nve_gpu.cpp
index 6612b8f65d..9392953398 100644
--- a/src/GPU/fix_nve_gpu.cpp
+++ b/src/GPU/fix_nve_gpu.cpp
@@ -37,7 +37,7 @@ using namespace FixConst;
 FixNVEGPU::FixNVEGPU(LAMMPS *lmp, int narg, char **arg) :
   FixNVE(lmp, narg, arg)
 {
-  _dtfm = 0;
+  _dtfm = nullptr;
   _nlocal_max = 0;
 }
 
@@ -57,7 +57,11 @@ void FixNVEGPU::setup(int vflag)
     _respa_on = 1;
   else
     _respa_on = 0;
-  if (atom->ntypes > 1) reset_dt();
+
+  // ensure that _dtfm array is initialized if the group is not "all"
+  // or there is more than one atom type as that re-ordeted array is used for
+  // per-type/per-atom masses and group membership detection.
+  if ((igroup != 0) || (atom->ntypes > 1)) reset_dt();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp
index de40cb2568..46d439b879 100644
--- a/src/MANYBODY/pair_eam_cd.cpp
+++ b/src/MANYBODY/pair_eam_cd.cpp
@@ -20,9 +20,6 @@
 
 #include "pair_eam_cd.h"
 
-#include <cmath>
-
-#include <cstring>
 #include "atom.h"
 #include "force.h"
 #include "comm.h"
@@ -31,11 +28,11 @@
 #include "error.h"
 #include "tokenizer.h"
 
-
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
-#define ASSERT(cond)
 #define MAXLINE 1024        // This sets the maximum line length in EAM input files.
 
 PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion)
@@ -298,7 +295,7 @@ void PairEAMCD::compute(int eflag, int vflag)
     // It will be replaced by the concentration at site i if atom i is either A or B.
 
     double x_i = -1.0;
-    double D_i, h_prime_i;
+    double D_i = 0.0, h_prime_i;
 
     // This if-clause is only required for ternary alloys.
 
@@ -307,7 +304,6 @@ void PairEAMCD::compute(int eflag, int vflag)
       // Compute local concentration at site i.
 
       x_i = rhoB[i]/rho[i];
-      ASSERT(x_i >= 0 && x_i<=1.0);
 
       if (cdeamVersion == 1) {
 
@@ -317,8 +313,6 @@ void PairEAMCD::compute(int eflag, int vflag)
         D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
       } else if (cdeamVersion == 2) {
         D_i = D_values[i];
-      } else {
-        ASSERT(false);
       }
     }
 
@@ -354,14 +348,11 @@ void PairEAMCD::compute(int eflag, int vflag)
 
         // This code line is required for ternary alloy.
 
-        if (jtype == speciesA || jtype == speciesB) {
-          ASSERT(rho[i] != 0.0);
-          ASSERT(rho[j] != 0.0);
+        if ((jtype == speciesA || jtype == speciesB) && rho[j] != 0.0) {
 
           // Compute local concentration at site j.
 
           x_j = rhoB[j]/rho[j];
-          ASSERT(x_j >= 0 && x_j<=1.0);
 
           double D_j=0.0;
           if (cdeamVersion == 1) {
@@ -372,8 +363,6 @@ void PairEAMCD::compute(int eflag, int vflag)
             D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
           } else if (cdeamVersion == 2) {
             D_j = D_values[j];
-          } else {
-            ASSERT(false);
           }
           double t2 = -rhoB[j];
           if (itype == speciesB) t2 += rho[j];
@@ -422,8 +411,6 @@ void PairEAMCD::compute(int eflag, int vflag)
             // Calculate h(x_ij) polynomial function.
 
             h = evalH(x_ij);
-          } else {
-            ASSERT(false);
           }
           fpair += h * phip;
           phi *= h;
@@ -460,7 +447,8 @@ void PairEAMCD::coeff(int narg, char **arg)
   // Make sure the EAM file is a CD-EAM binary alloy.
 
   if (setfl->nelements < 2)
-    error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");
+    error->all(FLERR,"The EAM file must contain at least 2 elements to be "
+                    "used with the eam/cd pair style.");
 
   // Read in the coefficients of the h polynomial from the end of the EAM file.
 
@@ -502,22 +490,28 @@ void PairEAMCD::read_h_coeff(char *filename)
     // Open potential file
 
     FILE *fptr;
-    char line[MAXLINE];
-    char nextline[MAXLINE];
     int convert_flag = unit_convert_flag;
     fptr = utils::open_potential(filename, lmp, &convert_flag);
     if (fptr == nullptr)
-      error->one(FLERR,"Cannot open EAMCD potential file {}",
-                                   filename);
+      error->one(FLERR,"Cannot open EAMCD potential file {}", filename);
 
     // h coefficients are stored at the end of the file.
-    // Skip to last line of file.
+    // Seek to end of file, read last part into a buffer and
+    // then skip over lines in buffer until reaching the end.
 
-    while (fgets(nextline, MAXLINE, fptr) != nullptr) {
-      strcpy(line, nextline);
-    }
+    platform::fseek(fptr, platform::END_OF_FILE);
+    platform::fseek(fptr, platform::ftell(fptr) - MAXLINE);
+    char *buf = new char[MAXLINE+1];
+    fread(buf, 1, MAXLINE, fptr);
+    buf[MAXLINE] = '\0';        // must 0-terminate buffer for string processing
+    Tokenizer lines(buf, "\n");
+    delete[] buf;
 
-    ValueTokenizer values(line);
+    std::string lastline;
+    while (lines.has_next())
+      lastline = lines.next();
+
+    ValueTokenizer values(lastline);
     int degree = values.next_int();
     nhcoeff = degree+1;
 
@@ -527,10 +521,8 @@ void PairEAMCD::read_h_coeff(char *filename)
     delete[] hcoeff;
     hcoeff = new double[nhcoeff];
 
-    int i = 0;
-    while (values.has_next()) {
-      hcoeff[i++] = values.next_double();
-    }
+    for (int i = 0; i < nhcoeff; ++i)
+      hcoeff[i] = values.next_double();
 
     // Close the potential file.
 
@@ -545,7 +537,6 @@ void PairEAMCD::read_h_coeff(char *filename)
   MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
@@ -572,7 +563,7 @@ int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
         buf[m++] = rhoB[j];
       }
       return m;
-    } else { ASSERT(false); return 0; }
+    } else return 0;
   } else if (communicationStage == 4) {
     for (i = 0; i < n; i++) {
       j = list[i];
@@ -604,8 +595,6 @@ void PairEAMCD::unpack_forward_comm(int n, int first, double *buf)
         rho[i] = buf[m++];
         rhoB[i] = buf[m++];
       }
-    } else {
-      ASSERT(false);
     }
   } else if (communicationStage == 4) {
     for (i = first; i < last; i++) {
@@ -636,7 +625,7 @@ int PairEAMCD::pack_reverse_comm(int n, int first, double *buf)
         buf[m++] = rhoB[i];
       }
       return m;
-    } else { ASSERT(false); return 0; }
+    } else return 0;
   } else if (communicationStage == 3) {
     for (i = first; i < last; i++) {
       buf[m++] = D_values[i];
@@ -666,8 +655,6 @@ void PairEAMCD::unpack_reverse_comm(int n, int *list, double *buf)
         rho[j] += buf[m++];
         rhoB[j] += buf[m++];
       }
-    } else {
-      ASSERT(false);
     }
   } else if (communicationStage == 3) {
     for (i = 0; i < n; i++) {
diff --git a/src/MANYBODY/pair_eam_cd.h b/src/MANYBODY/pair_eam_cd.h
index 77e909f48b..6846a6cd76 100644
--- a/src/MANYBODY/pair_eam_cd.h
+++ b/src/MANYBODY/pair_eam_cd.h
@@ -120,6 +120,7 @@ class PairEAMCD : public PairEAMAlloy {
     index.p = r * rdr + 1.0;
     index.m = static_cast<int>(index.p);
     index.m = index.m <= (nr - 1) ? index.m : (nr - 1);
+    index.m = index.m > 1 ? index.m : 1;
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;
@@ -132,6 +133,7 @@ class PairEAMCD : public PairEAMAlloy {
     index.p = rho * rdrho + 1.0;
     index.m = static_cast<int>(index.p);
     index.m = index.m <= (nrho - 1) ? index.m : (nrho - 1);
+    index.m = index.m > 1 ? index.m : 1;
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;
diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp
index a5ba07b779..2d7356db01 100644
--- a/src/MOLECULE/angle_table.cpp
+++ b/src/MOLECULE/angle_table.cpp
@@ -470,9 +470,9 @@ void AngleTable::compute_table(Table *tb)
 
   memory->create(tb->ang,tablength,"angle:ang");
   memory->create(tb->e,tablength,"angle:e");
-  memory->create(tb->de,tlm1,"angle:de");
+  memory->create(tb->de,tablength,"angle:de");
   memory->create(tb->f,tablength,"angle:f");
-  memory->create(tb->df,tlm1,"angle:df");
+  memory->create(tb->df,tablength,"angle:df");
   memory->create(tb->e2,tablength,"angle:e2");
   memory->create(tb->f2,tablength,"angle:f2");
 
@@ -488,6 +488,9 @@ void AngleTable::compute_table(Table *tb)
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
+  // get final elements from linear extrapolation
+  tb->de[tlm1] = 2.0*tb->de[tlm1-1] - tb->de[tlm1-2];
+  tb->df[tlm1] = 2.0*tb->df[tlm1-1] - tb->df[tlm1-2];
 
   double ep0 = - tb->f[0];
   double epn = - tb->f[tlm1];
@@ -575,7 +578,7 @@ void AngleTable::spline(double *x, double *y, int n,
   double p,qn,sig,un;
   double *u = new double[n];
 
-  if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
+  if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
   else {
     y2[0] = -0.5;
     u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@@ -587,7 +590,7 @@ void AngleTable::spline(double *x, double *y, int n,
     u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
     u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
   }
-  if (ypn > 0.99e30) qn = un = 0.0;
+  if (ypn > 0.99e300) qn = un = 0.0;
   else {
     qn = 0.5;
     un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
@@ -615,8 +618,7 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
   h = xa[khi]-xa[klo];
   a = (xa[khi]-x) / h;
   b = (x-xa[klo]) / h;
-  y = a*ya[klo] + b*ya[khi] +
-    ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+  y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
   return y;
 }
 
@@ -632,8 +634,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
 
   double fraction,a,b;
   const Table *tb = &tables[tabindex[type]];
-  int itable = static_cast<int> (x * tb->invdelta);
 
+  // invdelta is based on tablength-1
+  int itable = static_cast<int> (x * tb->invdelta);
   if (itable < 0) itable = 0;
   if (itable >= tablength) itable = tablength-1;
 
@@ -647,11 +650,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
     b = (x - tb->ang[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
     f = a * tb->f[itable] + b * tb->f[itable+1] +
-      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6;
   }
 }
 
@@ -681,7 +682,6 @@ void AngleTable::u_lookup(int type, double x, double &u)
     b = (x - tb->ang[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
   }
 }
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index db1314c76f..48d7377682 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -435,9 +435,9 @@ void BondTable::compute_table(Table *tb)
 
   memory->create(tb->r,tablength,"bond:r");
   memory->create(tb->e,tablength,"bond:e");
-  memory->create(tb->de,tlm1,"bond:de");
+  memory->create(tb->de,tablength,"bond:de");
   memory->create(tb->f,tablength,"bond:f");
-  memory->create(tb->df,tlm1,"bond:df");
+  memory->create(tb->df,tablength,"bond:df");
   memory->create(tb->e2,tablength,"bond:e2");
   memory->create(tb->f2,tablength,"bond:f2");
 
@@ -453,6 +453,9 @@ void BondTable::compute_table(Table *tb)
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
+  // get final elements from linear extrapolation
+  tb->de[tlm1] = 2.0*tb->de[tlm1-1] - tb->de[tlm1-2];
+  tb->df[tlm1] = 2.0*tb->df[tlm1-1] - tb->df[tlm1-2];
 
   double ep0 = - tb->f[0];
   double epn = - tb->f[tlm1];
@@ -538,7 +541,7 @@ void BondTable::spline(double *x, double *y, int n,
   double p,qn,sig,un;
   double *u = new double[n];
 
-  if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
+  if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
   else {
     y2[0] = -0.5;
     u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@@ -550,7 +553,7 @@ void BondTable::spline(double *x, double *y, int n,
     u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
     u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
   }
-  if (ypn > 0.99e30) qn = un = 0.0;
+  if (ypn > 0.99e300) qn = un = 0.0;
   else {
     qn = 0.5;
     un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
@@ -578,8 +581,7 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x)
   h = xa[khi]-xa[klo];
   a = (xa[khi]-x) / h;
   b = (x-xa[klo]) / h;
-  y = a*ya[klo] + b*ya[khi] +
-    ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+  y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
   return y;
 }
 
@@ -598,11 +600,9 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
   const Table *tb = &tables[tabindex[type]];
   const int itable = static_cast<int> ((x - tb->lo) * tb->invdelta);
   if (itable < 0)
-    error->one(FLERR,"Bond length < table inner cutoff: "
-               "type {} length {:.8}",type,x);
+    error->one(FLERR,"Bond length < table inner cutoff: type {} length {:.8}",type,x);
   else if (itable >= tablength)
-    error->one(FLERR,"Bond length > table outer cutoff: "
-               "type {} length {:.8}",type,x);
+    error->one(FLERR,"Bond length > table outer cutoff: type {} length {:.8}",type,x);
 
   if (tabstyle == LINEAR) {
     fraction = (x - tb->r[itable]) * tb->invdelta;
@@ -614,10 +614,8 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
     b = (x - tb->r[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
     f = a * tb->f[itable] + b * tb->f[itable+1] +
-      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6;
   }
 }
diff --git a/src/MOLECULE/dihedral_table.cpp b/src/MOLECULE/dihedral_table.cpp
index a91324dd98..dbca4a85c1 100644
--- a/src/MOLECULE/dihedral_table.cpp
+++ b/src/MOLECULE/dihedral_table.cpp
@@ -189,11 +189,8 @@ static int solve_cyc_tridiag( const double diag[], size_t d_stride,
    spline and splint routines modified from Numerical Recipes
 ------------------------------------------------------------------------- */
 
-static int cyc_spline(double const *xa,
-                      double const *ya,
-                      int n,
-                      double period,
-                      double *y2a, bool warn)
+static int cyc_spline(double const *xa, double const *ya, int n,
+                      double period, double *y2a, bool warn)
 {
 
   double *diag    = new double[n];
@@ -264,12 +261,8 @@ static int cyc_spline(double const *xa,
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splint(double const *xa,
-                         double const *ya,
-                         double const *y2a,
-                         int n,
-                         double period,
-                         double x)
+static double cyc_splint(double const *xa, double const *ya, double const *y2a,
+                         int n, double period, double x)
 {
   int klo = -1;
   int khi = n;
@@ -302,11 +295,8 @@ static double cyc_splint(double const *xa,
 } // cyc_splint()
 
 
-static double cyc_lin(double const *xa,
-                      double const *ya,
-                      int n,
-                      double period,
-                      double x)
+static double cyc_lin(double const *xa, double const *ya,
+                      int n, double period, double x)
 {
   int klo = -1;
   int khi = n;
@@ -337,21 +327,14 @@ static double cyc_lin(double const *xa,
 
 } // cyc_lin()
 
-
-
-
 // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
 //           with n control points at xa[], ya[], with parameters y2a[].
 //           The xa[] must be monotonically increasing and their
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splintD(double const *xa,
-                          double const *ya,
-                          double const *y2a,
-                          int n,
-                          double period,
-                          double x)
+static double cyc_splintD(double const *xa, double const *ya, double const *y2a,
+                          int n, double period, double x)
 {
   int klo = -1;
   int khi = n; // (not n-1)
@@ -829,9 +812,9 @@ void DihedralTable::coeff(int narg, char **arg)
   // We also want the angles to be sorted in increasing order.
   // This messy code fixes these problems with the user's data:
   {
-    double *phifile_tmp = new double [tb->ninput];  //temporary arrays
-    double *ffile_tmp = new double [tb->ninput];  //used for sorting
-    double *efile_tmp = new double [tb->ninput];
+    double *phifile_tmp = new double[tb->ninput];  //temporary arrays
+    double *ffile_tmp = new double[tb->ninput];  //used for sorting
+    double *efile_tmp = new double[tb->ninput];
 
     // After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
     // If so, find the discontinuity:
@@ -1184,8 +1167,7 @@ void DihedralTable::compute_table(Table *tb)
       if (! tb->f_unspecified)
         tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi);
     }
-  } // if (tabstyle == SPLINE)
-  else if (tabstyle == LINEAR) {
+  } else if (tabstyle == LINEAR) {
     if (! tb->f_unspecified) {
       for (int i = 0; i < tablength; i++) {
         double phi = i*tb->delta;
@@ -1193,8 +1175,7 @@ void DihedralTable::compute_table(Table *tb)
         tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
         tb->f[i]= cyc_lin(tb->phifile,tb->ffile,tb->ninput,MY_2PI,phi);
       }
-    }
-    else {
+    } else {
       for (int i = 0; i < tablength; i++) {
         double phi = i*tb->delta;
         tb->phi[i] = phi;
@@ -1269,8 +1250,7 @@ void DihedralTable::param_extract(Table *tb, char *line)
       //else if (word == "EQ") {
       //  tb->theta0 = values.next_double();
       //}
-      else error->one(FLERR,"Invalid keyword in dihedral angle "
-                                        "table parameters ({})", word);
+      else error->one(FLERR,"Invalid keyword in dihedral angle table parameters ({})", word);
     }
   } catch (TokenizerException &e) {
     error->one(FLERR, e.what());
diff --git a/src/MOLFILE/molfile_interface.cpp b/src/MOLFILE/molfile_interface.cpp
index 9ce3822082..5fd398570e 100644
--- a/src/MOLFILE/molfile_interface.cpp
+++ b/src/MOLFILE/molfile_interface.cpp
@@ -75,25 +75,25 @@ extern "C" {
 
   /* corresponding table of masses. */
   static const float pte_mass[] = {
-    /* X  */ 0.00000, 1.00794, 4.00260, 6.941, 9.012182, 10.811,
-    /* C  */ 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797,
-    /* Na */ 22.989770, 24.3050, 26.981538, 28.0855, 30.973761,
-    /* S  */ 32.065, 35.453, 39.948, 39.0983, 40.078, 44.955910,
-    /* Ti */ 47.867, 50.9415, 51.9961, 54.938049, 55.845, 58.9332,
-    /* Ni */ 58.6934, 63.546, 65.409, 69.723, 72.64, 74.92160,
-    /* Se */ 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585,
-    /* Zr */ 91.224, 92.90638, 95.94, 98.0, 101.07, 102.90550,
-    /* Pd */ 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760,
-    /* Te */ 127.60, 126.90447, 131.293, 132.90545, 137.327,
-    /* La */ 138.9055, 140.116, 140.90765, 144.24, 145.0, 150.36,
-    /* Eu */ 151.964, 157.25, 158.92534, 162.500, 164.93032,
-    /* Er */ 167.259, 168.93421, 173.04, 174.967, 178.49, 180.9479,
-    /* W  */ 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655,
-    /* Hg */ 200.59, 204.3833, 207.2, 208.98038, 209.0, 210.0, 222.0,
-    /* Fr */ 223.0, 226.0, 227.0, 232.0381, 231.03588, 238.02891,
-    /* Np */ 237.0, 244.0, 243.0, 247.0, 247.0, 251.0, 252.0, 257.0,
-    /* Md */ 258.0, 259.0, 262.0, 261.0, 262.0, 266.0, 264.0, 269.0,
-    /* Mt */ 268.0, 271.0, 272.0
+    /* X  */ 0.00000f, 1.00794f, 4.00260f, 6.941f, 9.012182f, 10.811f,
+    /* C  */ 12.0107f, 14.0067f, 15.9994f, 18.9984032f, 20.1797f,
+    /* Na */ 22.989770f, 24.3050f, 26.981538f, 28.0855f, 30.973761f,
+    /* S  */ 32.065f, 35.453f, 39.948f, 39.0983f, 40.078f, 44.955910f,
+    /* Ti */ 47.867f, 50.9415f, 51.9961f, 54.938049f, 55.845f, 58.9332f,
+    /* Ni */ 58.6934f, 63.546f, 65.409f, 69.723f, 72.64f, 74.92160f,
+    /* Se */ 78.96f, 79.904f, 83.798f, 85.4678f, 87.62f, 88.90585f,
+    /* Zr */ 91.224f, 92.90638f, 95.94f, 98.0f, 101.07f, 102.90550f,
+    /* Pd */ 106.42f, 107.8682f, 112.411f, 114.818f, 118.710f, 121.760f,
+    /* Te */ 127.60f, 126.90447f, 131.293f, 132.90545f, 137.327f,
+    /* La */ 138.9055f, 140.116f, 140.90765f, 144.24f, 145.0f, 150.36f,
+    /* Eu */ 151.964f, 157.25f, 158.92534f, 162.500f, 164.93032f,
+    /* Er */ 167.259f, 168.93421f, 173.04f, 174.967f, 178.49f, 180.9479f,
+    /* W  */ 183.84f, 186.207f, 190.23f, 192.217f, 195.078f, 196.96655f,
+    /* Hg */ 200.59f, 204.3833f, 207.2f, 208.98038f, 209.0f, 210.0f, 222.0f,
+    /* Fr */ 223.0f, 226.0f, 227.0f, 232.0381f, 231.03588f, 238.02891f,
+    /* Np */ 237.0f, 244.0f, 243.0f, 247.0f, 247.0f, 251.0f, 252.0f, 257.0f,
+    /* Md */ 258.0f, 259.0f, 262.0f, 261.0f, 262.0f, 266.0f, 264.0f, 269.0f,
+    /* Mt */ 268.0f, 271.0f, 272.0f
   };
 
   /*
@@ -107,25 +107,25 @@ extern "C" {
    * Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES) by default.
    */
   static const float pte_vdw_radius[] = {
-    /* X  */ 1.5, 1.2, 1.4, 1.82, 2.0, 2.0,
-    /* C  */ 1.7, 1.55, 1.52, 1.47, 1.54,
-    /* Na */ 1.36, 1.18, 2.0, 2.1, 1.8,
-    /* S  */ 1.8, 2.27, 1.88, 1.76, 1.37, 2.0,
-    /* Ti */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Ni */ 1.63, 1.4, 1.39, 1.07, 2.0, 1.85,
-    /* Se */ 1.9, 1.85, 2.02, 2.0, 2.0, 2.0,
-    /* Zr */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Pd */ 1.63, 1.72, 1.58, 1.93, 2.17, 2.0,
-    /* Te */ 2.06, 1.98, 2.16, 2.1, 2.0,
-    /* La */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Eu */ 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Er */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* W  */ 2.0, 2.0, 2.0, 2.0, 1.72, 1.66,
-    /* Hg */ 1.55, 1.96, 2.02, 2.0, 2.0, 2.0, 2.0,
-    /* Fr */ 2.0, 2.0, 2.0, 2.0, 2.0, 1.86,
-    /* Np */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Md */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Mt */ 2.0, 2.0, 2.0
+    /* X  */ 1.5f, 1.2f, 1.4f, 1.82f, 2.0f, 2.0f,
+    /* C  */ 1.7f, 1.55f, 1.52f, 1.47f, 1.54f,
+    /* Na */ 1.36f, 1.18f, 2.0f, 2.1f, 1.8f,
+    /* S  */ 1.8f, 2.27f, 1.88f, 1.76f, 1.37f, 2.0f,
+    /* Ti */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Ni */ 1.63f, 1.4f, 1.39f, 1.07f, 2.0f, 1.85f,
+    /* Se */ 1.9f, 1.85f, 2.02f, 2.0f, 2.0f, 2.0f,
+    /* Zr */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Pd */ 1.63f, 1.72f, 1.58f, 1.93f, 2.17f, 2.0f,
+    /* Te */ 2.06f, 1.98f, 2.16f, 2.1f, 2.0f,
+    /* La */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Eu */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Er */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* W  */ 2.0f, 2.0f, 2.0f, 2.0f, 1.72f, 1.66f,
+    /* Hg */ 1.55f, 1.96f, 2.02f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Fr */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 1.86f,
+    /* Np */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Md */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Mt */ 2.0f, 2.0f, 2.0f
   };
 
   /* lookup functions */
diff --git a/src/OPENMP/angle_table_omp.cpp b/src/OPENMP/angle_table_omp.cpp
index 892f9295a5..cca34a67f7 100644
--- a/src/OPENMP/angle_table_omp.cpp
+++ b/src/OPENMP/angle_table_omp.cpp
@@ -16,15 +16,16 @@
    Contributing author: Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
-#include "omp_compat.h"
 #include "angle_table_omp.h"
-#include <cmath>
+
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 
+#include <cmath>
 
+#include "omp_compat.h"
 #include "suffix.h"
 using namespace LAMMPS_NS;
 
diff --git a/src/OPENMP/bond_table_omp.cpp b/src/OPENMP/bond_table_omp.cpp
index faadca456a..dcc13c85c9 100644
--- a/src/OPENMP/bond_table_omp.cpp
+++ b/src/OPENMP/bond_table_omp.cpp
@@ -16,16 +16,16 @@
    Contributing author: Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
-#include "omp_compat.h"
 #include "bond_table_omp.h"
+
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 
-
 #include <cmath>
 
+#include "omp_compat.h"
 #include "suffix.h"
 using namespace LAMMPS_NS;
 
diff --git a/src/error.cpp b/src/error.cpp
index 5338f41cca..e2162cf661 100644
--- a/src/error.cpp
+++ b/src/error.cpp
@@ -196,7 +196,7 @@ void Error::one(const std::string &file, int line, const std::string &str)
   MPI_Comm_rank(world,&me);
 
   if (input && input->line) lastcmd = input->line;
-  std::string mesg = fmt::format("ERROR on proc {}: {} ({}:{})\n",
+  std::string mesg = fmt::format("ERROR on proc {}: {} ({}:{})\nLast command: {}\n",
                                  me,str,truncpath(file),line,lastcmd);
   utils::logmesg(lmp,mesg);
 
diff --git a/src/output.cpp b/src/output.cpp
index 4cea582283..a8fda4173d 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -344,8 +344,6 @@ void Output::setup(int memflag)
    //   can't remove an uneeded addstep from a compute, b/c don't know
    //     what other command may have added it
 
-   int writeflag;
-
    if (next_dump_any == ntimestep) {
 
      for (int idump = 0; idump < ndump; idump++) {