Cleaning up SPH package, minor bug fixes

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");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -81,8 +82,7 @@ void ComputeSPHEAtom::compute_peratom()
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       evector[i] = esph[i];
-      }
+    } else {
-      else {
       evector[i] = 0.0;
     }
   }

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,8 +73,6 @@ 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;
@@ -82,8 +80,7 @@ void ComputeSPHRhoAtom::compute_peratom()
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       rhoVector[i] = rho[i];
-      }
+    } else {
-      else {
       rhoVector[i] = 0.0;
     }
   }

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");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -84,8 +85,7 @@ void ComputeSPHTAtom::compute_peratom()
       if (cv[i] > 0.0) {
         tvector[i] = esph[i] / cv[i];
       }
-      }
+    } else {
-      else {
       tvector[i] = 0.0;
     }
   }

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) {
+  Fix(lmp, narg, arg)
-
+{
-  if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
+  if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1))
-    error->all(FLERR,
+    error->all(FLERR, "Fix sph requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph");
-        "Fix sph command requires atom_style with both energy and density");
 
   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;
 }

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,
+    error->all(FLERR, "Fix sph/stationary requires atom attributes energy and density, e.g. in atom_style sph");
-        "Fix sph/stationary command requires atom_style with both energy and density, e.g. meso");
 
   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;
 }

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;

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;

src/SPH/pair_sph_idealgas.cpp

@@ -13,14 +13,15 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_idealgas.h"
-#include <cmath>
-#include "atom.h"
-#include "force.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
-#include "domain.h"
 
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory.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;
-}

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;

src/SPH/pair_sph_lj.cpp

@@ -13,14 +13,15 @@
  ------------------------------------------------------------------------- */
 
 #include "pair_sph_lj.h"
-#include <cmath>
-#include "atom.h"
-#include "force.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
-#include "domain.h"
 
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory.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);

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:

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;
 

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 "atom.h"
+#include "comm.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");
   }

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 "atom.h"
+#include "comm.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];
+        velx = vxtmp - v[j][0];
-        vely=vytmp - v[j][1];
+        vely = vytmp - v[j][1];
-        velz=vztmp - v[j][2];
+        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");
   }