From e02ad1a3b208797cc9c5b3a08c24cabb5065b7a6 Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 17:02:22 -0600
Subject: [PATCH 1/7] Cleaning up SPH package, minor bug fixes

---
 src/SPH/compute_sph_e_atom.cpp        | 30 ++++++-------
 src/SPH/compute_sph_rho_atom.cpp      | 33 +++++++-------
 src/SPH/compute_sph_t_atom.cpp        | 32 +++++++-------
 src/SPH/fix_sph.cpp                   | 32 ++++++++------
 src/SPH/fix_sph_stationary.cpp        | 27 ++++++------
 src/SPH/fix_sph_stationary.h          |  3 --
 src/SPH/pair_sph_heatconduction.cpp   | 33 ++++++++------
 src/SPH/pair_sph_idealgas.cpp         | 46 +++++++++-----------
 src/SPH/pair_sph_idealgas.h           |  1 -
 src/SPH/pair_sph_lj.cpp               | 46 ++++++++++----------
 src/SPH/pair_sph_lj.h                 |  2 -
 src/SPH/pair_sph_rhosum.cpp           | 56 ++++++++++--------------
 src/SPH/pair_sph_taitwater.cpp        | 58 +++++++++----------------
 src/SPH/pair_sph_taitwater_morris.cpp | 62 +++++++++++----------------
 14 files changed, 215 insertions(+), 246 deletions(-)

diff --git a/src/SPH/compute_sph_e_atom.cpp b/src/SPH/compute_sph_e_atom.cpp
index 05e752c49bb..8bb0622acd2 100644
--- a/src/SPH/compute_sph_e_atom.cpp
+++ b/src/SPH/compute_sph_e_atom.cpp
@@ -13,13 +13,15 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_sph_e_atom.h"
-#include <cstring>
+
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -31,7 +33,7 @@ ComputeSPHEAtom::ComputeSPHEAtom(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3)
     error->all(FLERR,"Number of arguments for compute sph/e/atom command != 3");
   if (atom->esph_flag != 1)
-    error->all(FLERR,"Compute sph/e/atom command requires atom_style sph)");
+    error->all(FLERR,"Compute sph/e/atom requires atom attribut energy, e.g. in atom_style sph");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -51,12 +53,11 @@ ComputeSPHEAtom::~ComputeSPHEAtom()
 
 void ComputeSPHEAtom::init()
 {
-
   int count = 0;
   for (int i = 0; i < modify->ncompute; i++)
-    if (strcmp(modify->compute[i]->style,"evector/atom") == 0) count++;
+    if (strcmp(modify->compute[i]->style,"sph/e/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one compute evector/atom");
+    error->warning(FLERR,"More than one compute sph/e/atom");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -78,14 +79,13 @@ void ComputeSPHEAtom::compute_peratom()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-              evector[i] = esph[i];
-      }
-      else {
-              evector[i] = 0.0;
-      }
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      evector[i] = esph[i];
+    } else {
+      evector[i] = 0.0;
     }
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/SPH/compute_sph_rho_atom.cpp b/src/SPH/compute_sph_rho_atom.cpp
index 9c596bdfc30..54d12373689 100644
--- a/src/SPH/compute_sph_rho_atom.cpp
+++ b/src/SPH/compute_sph_rho_atom.cpp
@@ -13,13 +13,15 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_sph_rho_atom.h"
-#include <cstring>
+
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -30,8 +32,7 @@ ComputeSPHRhoAtom::ComputeSPHRhoAtom(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 3) error->all(FLERR,"Illegal compute sph/rho/atom command");
   if (atom->rho_flag != 1)
-    error->all(FLERR,"Compute sph/rho/atom command requires atom_style sph");
-
+    error->all(FLERR, "Compute sph/rho/atom requires atom attribute density, e.g. in atom_style sph");
   peratom_flag = 1;
   size_peratom_cols = 0;
 
@@ -50,12 +51,11 @@ ComputeSPHRhoAtom::~ComputeSPHRhoAtom()
 
 void ComputeSPHRhoAtom::init()
 {
-
   int count = 0;
   for (int i = 0; i < modify->ncompute; i++)
-    if (strcmp(modify->compute[i]->style,"rhoVector/atom") == 0) count++;
+    if (strcmp(modify->compute[i]->style,"sph/rho/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one compute rhoVector/atom");
+    error->warning(FLERR,"More than one compute sph/rho/atom");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -73,20 +73,17 @@ void ComputeSPHRhoAtom::compute_peratom()
     vector_atom = rhoVector;
   }
 
-  // compute kinetic energy for each atom in group
-
   double *rho = atom->rho;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-              rhoVector[i] = rho[i];
-      }
-      else {
-              rhoVector[i] = 0.0;
-      }
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      rhoVector[i] = rho[i];
+    } else {
+      rhoVector[i] = 0.0;
     }
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/SPH/compute_sph_t_atom.cpp b/src/SPH/compute_sph_t_atom.cpp
index 4a830134346..e795c32a5ea 100644
--- a/src/SPH/compute_sph_t_atom.cpp
+++ b/src/SPH/compute_sph_t_atom.cpp
@@ -13,13 +13,15 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_sph_t_atom.h"
-#include <cstring>
+
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -31,7 +33,7 @@ ComputeSPHTAtom::ComputeSPHTAtom(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3)
     error->all(FLERR,"Number of arguments for compute sph/t/atom command != 3");
   if ((atom->esph_flag != 1) || (atom->cv_flag != 1))
-    error->all(FLERR,"Compute sph/t/atom command requires atom_style sph");
+    error->all(FLERR, "Compute sph/t/atom requires atom attributes energy and specific heat, e.g. in atom_style sph");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -51,12 +53,11 @@ ComputeSPHTAtom::~ComputeSPHTAtom()
 
 void ComputeSPHTAtom::init()
 {
-
   int count = 0;
   for (int i = 0; i < modify->ncompute; i++)
-    if (strcmp(modify->compute[i]->style,"meso/t/atom") == 0) count++;
+    if (strcmp(modify->compute[i]->style,"sph/t/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one compute meso/t/atom");
+    error->warning(FLERR,"More than one compute sph/t/atom");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -79,16 +80,15 @@ void ComputeSPHTAtom::compute_peratom()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-              if (cv[i] > 0.0) {
-                      tvector[i] = esph[i] / cv[i];
-              }
-      }
-      else {
-              tvector[i] = 0.0;
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (cv[i] > 0.0) {
+        tvector[i] = esph[i] / cv[i];
       }
+    } else {
+      tvector[i] = 0.0;
     }
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/SPH/fix_sph.cpp b/src/SPH/fix_sph.cpp
index 0dafcee7c00..96479be9281 100644
--- a/src/SPH/fix_sph.cpp
+++ b/src/SPH/fix_sph.cpp
@@ -13,10 +13,13 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_sph.h"
+
 #include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "error.h"
 #include "force.h"
 #include "update.h"
-#include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -24,11 +27,10 @@ using namespace FixConst;
 /* ---------------------------------------------------------------------- */
 
 FixSPH::FixSPH(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg) {
-
-  if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
-    error->all(FLERR,
-        "Fix sph command requires atom_style with both energy and density");
+  Fix(lmp, narg, arg)
+{
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1))
+    error->all(FLERR, "Fix sph requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph");
 
   if (narg != 3)
     error->all(FLERR,"Illegal number of arguments for fix sph command");
@@ -38,7 +40,8 @@ FixSPH::FixSPH(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-int FixSPH::setmask() {
+int FixSPH::setmask()
+{
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= FINAL_INTEGRATE;
@@ -48,11 +51,14 @@ int FixSPH::setmask() {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPH::init() {
+void FixSPH::init()
+{
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void FixSPH::setup_pre_force(int /*vflag*/)
 {
   // set vest equal to v
@@ -76,7 +82,8 @@ void FixSPH::setup_pre_force(int /*vflag*/)
  allow for both per-type and per-atom mass
  ------------------------------------------------------------------------- */
 
-void FixSPH::initial_integrate(int /*vflag*/) {
+void FixSPH::initial_integrate(int /*vflag*/)
+{
   // update v and x and rho and e of atoms in group
 
   double **x = atom->x;
@@ -129,8 +136,8 @@ void FixSPH::initial_integrate(int /*vflag*/) {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPH::final_integrate() {
-
+void FixSPH::final_integrate()
+{
   // update v, rho, and e of atoms in group
 
   double **v = atom->v;
@@ -169,7 +176,8 @@ void FixSPH::final_integrate() {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPH::reset_dt() {
+void FixSPH::reset_dt()
+{
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 }
diff --git a/src/SPH/fix_sph_stationary.cpp b/src/SPH/fix_sph_stationary.cpp
index 673a7d50bfa..877bb5b17b7 100644
--- a/src/SPH/fix_sph_stationary.cpp
+++ b/src/SPH/fix_sph_stationary.cpp
@@ -13,10 +13,11 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_sph_stationary.h"
+
 #include "atom.h"
+#include "error.h"
 #include "force.h"
 #include "update.h"
-#include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -24,11 +25,10 @@ using namespace FixConst;
 /* ---------------------------------------------------------------------- */
 
 FixSPHStationary::FixSPHStationary(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg) {
-
+  Fix(lmp, narg, arg)
+{
   if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
-    error->all(FLERR,
-        "Fix sph/stationary command requires atom_style with both energy and density, e.g. meso");
+    error->all(FLERR, "Fix sph/stationary requires atom attributes energy and density, e.g. in atom_style sph");
 
   if (narg != 3)
     error->all(FLERR,"Illegal number of arguments for fix sph/stationary command");
@@ -38,7 +38,8 @@ FixSPHStationary::FixSPHStationary(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-int FixSPHStationary::setmask() {
+int FixSPHStationary::setmask()
+{
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= FINAL_INTEGRATE;
@@ -47,7 +48,8 @@ int FixSPHStationary::setmask() {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPHStationary::init() {
+void FixSPHStationary::init()
+{
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 }
@@ -56,8 +58,8 @@ void FixSPHStationary::init() {
  allow for both per-type and per-atom mass
  ------------------------------------------------------------------------- */
 
-void FixSPHStationary::initial_integrate(int /*vflag*/) {
-
+void FixSPHStationary::initial_integrate(int /*vflag*/)
+{
   double *rho = atom->rho;
   double *drho = atom->drho;
   double *esph = atom->esph;
@@ -80,8 +82,8 @@ void FixSPHStationary::initial_integrate(int /*vflag*/) {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPHStationary::final_integrate() {
-
+void FixSPHStationary::final_integrate()
+{
   double *esph = atom->esph;
   double *desph = atom->desph;
   double *rho = atom->rho;
@@ -101,7 +103,8 @@ void FixSPHStationary::final_integrate() {
 
 /* ---------------------------------------------------------------------- */
 
-void FixSPHStationary::reset_dt() {
+void FixSPHStationary::reset_dt()
+{
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 }
diff --git a/src/SPH/fix_sph_stationary.h b/src/SPH/fix_sph_stationary.h
index 51e336a038c..78ece5e5d4b 100644
--- a/src/SPH/fix_sph_stationary.h
+++ b/src/SPH/fix_sph_stationary.h
@@ -33,9 +33,6 @@ class FixSPHStationary : public Fix {
   void final_integrate() override;
   void reset_dt() override;
 
- private:
-  class NeighList *list;
-
  protected:
   double dtv, dtf;
   double *step_respa;
diff --git a/src/SPH/pair_sph_heatconduction.cpp b/src/SPH/pair_sph_heatconduction.cpp
index cb27982fa3e..145d9cadeec 100644
--- a/src/SPH/pair_sph_heatconduction.cpp
+++ b/src/SPH/pair_sph_heatconduction.cpp
@@ -13,14 +13,15 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_heatconduction.h"
-#include <cmath>
+
 #include "atom.h"
+#include "domain.h"
+#include "error.h"
 #include "force.h"
 #include "memory.h"
-#include "error.h"
 #include "neigh_list.h"
-#include "domain.h"
 
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -28,12 +29,16 @@ using namespace LAMMPS_NS;
 
 PairSPHHeatConduction::PairSPHHeatConduction(LAMMPS *lmp) : Pair(lmp)
 {
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
+    error->all(FLERR, "Pair sph/heatconduction requires atom attributes energy and density, e.g. in atom_style sph");
+
   restartinfo = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHHeatConduction::~PairSPHHeatConduction() {
+PairSPHHeatConduction::~PairSPHHeatConduction()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -44,7 +49,8 @@ PairSPHHeatConduction::~PairSPHHeatConduction() {
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHHeatConduction::compute(int eflag, int vflag) {
+void PairSPHHeatConduction::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, inum, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz;
 
@@ -124,7 +130,6 @@ void PairSPHHeatConduction::compute(int eflag, int vflag) {
         if (newton_pair || j < nlocal) {
           desph[j] -= deltaE;
         }
-
       }
     }
   }
@@ -134,7 +139,8 @@ void PairSPHHeatConduction::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHHeatConduction::allocate() {
+void PairSPHHeatConduction::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -152,7 +158,8 @@ void PairSPHHeatConduction::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHHeatConduction::settings(int narg, char **/*arg*/) {
+void PairSPHHeatConduction::settings(int narg, char **/*arg*/)
+{
   if (narg != 0)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/heatconduction");
@@ -162,7 +169,8 @@ void PairSPHHeatConduction::settings(int narg, char **/*arg*/) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHHeatConduction::coeff(int narg, char **arg) {
+void PairSPHHeatConduction::coeff(int narg, char **arg)
+{
   if (narg != 4)
     error->all(FLERR,"Incorrect number of args for pair_style sph/heatconduction coefficients");
   if (!allocated)
@@ -178,7 +186,6 @@ void PairSPHHeatConduction::coeff(int narg, char **arg) {
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
-      //printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
       alpha[i][j] = alpha_one;
       setflag[i][j] = 1;
@@ -194,7 +201,8 @@ void PairSPHHeatConduction::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHHeatConduction::init_one(int i, int j) {
+double PairSPHHeatConduction::init_one(int i, int j)
+{
 
   if (setflag[i][j] == 0) {
     error->all(FLERR,"All pair sph/heatconduction coeffs are not set");
@@ -209,7 +217,8 @@ double PairSPHHeatConduction::init_one(int i, int j) {
 /* ---------------------------------------------------------------------- */
 
 double PairSPHHeatConduction::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/,
-    double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) {
+    double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce)
+{
   fforce = 0.0;
 
   return 0.0;
diff --git a/src/SPH/pair_sph_idealgas.cpp b/src/SPH/pair_sph_idealgas.cpp
index 114d323bbd3..97cbda2ab34 100644
--- a/src/SPH/pair_sph_idealgas.cpp
+++ b/src/SPH/pair_sph_idealgas.cpp
@@ -13,14 +13,15 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_idealgas.h"
-#include <cmath>
+
 #include "atom.h"
+#include "domain.h"
+#include "error.h"
 #include "force.h"
-#include "neigh_list.h"
 #include "memory.h"
-#include "error.h"
-#include "domain.h"
+#include "neigh_list.h"
 
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -28,12 +29,17 @@ using namespace LAMMPS_NS;
 
 PairSPHIdealGas::PairSPHIdealGas(LAMMPS *lmp) : Pair(lmp)
 {
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1))
+    error->all(FLERR, "Pair sph/idealgas requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph");
+
   restartinfo = 0;
+  single_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHIdealGas::~PairSPHIdealGas() {
+PairSPHIdealGas::~PairSPHIdealGas()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -45,7 +51,8 @@ PairSPHIdealGas::~PairSPHIdealGas() {
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHIdealGas::compute(int eflag, int vflag) {
+void PairSPHIdealGas::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, inum, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
 
@@ -160,10 +167,6 @@ void PairSPHIdealGas::compute(int eflag, int vflag) {
 
         if (evflag)
           ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz);
-
-        if (evflag)
-          ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely,
-              delz);
       }
     }
   }
@@ -175,7 +178,8 @@ void PairSPHIdealGas::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHIdealGas::allocate() {
+void PairSPHIdealGas::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -194,7 +198,8 @@ void PairSPHIdealGas::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHIdealGas::settings(int narg, char **/*arg*/) {
+void PairSPHIdealGas::settings(int narg, char **/*arg*/)
+{
   if (narg != 0)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/idealgas");
@@ -204,7 +209,8 @@ void PairSPHIdealGas::settings(int narg, char **/*arg*/) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHIdealGas::coeff(int narg, char **arg) {
+void PairSPHIdealGas::coeff(int narg, char **arg)
+{
   if (narg != 4)
     error->all(FLERR,"Incorrect number of args for pair_style sph/idealgas coefficients");
   if (!allocated)
@@ -221,7 +227,6 @@ void PairSPHIdealGas::coeff(int narg, char **arg) {
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       viscosity[i][j] = viscosity_one;
-      //printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -236,8 +241,8 @@ void PairSPHIdealGas::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHIdealGas::init_one(int i, int j) {
-
+double PairSPHIdealGas::init_one(int i, int j)
+{
   if (setflag[i][j] == 0) {
       error->all(FLERR,"All pair sph/idealgas coeffs are not set");
   }
@@ -246,12 +251,3 @@ double PairSPHIdealGas::init_one(int i, int j) {
 
   return cut[i][j];
 }
-
-/* ---------------------------------------------------------------------- */
-
-double PairSPHIdealGas::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/,
-    double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) {
-  fforce = 0.0;
-
-  return 0.0;
-}
diff --git a/src/SPH/pair_sph_idealgas.h b/src/SPH/pair_sph_idealgas.h
index 3d61d766165..c3c11ac5bdf 100644
--- a/src/SPH/pair_sph_idealgas.h
+++ b/src/SPH/pair_sph_idealgas.h
@@ -32,7 +32,6 @@ class PairSPHIdealGas : public Pair {
   void settings(int, char **) override;
   void coeff(int, char **) override;
   double init_one(int, int) override;
-  double single(int, int, int, int, double, double, double, double &) override;
 
  protected:
   double **cut, **viscosity;
diff --git a/src/SPH/pair_sph_lj.cpp b/src/SPH/pair_sph_lj.cpp
index d59a3ad9926..4bdefca1a6d 100644
--- a/src/SPH/pair_sph_lj.cpp
+++ b/src/SPH/pair_sph_lj.cpp
@@ -13,14 +13,15 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_lj.h"
-#include <cmath>
+
 #include "atom.h"
+#include "domain.h"
+#include "error.h"
 #include "force.h"
-#include "neigh_list.h"
 #include "memory.h"
-#include "error.h"
-#include "domain.h"
+#include "neigh_list.h"
 
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -28,12 +29,17 @@ using namespace LAMMPS_NS;
 
 PairSPHLJ::PairSPHLJ(LAMMPS *lmp) : Pair(lmp)
 {
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->cv_flag != 1) || (atom->vest_flag != 1))
+    error->all(FLERR, "Pair sph/lj requires atom attributes energy, density, specific heat, and velocity estimates, e.g. in atom_style sph");
+
   restartinfo = 0;
+  single_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHLJ::~PairSPHLJ() {
+PairSPHLJ::~PairSPHLJ()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -45,7 +51,8 @@ PairSPHLJ::~PairSPHLJ() {
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHLJ::compute(int eflag, int vflag) {
+void PairSPHLJ::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, inum, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
 
@@ -182,7 +189,8 @@ void PairSPHLJ::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHLJ::allocate() {
+void PairSPHLJ::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -201,7 +209,8 @@ void PairSPHLJ::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHLJ::settings(int narg, char **/*arg*/) {
+void PairSPHLJ::settings(int narg, char **/*arg*/)
+{
   if (narg != 0)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/lj");
@@ -211,7 +220,8 @@ void PairSPHLJ::settings(int narg, char **/*arg*/) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHLJ::coeff(int narg, char **arg) {
+void PairSPHLJ::coeff(int narg, char **arg)
+{
   if (narg != 4)
     error->all(FLERR,
         "Incorrect args for pair_style sph/lj coefficients");
@@ -229,7 +239,6 @@ void PairSPHLJ::coeff(int narg, char **arg) {
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       viscosity[i][j] = viscosity_one;
-      printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -244,8 +253,8 @@ void PairSPHLJ::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHLJ::init_one(int i, int j) {
-
+double PairSPHLJ::init_one(int i, int j)
+{
   if (setflag[i][j] == 0) {
     error->all(FLERR,"All pair sph/lj coeffs are not set");
   }
@@ -256,16 +265,6 @@ double PairSPHLJ::init_one(int i, int j) {
   return cut[i][j];
 }
 
-/* ---------------------------------------------------------------------- */
-
-double PairSPHLJ::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/,
-    double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) {
-  fforce = 0.0;
-
-  return 0.0;
-}
-
-
 /*double PairSPHLJ::LJEOS2(double rho, double e, double cv) {
 
 
@@ -297,7 +296,8 @@ double PairSPHLJ::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/,
    Journal of Chemical Physics 73 pp. 5401-5403 (1980)
 */
 
-void PairSPHLJ::LJEOS2(double rho, double e, double cv, double *p, double *c) {
+void PairSPHLJ::LJEOS2(double rho, double e, double cv, double *p, double *c)
+{
   double T = e/cv;
   double beta = 1.0 / T;
   double beta_sqrt = sqrt(beta);
diff --git a/src/SPH/pair_sph_lj.h b/src/SPH/pair_sph_lj.h
index 47c34159e1d..eb93a14de9f 100644
--- a/src/SPH/pair_sph_lj.h
+++ b/src/SPH/pair_sph_lj.h
@@ -32,8 +32,6 @@ class PairSPHLJ : public Pair {
   void settings(int, char **) override;
   void coeff(int, char **) override;
   double init_one(int, int) override;
-  double single(int, int, int, int, double, double, double, double &) override;
-  //double LJEOS(int);
   void LJEOS2(double, double, double, double *, double *);
 
  protected:
diff --git a/src/SPH/pair_sph_rhosum.cpp b/src/SPH/pair_sph_rhosum.cpp
index 34adf3ce878..4fd96319eed 100644
--- a/src/SPH/pair_sph_rhosum.cpp
+++ b/src/SPH/pair_sph_rhosum.cpp
@@ -29,6 +29,9 @@ using namespace LAMMPS_NS;
 
 PairSPHRhoSum::PairSPHRhoSum(LAMMPS *lmp) : Pair(lmp)
 {
+  if (atom->rho_flag != 1)
+    error->all(FLERR, "Pair sph/rhosum requires atom attribute density, e.g. in atom_style sph");
+
   restartinfo = 0;
 
   // set comm size needed by this Pair
@@ -39,11 +42,11 @@ PairSPHRhoSum::PairSPHRhoSum(LAMMPS *lmp) : Pair(lmp)
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHRhoSum::~PairSPHRhoSum() {
+PairSPHRhoSum::~PairSPHRhoSum()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
-
     memory->destroy(cut);
   }
 }
@@ -52,14 +55,16 @@ PairSPHRhoSum::~PairSPHRhoSum() {
  init specific to this pair style
  ------------------------------------------------------------------------- */
 
-void PairSPHRhoSum::init_style() {
+void PairSPHRhoSum::init_style()
+{
   // need a full neighbor list
   neighbor->add_request(this, NeighConst::REQ_FULL);
 }
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHRhoSum::compute(int eflag, int vflag) {
+void PairSPHRhoSum::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz;
   double rsq, imass, h, ih, ihsq;
@@ -75,25 +80,6 @@ void PairSPHRhoSum::compute(int eflag, int vflag) {
   int *type = atom->type;
   double *mass = atom->mass;
 
-  // check consistency of pair coefficients
-
-  if (first) {
-    for (i = 1; i <= atom->ntypes; i++) {
-      for (j = 1; i <= atom->ntypes; i++) {
-        if (cutsq[i][j] > 0.0) {
-          if (!setflag[i][i] || !setflag[j][j]) {
-            if (comm->me == 0) {
-              printf(
-                  "SPH particle types %d and %d interact, but not all of their single particle properties are set.\n",
-                  i, j);
-            }
-          }
-        }
-      }
-    }
-    first = 0;
-  }
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -186,7 +172,6 @@ void PairSPHRhoSum::compute(int eflag, int vflag) {
 
             rho[i] += mass[jtype] * wf;
           }
-
         }
       }
     }
@@ -200,7 +185,8 @@ void PairSPHRhoSum::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHRhoSum::allocate() {
+void PairSPHRhoSum::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -210,7 +196,6 @@ void PairSPHRhoSum::allocate() {
       setflag[i][j] = 0;
 
   memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
-
   memory->create(cut, n + 1, n + 1, "pair:cut");
 }
 
@@ -218,7 +203,8 @@ void PairSPHRhoSum::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHRhoSum::settings(int narg, char **arg) {
+void PairSPHRhoSum::settings(int narg, char **arg)
+{
   if (narg != 1)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/rhosum");
@@ -229,7 +215,8 @@ void PairSPHRhoSum::settings(int narg, char **arg) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHRhoSum::coeff(int narg, char **arg) {
+void PairSPHRhoSum::coeff(int narg, char **arg)
+{
   if (narg != 3)
     error->all(FLERR,"Incorrect number of args for sph/rhosum coefficients");
   if (!allocated)
@@ -244,7 +231,6 @@ void PairSPHRhoSum::coeff(int narg, char **arg) {
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
-      //printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -259,7 +245,8 @@ void PairSPHRhoSum::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHRhoSum::init_one(int i, int j) {
+double PairSPHRhoSum::init_one(int i, int j)
+{
   if (setflag[i][j] == 0) {
     error->all(FLERR,"All pair sph/rhosum coeffs are not set");
   }
@@ -272,7 +259,8 @@ double PairSPHRhoSum::init_one(int i, int j) {
 /* ---------------------------------------------------------------------- */
 
 double PairSPHRhoSum::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/,
-    double /*factor_coul*/, double /*factor_lj*/, double &fforce) {
+    double /*factor_coul*/, double /*factor_lj*/, double &fforce)
+{
   fforce = 0.0;
 
   return 0.0;
@@ -281,7 +269,8 @@ double PairSPHRhoSum::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/,
 /* ---------------------------------------------------------------------- */
 
 int PairSPHRhoSum::pack_forward_comm(int n, int *list, double *buf,
-                                     int /*pbc_flag*/, int * /*pbc*/) {
+                                     int /*pbc_flag*/, int * /*pbc*/)
+{
   int i, j, m;
   double *rho = atom->rho;
 
@@ -295,7 +284,8 @@ int PairSPHRhoSum::pack_forward_comm(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHRhoSum::unpack_forward_comm(int n, int first, double *buf) {
+void PairSPHRhoSum::unpack_forward_comm(int n, int first, double *buf)
+{
   int i, m, last;
   double *rho = atom->rho;
 
diff --git a/src/SPH/pair_sph_taitwater.cpp b/src/SPH/pair_sph_taitwater.cpp
index 9a5991718de..d97492de638 100644
--- a/src/SPH/pair_sph_taitwater.cpp
+++ b/src/SPH/pair_sph_taitwater.cpp
@@ -13,15 +13,16 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_taitwater.h"
-#include <cmath>
+
 #include "atom.h"
-#include "force.h"
 #include "comm.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
 #include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "neigh_list.h"
 
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -29,6 +30,9 @@ using namespace LAMMPS_NS;
 
 PairSPHTaitwater::PairSPHTaitwater(LAMMPS *lmp) : Pair(lmp)
 {
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1))
+    error->all(FLERR, "Pair sph/taitwater requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph");
+
   restartinfo = 0;
   single_enable = 0;
   first = 1;
@@ -36,7 +40,8 @@ PairSPHTaitwater::PairSPHTaitwater(LAMMPS *lmp) : Pair(lmp)
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHTaitwater::~PairSPHTaitwater() {
+PairSPHTaitwater::~PairSPHTaitwater()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -51,7 +56,8 @@ PairSPHTaitwater::~PairSPHTaitwater() {
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHTaitwater::compute(int eflag, int vflag) {
+void PairSPHTaitwater::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, inum, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
 
@@ -72,25 +78,6 @@ void PairSPHTaitwater::compute(int eflag, int vflag) {
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
 
-  // check consistency of pair coefficients
-
-  if (first) {
-    for (i = 1; i <= atom->ntypes; i++) {
-      for (j = 1; i <= atom->ntypes; i++) {
-        if (cutsq[i][j] > 1.e-32) {
-          if (!setflag[i][i] || !setflag[j][j]) {
-            if (comm->me == 0) {
-              printf(
-                  "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n",
-                  i, j, sqrt(cutsq[i][j]));
-            }
-          }
-        }
-      }
-    }
-    first = 0;
-  }
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -201,7 +188,8 @@ void PairSPHTaitwater::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwater::allocate() {
+void PairSPHTaitwater::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -223,7 +211,8 @@ void PairSPHTaitwater::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwater::settings(int narg, char **/*arg*/) {
+void PairSPHTaitwater::settings(int narg, char **/*arg*/)
+{
   if (narg != 0)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/taitwater");
@@ -233,7 +222,8 @@ void PairSPHTaitwater::settings(int narg, char **/*arg*/) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwater::coeff(int narg, char **arg) {
+void PairSPHTaitwater::coeff(int narg, char **arg)
+{
   if (narg != 6)
     error->all(FLERR,
         "Incorrect args for pair_style sph/taitwater coefficients");
@@ -257,14 +247,8 @@ void PairSPHTaitwater::coeff(int narg, char **arg) {
     B[i] = B_one;
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       viscosity[i][j] = viscosity_one;
-      //printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
-
       setflag[i][j] = 1;
-
-      //cut[j][i] = cut[i][j];
-      //viscosity[j][i] = viscosity[i][j];
-      //setflag[j][i] = 1;
       count++;
     }
   }
@@ -277,8 +261,8 @@ void PairSPHTaitwater::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHTaitwater::init_one(int i, int j) {
-
+double PairSPHTaitwater::init_one(int i, int j)
+{
   if (setflag[i][j] == 0) {
     error->all(FLERR,"All pair sph/taitwater coeffs are set");
   }
diff --git a/src/SPH/pair_sph_taitwater_morris.cpp b/src/SPH/pair_sph_taitwater_morris.cpp
index 50fcd270f6c..2209cc4ecf2 100644
--- a/src/SPH/pair_sph_taitwater_morris.cpp
+++ b/src/SPH/pair_sph_taitwater_morris.cpp
@@ -13,15 +13,16 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_taitwater_morris.h"
-#include <cmath>
+
 #include "atom.h"
-#include "force.h"
 #include "comm.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
 #include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "neigh_list.h"
 
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -29,6 +30,9 @@ using namespace LAMMPS_NS;
 
 PairSPHTaitwaterMorris::PairSPHTaitwaterMorris(LAMMPS *lmp) : Pair(lmp)
 {
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1))
+    error->all(FLERR, "Pair sph/taitwater/morris requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph");
+
   restartinfo = 0;
   first = 1;
   single_enable = 0;
@@ -36,7 +40,8 @@ PairSPHTaitwaterMorris::PairSPHTaitwaterMorris(LAMMPS *lmp) : Pair(lmp)
 
 /* ---------------------------------------------------------------------- */
 
-PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() {
+PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris()
+{
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -51,7 +56,8 @@ PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() {
 
 /* ---------------------------------------------------------------------- */
 
-void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
+void PairSPHTaitwaterMorris::compute(int eflag, int vflag)
+{
   int i, j, ii, jj, inum, jnum, itype, jtype;
   double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
 
@@ -72,25 +78,6 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
 
-  // check consistency of pair coefficients
-
-  if (first) {
-    for (i = 1; i <= atom->ntypes; i++) {
-      for (j = 1; i <= atom->ntypes; i++) {
-        if (cutsq[i][j] > 1.e-32) {
-          if (!setflag[i][i] || !setflag[j][j]) {
-            if (comm->me == 0) {
-              printf(
-                  "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n",
-                  i, j, sqrt(cutsq[i][j]));
-            }
-          }
-        }
-      }
-    }
-    first = 0;
-  }
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -152,9 +139,9 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
         fj = tmp * tmp * tmp;
         fj = B[jtype] * (fj * fj * tmp - 1.0) / (rho[j] * rho[j]);
 
-        velx=vxtmp - v[j][0];
-        vely=vytmp - v[j][1];
-        velz=vztmp - v[j][2];
+        velx = vxtmp - v[j][0];
+        vely = vytmp - v[j][1];
+        velz = vztmp - v[j][2];
 
         // dot product of velocity delta and distance vector
         delVdotDelR = delx * velx + dely * vely + delz * velz;
@@ -169,8 +156,6 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
         fpair = -imass * jmass * (fi + fj) * wfd;
         deltaE = -0.5 *(fpair * delVdotDelR + fvisc * (velx*velx + vely*vely + velz*velz));
 
-       // printf("testvar= %f, %f \n", delx, dely);
-
         f[i][0] += delx * fpair + velx * fvisc;
         f[i][1] += dely * fpair + vely * fvisc;
         f[i][2] += delz * fpair + velz * fvisc;
@@ -189,6 +174,7 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
           drho[j] += imass * delVdotDelR * wfd;
         }
 
+        // viscous forces do not contribute to virial
         if (evflag)
           ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz);
       }
@@ -202,7 +188,8 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) {
  allocate all arrays
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwaterMorris::allocate() {
+void PairSPHTaitwaterMorris::allocate()
+{
   allocated = 1;
   int n = atom->ntypes;
 
@@ -224,7 +211,8 @@ void PairSPHTaitwaterMorris::allocate() {
  global settings
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/) {
+void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/)
+{
   if (narg != 0)
     error->all(FLERR,
         "Illegal number of arguments for pair_style sph/taitwater/morris");
@@ -234,7 +222,8 @@ void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/) {
  set coeffs for one or more type pairs
  ------------------------------------------------------------------------- */
 
-void PairSPHTaitwaterMorris::coeff(int narg, char **arg) {
+void PairSPHTaitwaterMorris::coeff(int narg, char **arg)
+{
   if (narg != 6)
     error->all(FLERR,
         "Incorrect args for pair_style sph/taitwater/morris coefficients");
@@ -258,7 +247,6 @@ void PairSPHTaitwaterMorris::coeff(int narg, char **arg) {
     B[i] = B_one;
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       viscosity[i][j] = viscosity_one;
-      //printf("setting cut[%d][%d] = %f\n", i, j, cut_one);
       cut[i][j] = cut_one;
 
       setflag[i][j] = 1;
@@ -274,8 +262,8 @@ void PairSPHTaitwaterMorris::coeff(int narg, char **arg) {
  init for one type pair i,j and corresponding j,i
  ------------------------------------------------------------------------- */
 
-double PairSPHTaitwaterMorris::init_one(int i, int j) {
-
+double PairSPHTaitwaterMorris::init_one(int i, int j)
+{
   if (setflag[i][j] == 0) {
     error->all(FLERR,"All pair sph/taitwater/morris coeffs are not set");
   }

From 6abbbdba6c98f9ba77b56dc5820053dc6ec01310 Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 17:10:23 -0600
Subject: [PATCH 2/7] Patching vremap in SPH

---
 src/SPH/fix_sph.cpp | 42 +++++++++++++++++++++++++++++++++++++++---
 src/SPH/fix_sph.h   |  2 ++
 2 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/src/SPH/fix_sph.cpp b/src/SPH/fix_sph.cpp
index 96479be9281..876d7ac4fb5 100644
--- a/src/SPH/fix_sph.cpp
+++ b/src/SPH/fix_sph.cpp
@@ -61,6 +61,10 @@ void FixSPH::init()
 
 void FixSPH::setup_pre_force(int /*vflag*/)
 {
+  remap_v_flag = domain->deform_vremap;
+  if (remap_v_flag && (!comm->ghost_velocity))
+    error->all(FLERR, "Fix sph requires ghost atoms store velocity when deforming with remap v");
+
   // set vest equal to v
   double **v = atom->v;
   double **vest = atom->vest;
@@ -119,9 +123,16 @@ void FixSPH::initial_integrate(int /*vflag*/)
       rho[i] += dtf * drho[i]; // ... and density
 
       // extrapolate velocity for use with velocity-dependent potentials, e.g. SPH
-      vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0];
-      vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1];
-      vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2];
+      // if velocities are remapped, perform this extrapolation after communication
+      if (remap_v_flag) {
+        vest[i][0] = dtfm * f[i][0];
+        vest[i][1] = dtfm * f[i][1];
+        vest[i][2] = dtfm * f[i][2];
+      } else {
+        vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0];
+        vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1];
+        vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2];
+      }
 
       v[i][0] += dtfm * f[i][0];
       v[i][1] += dtfm * f[i][1];
@@ -136,6 +147,31 @@ void FixSPH::initial_integrate(int /*vflag*/)
 
 /* ---------------------------------------------------------------------- */
 
+void FixSPH::pre_force(int /*vflag*/)
+{
+  // if velocities are remapped, calculate estimates here
+  // note that vest currently stores dtfm * force
+  if (!remap_v_flag) return;
+
+  double **v = atom->v;
+  double **vest = atom->vest;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup)
+    nlocal = atom->nfirst;
+
+  int nall = nlocal + atom->nghost;
+  for (int i = 0; i < nall; i++) {
+    if (mask[i] & groupbit) {
+      vest[i][0] += v[i][0];
+      vest[i][1] += v[i][1];
+      vest[i][2] += v[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixSPH::final_integrate()
 {
   // update v, rho, and e of atoms in group
diff --git a/src/SPH/fix_sph.h b/src/SPH/fix_sph.h
index 415e38bd214..7676844dd94 100644
--- a/src/SPH/fix_sph.h
+++ b/src/SPH/fix_sph.h
@@ -30,6 +30,7 @@ class FixSPH : public Fix {
   int setmask() override;
   void init() override;
   void setup_pre_force(int) override;
+  void pre_force(int) override;
   void initial_integrate(int) override;
   void final_integrate() override;
   void reset_dt() override;
@@ -41,6 +42,7 @@ class FixSPH : public Fix {
   double dtv, dtf;
   double *step_respa;
   int mass_require;
+  int remap_v_flag;
 
   class Pair *pair;
 };

From 49b377fc3db56bac056cff1009c55c896c23f644 Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 17:51:48 -0600
Subject: [PATCH 3/7] Adding errors to unpatched uses of vest with vremap

---
 src/DPD-MESO/fix_mvv_dpd.cpp            | 5 +++++
 src/DPD-MESO/fix_mvv_edpd.cpp           | 5 +++++
 src/DPD-MESO/fix_mvv_tdpd.cpp           | 5 +++++
 src/DPD-SMOOTH/fix_meso_move.cpp        | 7 +++++++
 src/DPD-SMOOTH/fix_rigid_meso.cpp       | 5 +++++
 src/MACHDYN/fix_smd_integrate_tlsph.cpp | 5 +++++
 src/MACHDYN/fix_smd_integrate_ulsph.cpp | 5 +++++
 7 files changed, 37 insertions(+)

diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp
index d51000b15b4..6def3305d8b 100644
--- a/src/DPD-MESO/fix_mvv_dpd.cpp
+++ b/src/DPD-MESO/fix_mvv_dpd.cpp
@@ -65,6 +65,11 @@ void FixMvvDPD::init()
   if (!atom->vest_flag)
     error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/dpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^mdpd",0) && !force->pair_match("^dpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!(force->pair_match("^mdpd",0,1) || force->pair_match("^dpd",0),1)) {
diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp
index 7ac0dc3de7d..389a8c97bf1 100644
--- a/src/DPD-MESO/fix_mvv_edpd.cpp
+++ b/src/DPD-MESO/fix_mvv_edpd.cpp
@@ -73,6 +73,11 @@ void FixMvvEDPD::init()
 {
   if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/edpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^edpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!force->pair_match("^edpd",0,1)) {
diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp
index 122fd563652..1babbb8c5fc 100644
--- a/src/DPD-MESO/fix_mvv_tdpd.cpp
+++ b/src/DPD-MESO/fix_mvv_tdpd.cpp
@@ -71,6 +71,11 @@ void FixMvvTDPD::init()
 {
   if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/tdpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^tdpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!force->pair_match("^tdpd",0,1)) {
diff --git a/src/DPD-SMOOTH/fix_meso_move.cpp b/src/DPD-SMOOTH/fix_meso_move.cpp
index 2c3213c3cd2..89a7814f6a5 100644
--- a/src/DPD-SMOOTH/fix_meso_move.cpp
+++ b/src/DPD-SMOOTH/fix_meso_move.cpp
@@ -350,7 +350,14 @@ void FixMesoMove::init () {
 }
 
 void FixMesoMove::setup_pre_force (int /*vflag*/) {
+
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix meso/move cannot be used with velocity remapping");
+
   // set vest equal to v
+
   double **v = atom->v;
   double **vest = atom->vest;
   int *mask = atom->mask;
diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp
index e0ad501e023..bfceb3295e8 100644
--- a/src/DPD-SMOOTH/fix_rigid_meso.cpp
+++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp
@@ -92,6 +92,11 @@ void FixRigidMeso::setup (int vflag) {
     conjqm[ibody][2] *= 2.0;
     conjqm[ibody][3] *= 2.0;
   }
+
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix rigid/meso cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
index f138a3c387b..97cce5524db 100644
--- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
@@ -115,6 +115,11 @@ void FixSMDIntegrateTlsph::init() {
         dtv = update->dt;
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
+
+        // Cannot use vremap since its effects aren't       propagated to vest
+        //   see RHEO or SPH packages for examples patches
+        if (domain->deform_vremap)
+          error->all(FLERR, "Fix smd/integrate_tlsph cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
index 8ef5c65b676..395946ba208 100644
--- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
@@ -144,6 +144,11 @@ void FixSMDIntegrateUlsph::init() {
         dtv = update->dt;
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
+
+        // Cannot use vremap since its effects aren't       propagated to vest
+        //   see RHEO or SPH packages for examples patches
+        if (domain->deform_vremap)
+          error->all(FLERR, "Fix smd/integrate_ulsph cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------

From d43e87bce1acfb33b5b4bb913c5bfc89a0039b1a Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 17:54:52 -0600
Subject: [PATCH 4/7] Missing word

---
 src/DPD-MESO/fix_mvv_dpd.cpp            | 2 +-
 src/DPD-MESO/fix_mvv_edpd.cpp           | 2 +-
 src/DPD-MESO/fix_mvv_tdpd.cpp           | 2 +-
 src/DPD-SMOOTH/fix_meso_move.cpp        | 2 +-
 src/DPD-SMOOTH/fix_rigid_meso.cpp       | 2 +-
 src/MACHDYN/fix_smd_integrate_tlsph.cpp | 4 ++--
 src/MACHDYN/fix_smd_integrate_ulsph.cpp | 4 ++--
 7 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp
index 6def3305d8b..5e91b4bad91 100644
--- a/src/DPD-MESO/fix_mvv_dpd.cpp
+++ b/src/DPD-MESO/fix_mvv_dpd.cpp
@@ -66,7 +66,7 @@ void FixMvvDPD::init()
     error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd");
 
   // Cannot use vremap since its effects aren't propagated to vest
-  //   see RHEO or SPH packages for examples patches
+  //   see RHEO or SPH packages for examples of patches
   if (domain->deform_vremap)
     error->all(FLERR, "Fix mvv/dpd cannot be used with velocity remapping");
 
diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp
index 389a8c97bf1..d34e9b5e5ab 100644
--- a/src/DPD-MESO/fix_mvv_edpd.cpp
+++ b/src/DPD-MESO/fix_mvv_edpd.cpp
@@ -74,7 +74,7 @@ void FixMvvEDPD::init()
   if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd");
 
   // Cannot use vremap since its effects aren't propagated to vest
-  //   see RHEO or SPH packages for examples patches
+  //   see RHEO or SPH packages for examples of patches
   if (domain->deform_vremap)
     error->all(FLERR, "Fix mvv/edpd cannot be used with velocity remapping");
 
diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp
index 1babbb8c5fc..a34c86b6dda 100644
--- a/src/DPD-MESO/fix_mvv_tdpd.cpp
+++ b/src/DPD-MESO/fix_mvv_tdpd.cpp
@@ -72,7 +72,7 @@ void FixMvvTDPD::init()
   if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd");
 
   // Cannot use vremap since its effects aren't propagated to vest
-  //   see RHEO or SPH packages for examples patches
+  //   see RHEO or SPH packages for examples of patches
   if (domain->deform_vremap)
     error->all(FLERR, "Fix mvv/tdpd cannot be used with velocity remapping");
 
diff --git a/src/DPD-SMOOTH/fix_meso_move.cpp b/src/DPD-SMOOTH/fix_meso_move.cpp
index 89a7814f6a5..c9d4a414009 100644
--- a/src/DPD-SMOOTH/fix_meso_move.cpp
+++ b/src/DPD-SMOOTH/fix_meso_move.cpp
@@ -352,7 +352,7 @@ void FixMesoMove::init () {
 void FixMesoMove::setup_pre_force (int /*vflag*/) {
 
   // Cannot use vremap since its effects aren't propagated to vest
-  //   see RHEO or SPH packages for examples patches
+  //   see RHEO or SPH packages for examples of patches
   if (domain->deform_vremap)
     error->all(FLERR, "Fix meso/move cannot be used with velocity remapping");
 
diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp
index bfceb3295e8..54fd72f9247 100644
--- a/src/DPD-SMOOTH/fix_rigid_meso.cpp
+++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp
@@ -94,7 +94,7 @@ void FixRigidMeso::setup (int vflag) {
   }
 
   // Cannot use vremap since its effects aren't propagated to vest
-  //   see RHEO or SPH packages for examples patches
+  //   see RHEO or SPH packages for examples of patches
   if (domain->deform_vremap)
     error->all(FLERR, "Fix rigid/meso cannot be used with velocity remapping");
 }
diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
index 97cce5524db..e00c88bc225 100644
--- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
@@ -116,8 +116,8 @@ void FixSMDIntegrateTlsph::init() {
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
 
-        // Cannot use vremap since its effects aren't       propagated to vest
-        //   see RHEO or SPH packages for examples patches
+        // Cannot use vremap since its effects aren't propagated to vest
+        //   see RHEO or SPH packages for examples of patches
         if (domain->deform_vremap)
           error->all(FLERR, "Fix smd/integrate_tlsph cannot be used with velocity remapping");
 }
diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
index 395946ba208..c6530c989bc 100644
--- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
@@ -145,8 +145,8 @@ void FixSMDIntegrateUlsph::init() {
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
 
-        // Cannot use vremap since its effects aren't       propagated to vest
-        //   see RHEO or SPH packages for examples patches
+        // Cannot use vremap since its effects aren't propagated to vest
+        //   see RHEO or SPH packages for examples of patches
         if (domain->deform_vremap)
           error->all(FLERR, "Fix smd/integrate_ulsph cannot be used with velocity remapping");
 }

From 3deffb0dfd96fb7b57d735762da56bc6d71694ad Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 18:10:27 -0600
Subject: [PATCH 5/7] typo

---
 src/SPH/compute_sph_e_atom.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/SPH/compute_sph_e_atom.cpp b/src/SPH/compute_sph_e_atom.cpp
index 8bb0622acd2..b98fd974795 100644
--- a/src/SPH/compute_sph_e_atom.cpp
+++ b/src/SPH/compute_sph_e_atom.cpp
@@ -33,7 +33,7 @@ ComputeSPHEAtom::ComputeSPHEAtom(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3)
     error->all(FLERR,"Number of arguments for compute sph/e/atom command != 3");
   if (atom->esph_flag != 1)
-    error->all(FLERR,"Compute sph/e/atom requires atom attribut energy, e.g. in atom_style sph");
+    error->all(FLERR,"Compute sph/e/atom requires atom attribute energy, e.g. in atom_style sph");
 
   peratom_flag = 1;
   size_peratom_cols = 0;

From 68a6bc069330a2e6583af78e67f93be70e1c4bcf Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 18:15:34 -0600
Subject: [PATCH 6/7] Adding missing header file

---
 src/DPD-MESO/fix_mvv_dpd.cpp            |  1 +
 src/DPD-MESO/fix_mvv_edpd.cpp           |  1 +
 src/DPD-MESO/fix_mvv_tdpd.cpp           |  1 +
 src/DPD-SMOOTH/fix_rigid_meso.cpp       |  3 ++-
 src/MACHDYN/fix_smd_integrate_tlsph.cpp | 15 +++++++++------
 src/MACHDYN/fix_smd_integrate_ulsph.cpp | 13 ++++++++-----
 6 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp
index 5e91b4bad91..6b5440902bb 100644
--- a/src/DPD-MESO/fix_mvv_dpd.cpp
+++ b/src/DPD-MESO/fix_mvv_dpd.cpp
@@ -24,6 +24,7 @@
 #include "fix_mvv_dpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp
index d34e9b5e5ab..32a1608d6e9 100644
--- a/src/DPD-MESO/fix_mvv_edpd.cpp
+++ b/src/DPD-MESO/fix_mvv_edpd.cpp
@@ -33,6 +33,7 @@
 #include "fix_mvv_edpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp
index a34c86b6dda..d4e2d8a5c9c 100644
--- a/src/DPD-MESO/fix_mvv_tdpd.cpp
+++ b/src/DPD-MESO/fix_mvv_tdpd.cpp
@@ -29,6 +29,7 @@
 #include "fix_mvv_tdpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp
index 54fd72f9247..3e8b07e6dcf 100644
--- a/src/DPD-SMOOTH/fix_rigid_meso.cpp
+++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp
@@ -29,11 +29,12 @@
 ------------------------------------------------------------------------- */
 
 #include "fix_rigid_meso.h"
+
 #include "math_extra.h"
 #include "atom.h"
 #include "domain.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
index e00c88bc225..ee617c15775 100644
--- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp
@@ -24,15 +24,18 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_smd_integrate_tlsph.h"
-#include <cmath>
-#include <cstring>
-#include <Eigen/Eigen>
+
 #include "atom.h"
-#include "force.h"
-#include "update.h"
+#include "comm.h"
+#include "domain.h"
 #include "error.h"
+#include "force.h"
 #include "pair.h"
-#include "comm.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+#include <Eigen/Eigen>
 
 using namespace Eigen;
 using namespace LAMMPS_NS;
diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
index c6530c989bc..6e2934c2ec0 100644
--- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp
+++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp
@@ -24,15 +24,18 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_smd_integrate_ulsph.h"
-#include <cmath>
-#include <cstring>
-#include <Eigen/Eigen>
+
 #include "atom.h"
 #include "comm.h"
-#include "force.h"
-#include "update.h"
+#include "domain.h"
 #include "error.h"
+#include "force.h"
 #include "pair.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+#include <Eigen/Eigen>
 
 using namespace Eigen;
 using namespace LAMMPS_NS;

From a972a282a792e4ad6cfd23e746cf7289b4c59b80 Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Fri, 26 Jul 2024 17:43:04 -0600
Subject: [PATCH 7/7] Adding remap v restriction to doc files

---
 doc/src/fix_meso_move.rst           | 5 +++++
 doc/src/fix_mvv_dpd.rst             | 5 +++++
 doc/src/fix_rigid_meso.rst          | 5 +++++
 doc/src/fix_smd_integrate_tlsph.rst | 5 +++++
 doc/src/fix_smd_integrate_ulsph.rst | 5 +++++
 src/DPD-SMOOTH/fix_rigid_meso.cpp   | 2 +-
 6 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/doc/src/fix_meso_move.rst b/doc/src/fix_meso_move.rst
index 55d54b21077..281a405ab06 100644
--- a/doc/src/fix_meso_move.rst
+++ b/doc/src/fix_meso_move.rst
@@ -247,6 +247,11 @@ defined by the :doc:`atom_style sph <atom_style>` command.
 
 All particles in the group must be mesoscopic SPH/SDPD particles.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 
diff --git a/doc/src/fix_mvv_dpd.rst b/doc/src/fix_mvv_dpd.rst
index ff5b169f973..740785d5628 100644
--- a/doc/src/fix_mvv_dpd.rst
+++ b/doc/src/fix_mvv_dpd.rst
@@ -97,6 +97,11 @@ These fixes are part of the DPD-MESO package. They are only enabled if
 LAMMPS was built with that package. See the :doc:`Build package
 <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 
diff --git a/doc/src/fix_rigid_meso.rst b/doc/src/fix_rigid_meso.rst
index c8c994fd26c..933369788ed 100644
--- a/doc/src/fix_rigid_meso.rst
+++ b/doc/src/fix_rigid_meso.rst
@@ -353,6 +353,11 @@ defined by the :doc:`atom_style sph <atom_style>` command.
 
 All particles in the group must be mesoscopic SPH/SDPD particles.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 
diff --git a/doc/src/fix_smd_integrate_tlsph.rst b/doc/src/fix_smd_integrate_tlsph.rst
index 515b30d4ff5..e714c91a174 100644
--- a/doc/src/fix_smd_integrate_tlsph.rst
+++ b/doc/src/fix_smd_integrate_tlsph.rst
@@ -53,6 +53,11 @@ Restrictions
 This fix is part of the MACHDYN package.  It is only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 
diff --git a/doc/src/fix_smd_integrate_ulsph.rst b/doc/src/fix_smd_integrate_ulsph.rst
index 17dfdb7b174..60c185db5f4 100644
--- a/doc/src/fix_smd_integrate_ulsph.rst
+++ b/doc/src/fix_smd_integrate_ulsph.rst
@@ -61,6 +61,11 @@ Restrictions
 This fix is part of the MACHDYN package.  It is only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 
diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp
index 3e8b07e6dcf..38b9d0a09ca 100644
--- a/src/DPD-SMOOTH/fix_rigid_meso.cpp
+++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp
@@ -30,10 +30,10 @@
 
 #include "fix_rigid_meso.h"
 
-#include "math_extra.h"
 #include "atom.h"
 #include "domain.h"
 #include "error.h"
+#include "math_extra.h"
 #include "memory.h"
 
 using namespace LAMMPS_NS;