Various small cleanups to RHEO package flagged by compiler warnings

The individual changes are:
- remove of unused function parameters
- replace non-standard variable length arrays on the stack with static ones
- disable citation removed from the manual
- replace #defined constants with enum or static constexpr
- enable and apply clang-format

src/RHEO/atom_vec_rheo.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories

src/RHEO/atom_vec_rheo_thermal.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -47,8 +46,10 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
   // order of fields in a string does not matter
   // except: fields_data_atom & fields_data_vel must match data file
 
-  fields_grow = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
+  fields_grow = {"rheo_status", "rho",          "drho",     "temperature", "esph",
-  fields_copy = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
+                 "heatflow",    "conductivity", "pressure", "viscosity"};
+  fields_copy = {"rheo_status", "rho",          "drho",     "temperature", "esph",
+                 "heatflow",    "conductivity", "pressure", "viscosity"};
   fields_comm = {"rheo_status", "rho", "esph"};
   fields_comm_vel = {"rheo_status", "rho", "esph"};
   fields_reverse = {"drho", "heatflow"};
@@ -56,7 +57,8 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
   fields_border_vel = {"rheo_status", "rho", "esph"};
   fields_exchange = {"rheo_status", "rho", "esph"};
   fields_restart = {"rheo_status", "rho", "esph"};
-  fields_create = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
+  fields_create = {"rheo_status", "rho",          "drho",     "temperature", "esph",
+                   "heatflow",    "conductivity", "pressure", "viscosity"};
   fields_data_atom = {"id", "type", "rheo_status", "rho", "esph", "x"};
   fields_data_vel = {"id", "v"};
 

src/RHEO/bond_rheo_shell.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -37,7 +36,7 @@
 #include <cmath>
 #include <cstring>
 
-#define EPSILON 1e-10
+static constexpr double EPSILON = 1e-10;
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
@@ -45,8 +44,8 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
-  BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr),
+    BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr),
-  id_fix(nullptr), compute_surface(nullptr)
+    id_fix(nullptr), compute_surface(nullptr)
 {
   partial_flag = 1;
   comm_reverse = 1;
@@ -168,8 +167,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
     store_data();
   }
 
-  if (hybrid_flag)
+  if (hybrid_flag) fix_bond_history->compress_history();
-    fix_bond_history->compress_history();
 
   int i1, i2, itmp, n, type;
   double delx, dely, delz, delvx, delvy, delvz;
@@ -191,7 +189,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
 
   double **bondstore = fix_bond_history->bondstore;
 
-  if (atom->nmax > nmax_store){
+  if (atom->nmax > nmax_store) {
     nmax_store = atom->nmax;
     memory->destroy(dbond);
     memory->create(dbond, nmax_store, "rheo/shell:dbond");
@@ -226,8 +224,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
     r = sqrt(rsq);
 
     // If bond hasn't been set - zero data
-    if (t < EPSILON || std::isnan(t))
+    if (t < EPSILON || std::isnan(t)) t = store_bond(n, i1, i2);
-      t = store_bond(n, i1, i2);
 
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
@@ -297,15 +294,14 @@ void BondRHEOShell::compute(int eflag, int vflag)
   // Communicate changes in nbond
   if (newton_bond) comm->reverse_comm(this);
 
-  for(int i = 0; i < nlocal; i++) {
+  for (int i = 0; i < nlocal; i++) {
     nbond[i] += dbond[i];
 
     // If it has bonds, no shifting
     if (nbond[i] != 0) status[i] |= STATUS_NO_SHIFT;
   }
 
-  if (hybrid_flag)
+  if (hybrid_flag) fix_bond_history->uncompress_history();
-    fix_bond_history->uncompress_history();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -368,12 +364,13 @@ void BondRHEOShell::init_style()
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell");
   class FixRHEO *fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
-  if (!fix_rheo->surface_flag) error->all(FLERR,
+  if (!fix_rheo->surface_flag)
-      "Bond rheo/shell requires surface calculation in fix rheo");
+    error->all(FLERR, "Bond rheo/shell requires surface calculation in fix rheo");
   compute_surface = fix_rheo->compute_surface;
 
   fixes = modify->get_fix_by_style("^rheo/oxidation$");
-  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell");
+  if (fixes.size() == 0)
+    error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell");
   class FixRHEOOxidation *fix_rheo_oxidation = dynamic_cast<FixRHEOOxidation *>(fixes[0]);
 
   rsurf = fix_rheo_oxidation->rsurf;
@@ -402,7 +399,6 @@ void BondRHEOShell::settings(int narg, char **arg)
     error->all(FLERR, "Illegal bond rheo/shell command, must specify positive formation time");
 }
 
-
 /* ----------------------------------------------------------------------
    used to check bond communiction cutoff - not perfect, estimates based on local-local only
 ------------------------------------------------------------------------- */
@@ -464,13 +460,10 @@ void BondRHEOShell::write_restart_settings(FILE *fp)
 
 void BondRHEOShell::read_restart_settings(FILE *fp)
 {
-  if (comm->me == 0) {
+  if (comm->me == 0) { utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error); }
-    utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error);
-  }
   MPI_Bcast(&tform, 1, MPI_DOUBLE, 0, world);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf)
@@ -479,9 +472,7 @@ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf)
   m = 0;
   last = first + n;
 
-  for (i = first; i < last; i++) {
+  for (i = first; i < last; i++) { buf[m++] = dbond[i]; }
-    buf[m++] = dbond[i];
-  }
   return m;
 }
 
@@ -561,8 +552,8 @@ void BondRHEOShell::process_ineligibility(int i, int j)
         n = num_bond[i];
         bond_type[i][m] = bond_type[i][n - 1];
         bond_atom[i][m] = bond_atom[i][n - 1];
-        for (auto &ihistory: histories) {
+        for (auto &ihistory : histories) {
-          auto fix_bond_history2 = dynamic_cast<FixBondHistory *> (ihistory);
+          auto fix_bond_history2 = dynamic_cast<FixBondHistory *>(ihistory);
           fix_bond_history2->shift_history(i, m, n - 1);
           fix_bond_history2->delete_history(i, n - 1);
         }
@@ -579,8 +570,8 @@ void BondRHEOShell::process_ineligibility(int i, int j)
         n = num_bond[j];
         bond_type[j][m] = bond_type[j][n - 1];
         bond_atom[j][m] = bond_atom[j][n - 1];
-        for (auto &ihistory: histories) {
+        for (auto &ihistory : histories) {
-          auto fix_bond_history2 = dynamic_cast<FixBondHistory *> (ihistory);
+          auto fix_bond_history2 = dynamic_cast<FixBondHistory *>(ihistory);
           fix_bond_history2->shift_history(j, m, n - 1);
           fix_bond_history2->delete_history(j, n - 1);
         }

src/RHEO/compute_rheo_grad.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -21,15 +20,15 @@
 
 #include "atom.h"
 #include "comm.h"
-#include "compute_rheo_kernel.h"
 #include "compute_rheo_interface.h"
+#include "compute_rheo_kernel.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
 #include "force.h"
 #include "memory.h"
-#include "neighbor.h"
 #include "neigh_list.h"
+#include "neighbor.h"
 #include "update.h"
 
 #include <cmath>
@@ -37,24 +36,29 @@
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
 
-enum{COMMGRAD, COMMFIELD};
+enum { COMMGRAD, COMMFIELD };
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr),
+    Compute(lmp, narg, arg), gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr),
-  fix_rheo(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
+    fix_rheo(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
-  list(nullptr)
+    list(nullptr)
 {
-  if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
+  if (narg < 4) error->all(FLERR, "Illegal compute rheo/grad command");
 
   velocity_flag = energy_flag = rho_flag = eta_flag = 0;
   for (int iarg = 3; iarg < narg; iarg++) {
-    if (strcmp(arg[iarg], "velocity") == 0) velocity_flag = 1;
+    if (strcmp(arg[iarg], "velocity") == 0)
-    else if (strcmp(arg[iarg], "rho") == 0) rho_flag = 1;
+      velocity_flag = 1;
-    else if (strcmp(arg[iarg], "energy") == 0) energy_flag = 1;
+    else if (strcmp(arg[iarg], "rho") == 0)
-    else if (strcmp(arg[iarg], "viscosity") == 0) eta_flag = 1;
+      rho_flag = 1;
-    else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]);
+    else if (strcmp(arg[iarg], "energy") == 0)
+      energy_flag = 1;
+    else if (strcmp(arg[iarg], "viscosity") == 0)
+      eta_flag = 1;
+    else
+      error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]);
   }
 
   ncomm_grad = 0;
@@ -89,7 +93,6 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
 
   nmax_store = 0;
   grow_arrays(atom->nmax);
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -161,20 +164,16 @@ void ComputeRHEOGrad::compute_peratom()
 
   for (i = 0; i < nmax_store; i++) {
     if (velocity_flag) {
-      for (k = 0; k < dim * dim; k++)
+      for (k = 0; k < dim * dim; k++) gradv[i][k] = 0.0;
-        gradv[i][k] = 0.0;
     }
     if (rho_flag) {
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) gradr[i][k] = 0.0;
-        gradr[i][k] = 0.0;
     }
     if (energy_flag) {
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) grade[i][k] = 0.0;
-        grade[i][k] = 0.0;
     }
     if (eta_flag) {
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) gradn[i][k] = 0.0;
-        gradn[i][k] = 0.0;
     }
   }
 
@@ -216,10 +215,10 @@ void ComputeRHEOGrad::compute_peratom()
         if (interface_flag) {
           if (fluidi && (!fluidj)) {
             compute_interface->correct_v(vj, vi, j, i);
-            rhoj = compute_interface->correct_rho(j, i);
+            rhoj = compute_interface->correct_rho(j);
           } else if ((!fluidi) && fluidj) {
             compute_interface->correct_v(vi, vj, i, j);
-            rhoi = compute_interface->correct_rho(i, j);
+            rhoi = compute_interface->correct_rho(i);
           } else if ((!fluidi) && (!fluidj)) {
             rhoi = rho0[itype];
             rhoj = rho0[jtype];
@@ -248,34 +247,34 @@ void ComputeRHEOGrad::compute_peratom()
 
         for (a = 0; a < dim; a++) {
           for (b = 0; b < dim; b++) {
-            if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz
+            if (velocity_flag)    // uxx uxy uxz uyx uyy uyz uzx uzy uzz
               gradv[i][a * dim + b] -= vij[a] * Volj * dWij[b];
           }
 
-          if (rho_flag) // P,x  P,y  P,z
+          if (rho_flag)    // P,x  P,y  P,z
             gradr[i][a] -= drho * Volj * dWij[a];
 
-          if (energy_flag) // e,x  e,y  e,z
+          if (energy_flag)    // e,x  e,y  e,z
             grade[i][a] -= de * Volj * dWij[a];
 
-          if (eta_flag) // n,x  n,y  n,z
+          if (eta_flag)    // n,x  n,y  n,z
             gradn[i][a] -= deta * Volj * dWij[a];
         }
 
         if (newton || j < nlocal) {
           for (a = 0; a < dim; a++) {
             for (b = 0; b < dim; b++) {
-              if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz
+              if (velocity_flag)    // uxx uxy uxz uyx uyy uyz uzx uzy uzz
                 gradv[j][a * dim + b] += vij[a] * Voli * dWji[b];
             }
 
-            if (rho_flag) // P,x  P,y  P,z
+            if (rho_flag)    // P,x  P,y  P,z
               gradr[j][a] += drho * Voli * dWji[a];
 
-            if (energy_flag) // e,x  e,y  e,z
+            if (energy_flag)    // e,x  e,y  e,z
               grade[j][a] += de * Voli * dWji[a];
 
-            if (eta_flag) // n,x  n,y  n,z
+            if (eta_flag)    // n,x  n,y  n,z
               gradn[j][a] += deta * Voli * dWji[a];
           }
         }
@@ -295,7 +294,6 @@ void ComputeRHEOGrad::forward_gradients()
   comm->forward_comm(this);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 void ComputeRHEOGrad::forward_fields()
@@ -307,10 +305,9 @@ void ComputeRHEOGrad::forward_fields()
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
+int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
-                                        int pbc_flag, int *pbc)
 {
-  int i,j,k,m;
+  int i, j, k, m;
   int *mask = atom->mask;
   double *rho = atom->rho;
   double *energy = atom->esph;
@@ -333,38 +330,30 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
     if (comm_stage == COMMGRAD) {
 
       if (velocity_flag)
-        for (k = 0; k < dim * dim; k++)
+        for (k = 0; k < dim * dim; k++) buf[m++] = gradv[j][k];
-          buf[m++] = gradv[j][k];
 
       if (rho_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) buf[m++] = gradr[j][k];
-          buf[m++] = gradr[j][k];
 
       if (energy_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) buf[m++] = grade[j][k];
-          buf[m++] = grade[j][k];
 
       if (eta_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) buf[m++] = gradn[j][k];
-          buf[m++] = gradn[j][k];
 
     } else if (comm_stage == COMMFIELD) {
 
       if (velocity_flag) {
         if (remap_v_flag && pbc_flag && (mask[j] & deform_groupbit)) {
-          for (k = 0; k < dim; k++)
+          for (k = 0; k < dim; k++) buf[m++] = v[j][k] + dv[k];
-            buf[m++] = v[j][k] + dv[k];
         } else {
-          for (k = 0; k < dim; k++)
+          for (k = 0; k < dim; k++) buf[m++] = v[j][k];
-            buf[m++] = v[j][k];
         }
       }
 
-      if (rho_flag)
+      if (rho_flag) buf[m++] = rho[j];
-        buf[m++] = rho[j];
 
-      if (energy_flag)
+      if (energy_flag) buf[m++] = energy[j];
-        buf[m++] = energy[j];
     }
   }
   return m;
@@ -376,7 +365,8 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
   double *rho = atom->rho;
-  double *energy = atom->esph;  double **v = atom->v;
+  double *energy = atom->esph;
+  double **v = atom->v;
   int dim = domain->dimension;
 
   m = 0;
@@ -384,31 +374,24 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
   for (i = first; i < last; i++) {
     if (comm_stage == COMMGRAD) {
       if (velocity_flag)
-        for (k = 0; k < dim * dim; k++)
+        for (k = 0; k < dim * dim; k++) gradv[i][k] = buf[m++];
-          gradv[i][k] = buf[m++];
 
       if (rho_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) gradr[i][k] = buf[m++];
-          gradr[i][k] = buf[m++];
 
       if (energy_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) grade[i][k] = buf[m++];
-          grade[i][k] = buf[m++];
 
       if (eta_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) gradn[i][k] = buf[m++];
-          gradn[i][k] = buf[m++];
 
     } else if (comm_stage == COMMFIELD) {
       if (velocity_flag)
-        for (k = 0; k < dim; k++)
+        for (k = 0; k < dim; k++) v[i][k] = buf[m++];
-          v[i][k] = buf[m++];
 
-      if (rho_flag)
+      if (rho_flag) rho[i] = buf[m++];
-        rho[i] = buf[m++];
 
-      if (energy_flag)
+      if (energy_flag) energy[i] = buf[m++];
-        energy[i] = buf[m++];
     }
   }
 }
@@ -417,27 +400,23 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
 
 int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,k,m,last;
+  int i, k, m, last;
   int dim = domain->dimension;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     if (velocity_flag)
-      for (k = 0; k < dim * dim; k++)
+      for (k = 0; k < dim * dim; k++) buf[m++] = gradv[i][k];
-        buf[m++] = gradv[i][k];
 
     if (rho_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) buf[m++] = gradr[i][k];
-        buf[m++] = gradr[i][k];
 
     if (energy_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) buf[m++] = grade[i][k];
-        buf[m++] = grade[i][k];
 
     if (eta_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) buf[m++] = gradn[i][k];
-        buf[m++] = gradn[i][k];
   }
   return m;
 }
@@ -446,27 +425,23 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,k,j,m;
+  int i, k, j, m;
   int dim = domain->dimension;
 
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
     if (velocity_flag)
-      for (k = 0; k < dim * dim; k++)
+      for (k = 0; k < dim * dim; k++) gradv[j][k] += buf[m++];
-        gradv[j][k] += buf[m++];
 
     if (rho_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) gradr[j][k] += buf[m++];
-        gradr[j][k] += buf[m++];
 
     if (energy_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) grade[j][k] += buf[m++];
-        grade[j][k] += buf[m++];
 
     if (eta_flag)
-      for (k = 0; k < dim; k++)
+      for (k = 0; k < dim; k++) gradn[j][k] += buf[m++];
-        gradn[j][k] += buf[m++];
   }
 }
 
@@ -475,17 +450,13 @@ void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf)
 void ComputeRHEOGrad::grow_arrays(int nmax)
 {
   int dim = domain->dimension;
-  if (velocity_flag)
+  if (velocity_flag) memory->grow(gradv, nmax, dim * dim, "rheo:grad_v");
-    memory->grow(gradv, nmax, dim * dim, "rheo:grad_v");
 
-  if (rho_flag)
+  if (rho_flag) memory->grow(gradr, nmax, dim, "rheo:grad_rho");
-    memory->grow(gradr, nmax, dim, "rheo:grad_rho");
 
-  if (energy_flag)
+  if (energy_flag) memory->grow(grade, nmax, dim, "rheo:grad_energy");
-    memory->grow(grade, nmax, dim, "rheo:grad_energy");
 
-  if (eta_flag)
+  if (eta_flag) memory->grow(gradn, nmax, dim, "rheo:grad_eta");
-    memory->grow(gradn, nmax, dim, "rheo:grad_eta");
   nmax_store = nmax;
 }
 
@@ -496,17 +467,13 @@ double ComputeRHEOGrad::memory_usage()
   double bytes = 0.0;
   int dim = domain->dimension;
 
-  if (velocity_flag)
+  if (velocity_flag) bytes = (size_t) nmax_store * dim * dim * sizeof(double);
-    bytes = (size_t) nmax_store * dim * dim * sizeof(double);
 
-  if (rho_flag)
+  if (rho_flag) bytes = (size_t) nmax_store * dim * sizeof(double);
-    bytes = (size_t) nmax_store * dim * sizeof(double);
 
-  if (energy_flag)
+  if (energy_flag) bytes = (size_t) nmax_store * dim * sizeof(double);
-    bytes = (size_t) nmax_store * dim * sizeof(double);
 
-  if (eta_flag)
+  if (eta_flag) bytes = (size_t) nmax_store * dim * sizeof(double);
-    bytes = (size_t) nmax_store * dim * sizeof(double);
 
   return bytes;
 }

src/RHEO/compute_rheo_interface.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -21,18 +20,18 @@
 
 #include "atom.h"
 #include "comm.h"
-#include "domain.h"
 #include "compute_rheo_kernel.h"
+#include "domain.h"
 #include "error.h"
-#include "force.h"
 #include "fix_rheo.h"
 #include "fix_rheo_pressure.h"
+#include "force.h"
 #include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 
 #include <cmath>
 
@@ -45,12 +44,12 @@ static constexpr double EPSILON = 1e-1;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr),
+    Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), rho0(nullptr),
-  rho0(nullptr), norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr),
+    norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr),
-  compute_kernel(nullptr), fix_pressure(nullptr)
+    fix_pressure(nullptr)
 
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
+  if (narg != 3) error->all(FLERR, "Illegal compute rheo/interface command");
 
   comm_forward = 3;
   comm_reverse = 4;
@@ -166,20 +165,18 @@ void ComputeRHEOInterface::compute_peratom()
       if (rsq < cutsq) {
         jtype = type[j];
         fluidj = !(status[j] & PHASECHECK);
-        w = compute_kernel->calc_w_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq));
+        w = compute_kernel->calc_w_quintic(sqrt(rsq));
 
         norm[i] += w;
 
         status_match = 0;
-        if ((fluidi && fluidj) || ((!fluidi) && (!fluidj)))
+        if ((fluidi && fluidj) || ((!fluidi) && (!fluidj))) status_match = 1;
-          status_match = 1;
         if (status_match) {
           chi[i] += w;
         } else {
           if (!fluidi) {
             dot = 0;
-            for (a = 0; a < 3; a++)
+            for (a = 0; a < 3; a++) dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a];
-              dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a];
 
             rho[i] += w * (fix_pressure->calc_pressure(rho[j], jtype) - rho[j] * dot);
             normwf[i] += w;
@@ -193,8 +190,7 @@ void ComputeRHEOInterface::compute_peratom()
           } else {
             if (!fluidj) {
               dot = 0;
-              for (a = 0; a < 3; a++)
+              for (a = 0; a < 3; a++) dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a];
-                dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a];
 
               rho[j] += w * (fix_pressure->calc_pressure(rho[i], itype) + rho[i] * dot);
               normwf[j] += w;
@@ -228,7 +224,8 @@ void ComputeRHEOInterface::compute_peratom()
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
+int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
+                                            int * /*pbc*/)
 {
   int m = 0;
   double *rho = atom->rho;
@@ -293,7 +290,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
     int j = list[i];
     norm[j] += buf[m++];
     chi[j] += buf[m++];
-    if (status[j] & PHASECHECK){
+    if (status[j] & PHASECHECK) {
       normwf[j] += buf[m++];
       rho[j] += buf[m++];
     } else {
@@ -321,7 +318,7 @@ void ComputeRHEOInterface::correct_v(double *v_solid, double *v_fluid, int i_sol
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOInterface::correct_rho(int i_solid, int i_fluid)
+double ComputeRHEOInterface::correct_rho(int i_solid)
 {
   int itype = atom->type[i_solid];
   return MAX(rho0[itype], atom->rho[i_solid]);
@@ -347,28 +344,26 @@ void ComputeRHEOInterface::store_forces()
   if (fixlist.size() != 0) {
     for (const auto &fix : fixlist) {
       for (int i = 0; i < atom->nlocal; i++) {
-        if (rmass) minv = 1.0 / rmass[i];
+        if (rmass)
-        else minv = 1.0 / mass[type[i]];
+          minv = 1.0 / rmass[i];
-        if (mask[i] & fix->groupbit)
-          for (int a = 0; a < 3; a++)
-            fp_store[i][a] = f[i][a] * minv;
         else
-          for (int a = 0; a < 3; a++)
+          minv = 1.0 / mass[type[i]];
-            fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
+        if (mask[i] & fix->groupbit)
+          for (int a = 0; a < 3; a++) fp_store[i][a] = f[i][a] * minv;
+        else
+          for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
       }
     }
   } else {
     if (rmass) {
       for (int i = 0; i < atom->nlocal; i++) {
         minv = 1.0 / rmass[i];
-        for (int a = 0; a < 3; a++)
+        for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
-          fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
       }
     } else {
       for (int i = 0; i < atom->nlocal; i++) {
         minv = 1.0 / mass[type[i]];
-        for (int a = 0; a < 3; a++)
+        for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
-          fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
       }
     }
   }
@@ -388,4 +383,3 @@ double ComputeRHEOInterface::memory_usage()
   double bytes = 3 * nmax_store * sizeof(double);
   return bytes;
 }
-

src/RHEO/compute_rheo_interface.h

@@ -37,7 +37,7 @@ class ComputeRHEOInterface : public Compute {
   void unpack_reverse_comm(int, int *, double *) override;
   double memory_usage() override;
   void correct_v(double *, double *, int, int);
-  double correct_rho(int, int);
+  double correct_rho(int);
   void store_forces();
 
   double *chi, **fp_store;

src/RHEO/compute_rheo_kernel.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -30,31 +29,34 @@
 #include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 #include "pair.h"
 #include "update.h"
 #include "utils.h"
 
 #include <cmath>
+#include <gsl/gsl_cblas.h>
+#include <gsl/gsl_errno.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_cblas.h>
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
 using namespace MathConst;
 using namespace MathExtra;
 
+// max value of Mdim 1 + dim + dim * (dim + 1) / 2 with dim = 3
+static constexpr int MAX_MDIM = 12;
+
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), coordination(nullptr), fix_rheo(nullptr), C(nullptr), C0(nullptr),
+    Compute(lmp, narg, arg), coordination(nullptr), fix_rheo(nullptr), C(nullptr), C0(nullptr),
-  list(nullptr), compute_interface(nullptr)
+    list(nullptr), compute_interface(nullptr)
 {
-  if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command");
+  if (narg != 4) error->all(FLERR, "Illegal compute rheo/kernel command");
 
   kernel_style = utils::inumeric(FLERR, arg[3], false, lmp);
 
@@ -114,13 +116,12 @@ void ComputeRHEOKernel::init()
   cutinv = 1.0 / cut;
   cutsqinv = cutinv * cutinv;
 
-
   if (kernel_style != WENDLANDC4) {
     if (dim == 3) {
       pre_w = 1.0 / (120.0 * MY_PI) * 27.0 * cutsqinv * cutinv;
       pre_wp = pre_w * 3.0 * cutinv;
     } else {
-      pre_w = 7.0 / (478.0 * MY_PI)  * 9 * cutsqinv;
+      pre_w = 7.0 / (478.0 * MY_PI) * 9 * cutsqinv;
       pre_wp = pre_w * 3.0 * cutinv;
     }
   } else {
@@ -157,31 +158,25 @@ int ComputeRHEOKernel::check_corrections(int i)
 {
   // Skip if there were gsl errors for this atom
   if (gsl_error_flag)
-    if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end())
+    if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end()) return 0;
-      return 0;
 
   // Skip if undercoordinated
-  if (coordination[i] < zmin)
+  if (coordination[i] < zmin) return 0;
-    return 0;
 
   // Skip if corrections not yet calculated
-  if (!corrections_calculated)
+  if (!corrections_calculated) return 0;
-    return 0;
 
   return 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_w_self(int i, int j)
+double ComputeRHEOKernel::calc_w_self()
 {
-  double w;
   if (kernel_style == WENDLANDC4)
-    w = calc_w_wendlandc4(i, j, 0.0, 0.0, 0.0, 0.0);
+    return calc_w_wendlandc4(0.0);
   else
-    w = calc_w_quintic(i, j, 0.0, 0.0, 0.0, 0.0);
+    return calc_w_quintic(0.0);
-
-  return w;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -191,8 +186,7 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
   double w = 0.0;
   int corrections_i, corrections_j, corrections;
 
-  if (kernel_style == WENDLANDC4)
+  if (kernel_style == WENDLANDC4) return calc_w_wendlandc4(r);
-    return calc_w_wendlandc4(i, j, delx, dely, delz, r);
 
   if (kernel_style != QUINTIC) {
     corrections_i = check_corrections(i);
@@ -202,10 +196,14 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
     corrections = 0;
   }
 
-  if (!corrections) w = calc_w_quintic(i, j, delx, dely, delz, r);
+  if (!corrections)
-  else if (kernel_style == RK0) w = calc_w_rk0(i, j, delx, dely, delz, r);
+    w = calc_w_quintic(r);
-  else if (kernel_style == RK1) w = calc_w_rk1(i, j, delx, dely, delz, r);
+  else if (kernel_style == RK0)
-  else if (kernel_style == RK2) w = calc_w_rk2(i, j, delx, dely, delz, r);
+    w = calc_w_rk0(i, j, r);
+  else if (kernel_style == RK1)
+    w = calc_w_rk1(i, j, delx, dely, delz, r);
+  else if (kernel_style == RK2)
+    w = calc_w_rk2(i, j, delx, dely, delz, r);
 
   return w;
 }
@@ -217,8 +215,7 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
   double wp;
   int corrections_i, corrections_j;
 
-  if (kernel_style == WENDLANDC4)
+  if (kernel_style == WENDLANDC4) return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji);
-    return calc_dw_wendlandc4(i, j, delx, dely, delz, r, dWij, dWji);
 
   if (kernel_style != QUINTIC) {
     corrections_i = check_corrections(i);
@@ -226,15 +223,15 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
   }
 
   // Calc wp and default dW's, a bit inefficient but can redo later
-  wp = calc_dw_quintic(i, j, delx, dely, delz, r, dWij, dWji);
+  wp = calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
 
   // Overwrite if there are corrections
   if (kernel_style == RK1) {
-    if (corrections_i) calc_dw_rk1(i, j, delx, dely, delz, r, dWij);
+    if (corrections_i) calc_dw_rk1(i, delx, dely, delz, r, dWij);
-    if (corrections_j) calc_dw_rk1(j, i, -delx, -dely, -delz, r, dWji);
+    if (corrections_j) calc_dw_rk1(j, -delx, -dely, -delz, r, dWji);
   } else if (kernel_style == RK2) {
-    if (corrections_i) calc_dw_rk2(i, j, delx, dely, delz, r, dWij);
+    if (corrections_i) calc_dw_rk2(i, delx, dely, delz, r, dWij);
-    if (corrections_j) calc_dw_rk2(j, i, -delx, -dely, -delz, r, dWji);
+    if (corrections_j) calc_dw_rk2(j, -delx, -dely, -delz, r, dWji);
   }
 
   return wp;
@@ -242,30 +239,28 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, double delz, double r)
+double ComputeRHEOKernel::calc_w_quintic(double r)
 {
   double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s;
   s = r * 3.0 * cutinv;
 
-        if (s > 3.0) {
+  if (s > 3.0) { w = 0.0; }
-          w = 0.0;
-        }
 
-        if (s <= 3.0) {
+  if (s <= 3.0) {
-          tmp3 = 3.0 - s;
+    tmp3 = 3.0 - s;
-          tmp3sq = tmp3 * tmp3;
+    tmp3sq = tmp3 * tmp3;
-          w = tmp3sq * tmp3sq * tmp3;
+    w = tmp3sq * tmp3sq * tmp3;
-        }
+  }
-        if (s <= 2.0) {
+  if (s <= 2.0) {
-      tmp2 = 2.0 - s;
+    tmp2 = 2.0 - s;
-      tmp2sq = tmp2 * tmp2;
+    tmp2sq = tmp2 * tmp2;
-      w -= 6.0 * tmp2sq * tmp2sq * tmp2;
+    w -= 6.0 * tmp2sq * tmp2sq * tmp2;
-        }
+  }
-        if (s <= 1.0) {
+  if (s <= 1.0) {
-          tmp1 = 1.0 - s;
+    tmp1 = 1.0 - s;
-          tmp1sq = tmp1 * tmp1;
+    tmp1sq = tmp1 * tmp1;
-          w += 15.0 * tmp1sq * tmp1sq * tmp1;
+    w += 15.0 * tmp1sq * tmp1sq * tmp1;
-        }
+  }
 
   w *= pre_w;
 
@@ -277,15 +272,14 @@ double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely,
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2)
+double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r,
+                                          double *dW1, double *dW2)
 {
   double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv;
 
   s = r * 3.0 * cutinv;
 
-  if (s > 3.0) {
+  if (s > 3.0) { wp = 0.0; }
-    wp = 0.0;
-  }
   if (s <= 3.0) {
     tmp3 = 3.0 - s;
     tmp3sq = tmp3 * tmp3;
@@ -317,7 +311,7 @@ double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double dely, double delz, double r)
+double ComputeRHEOKernel::calc_w_wendlandc4(double r)
 {
   double w, tmp6, s;
   s = r * cutinv;
@@ -340,7 +334,8 @@ double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double de
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2)
+double ComputeRHEOKernel::calc_dw_wendlandc4(double delx, double dely, double delz, double r,
+                                             double *dW1, double *dW2)
 {
   double wp, tmp1, tmp5, tmp6, s, wprinv;
 
@@ -372,11 +367,11 @@ double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double d
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::calc_w_rk0(int i, int j, double delx, double dely, double delz, double r)
+double ComputeRHEOKernel::calc_w_rk0(int i, int j, double r)
 {
   double w;
 
-  w = calc_w_quintic(i, j, delx, dely, delz, r);
+  w = calc_w_quintic(r);
 
   Wij = C0[i] * w;
   Wji = C0[j] * w;
@@ -389,12 +384,12 @@ double ComputeRHEOKernel::calc_w_rk0(int i, int j, double delx, double dely, dou
 double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, double delz, double r)
 {
   int b;
-  double w, dx[3], H[Mdim];
+  double w, dx[3], H[MAX_MDIM];
 
   dx[0] = delx;
   dx[1] = dely;
   dx[2] = delz;
-  w = calc_w_quintic(i, j, delx, dely, delz, r);
+  w = calc_w_quintic(r);
 
   if (dim == 2) {
     H[0] = 1.0;
@@ -408,7 +403,7 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou
   }
   Wij = 0;
   for (b = 0; b < Mdim; b++) {
-    Wij += C[i][0][b] * H[b];  // C columns: 1 x y (z) xx yy (zz)
+    Wij += C[i][0][b] * H[b];    // C columns: 1 x y (z) xx yy (zz)
   }
   Wij *= w;
 
@@ -419,7 +414,7 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou
 
   Wji = 0;
   for (b = 0; b < Mdim; b++) {
-    Wji += C[j][0][b] * H[b];  // C columns: 1 x y (z) xx yy (zz)
+    Wji += C[j][0][b] * H[b];    // C columns: 1 x y (z) xx yy (zz)
   }
   Wji *= w;
 
@@ -431,11 +426,11 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou
 double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, double delz, double r)
 {
   int b;
-  double w, dx[3], H[Mdim];
+  double w, dx[3], H[MAX_MDIM];
   dx[0] = delx;
   dx[1] = dely;
   dx[2] = delz;
-  w = calc_w_quintic(i, j, delx, dely, delz, r);
+  w = calc_w_quintic(r);
 
   if (dim == 2) {
     H[0] = 1.0;
@@ -458,7 +453,7 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou
   }
   Wij = 0;
   for (b = 0; b < Mdim; b++) {
-    Wij += C[i][0][b] * H[b];  // C columns: 1 x y (z) xx yy (zz)
+    Wij += C[i][0][b] * H[b];    // C columns: 1 x y (z) xx yy (zz)
   }
   Wij *= w;
 
@@ -469,7 +464,7 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou
 
   Wji = 0;
   for (b = 0; b < Mdim; b++) {
-    Wji += C[j][0][b] * H[b];  // C columns: 1 x y (z) xx yy (zz)
+    Wji += C[j][0][b] * H[b];    // C columns: 1 x y (z) xx yy (zz)
   }
   Wji *= w;
 
@@ -478,15 +473,16 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou
 
 /* ---------------------------------------------------------------------- */
 
-void ComputeRHEOKernel::calc_dw_rk1(int i, int j, double delx, double dely, double delz, double r, double *dW)
+void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz, double r,
+                                    double *dW)
 {
   int a, b;
-  double w, dx[3], H[Mdim];
+  double w, dx[3], H[MAX_MDIM];
   dx[0] = delx;
   dx[1] = dely;
   dx[2] = delz;
 
-  w = calc_w_quintic(i, j, delx, dely, delz, r);
+  w = calc_w_quintic(r);
 
   //Populate correction basis
   if (dim == 2) {
@@ -506,24 +502,24 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, int j, double delx, double dely, doub
     dW[a] = 0.0;
     for (b = 0; b < Mdim; b++) {
       //First derivative kernels
-      dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z)
+      dW[a] += C[i][1 + a][b] * H[b];    // C columns: 1 x y (z)
     }
     dW[a] *= w;
   }
 }
 
-
 /* ---------------------------------------------------------------------- */
 
-void ComputeRHEOKernel::calc_dw_rk2(int i, int j, double delx, double dely, double delz, double r, double *dW)
+void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz, double r,
+                                    double *dW)
 {
   int a, b;
-  double w, dx[3], H[Mdim];
+  double w, dx[3], H[MAX_MDIM];
   dx[0] = delx;
   dx[1] = dely;
   dx[2] = delz;
 
-  w = calc_w_quintic(i, j, delx, dely, delz, r);
+  w = calc_w_quintic(r);
 
   //Populate correction basis
   if (dim == 2) {
@@ -552,7 +548,7 @@ void ComputeRHEOKernel::calc_dw_rk2(int i, int j, double delx, double dely, doub
     dW[a] = 0.0;
     for (b = 0; b < Mdim; b++) {
       //First derivative kernels
-      dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz)
+      dW[a] += C[i][1 + a][b] * H[b];    // C columns: 1 x y (z) xx yy (zz)
     }
     dW[a] *= w;
   }
@@ -619,14 +615,15 @@ void ComputeRHEOKernel::compute_peratom()
 
         if (rsq < cutsq) {
           r = sqrt(rsq);
-          w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r);
+          w = calc_w_quintic(r);
           rhoj = rho[j];
           if (interface_flag)
-            if (status[j] & PHASECHECK)
+            if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j);
-              rhoj = compute_interface->correct_rho(j, i);
 
-          if (rmass) vj = rmass[j] / rhoj;
+          if (rmass)
-          else vj = mass[type[j]] / rhoj;
+            vj = rmass[j] / rhoj;
+          else
+            vj = mass[type[j]] / rhoj;
           M += w * vj;
         }
       }
@@ -637,7 +634,7 @@ void ComputeRHEOKernel::compute_peratom()
   } else if (correction_order > 0) {
 
     // Moment matrix M and polynomial basis vector cut (1d for gsl compatibility)
-    double H[Mdim], M[Mdim * Mdim];
+    double H[MAX_MDIM], M[MAX_MDIM * MAX_MDIM];
 
     for (ii = 0; ii < inum; ii++) {
       i = ilist[ii];
@@ -650,9 +647,7 @@ void ComputeRHEOKernel::compute_peratom()
 
       // Zero upper-triangle M and cut (will be symmetric):
       for (a = 0; a < Mdim; a++) {
-        for (b = a; b < Mdim; b++) {
+        for (b = a; b < Mdim; b++) { M[a * Mdim + b] = 0; }
-          M[a * Mdim + b] = 0;
-        }
       }
 
       for (jj = 0; jj < jnum; jj++) {
@@ -667,15 +662,16 @@ void ComputeRHEOKernel::compute_peratom()
 
         if (rsq < cutsq) {
           r = sqrt(rsq);
-          w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r);
+          w = calc_w_quintic(r);
 
           rhoj = rho[j];
           if (interface_flag)
-            if (status[j] & PHASECHECK)
+            if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j);
-              rhoj = compute_interface->correct_rho(j, i);
 
-          if (rmass) vj = rmass[j] / rhoj;
+          if (rmass)
-          else vj = mass[type[j]] / rhoj;
+            vj = rmass[j] / rhoj;
+          else
+            vj = mass[type[j]] / rhoj;
 
           //Populate the H-vector of polynomials (2D)
           if (dim == 2) {
@@ -704,18 +700,14 @@ void ComputeRHEOKernel::compute_peratom()
 
           // Populate the upper triangle
           for (a = 0; a < Mdim; a++) {
-            for (b = a; b < Mdim; b++) {
+            for (b = a; b < Mdim; b++) { M[a * Mdim + b] += H[a] * H[b] * w * vj; }
-              M[a * Mdim + b] += H[a] * H[b] * w * vj;
-            }
           }
         }
       }
 
       // Populate the lower triangle from the symmetric entries of M:
       for (a = 0; a < Mdim; a++) {
-        for (b = a; b < Mdim; b++) {
+        for (b = a; b < Mdim; b++) { M[b * Mdim + a] = M[a * Mdim + b]; }
-          M[b * Mdim + a] = M[a * Mdim + b];
-        }
       }
 
       // Skip if undercoordinated
@@ -723,7 +715,7 @@ void ComputeRHEOKernel::compute_peratom()
 
       // Use gsl to get Minv, use Cholesky decomposition since the
       // polynomials are independent, M is symmetrix & positive-definite
-      gM = gsl_matrix_view_array(M,Mdim,Mdim);
+      gM = gsl_matrix_view_array(M, Mdim, Mdim);
       gsl_error = gsl_linalg_cholesky_decomp1(&gM.matrix);
 
       if (gsl_error) {
@@ -738,7 +730,7 @@ void ComputeRHEOKernel::compute_peratom()
         continue;
       }
 
-      gsl_linalg_cholesky_invert(&gM.matrix);   //M is now M^-1
+      gsl_linalg_cholesky_invert(&gM.matrix);    //M is now M^-1
 
       // Correction coefficients are columns of M^-1 multiplied by an appropriate coefficient
       // Solve the linear system several times to get coefficientns
@@ -759,16 +751,15 @@ void ComputeRHEOKernel::compute_peratom()
 
       // Pack coefficients into C
       for (a = 0; a < Mdim; a++) {
-        C[i][0][a] = M[a * Mdim + 0]; // all rows of column 0
+        C[i][0][a] = M[a * Mdim + 0];    // all rows of column 0
         for (b = 0; b < dim; b++) {
           //First derivatives
           C[i][1 + b][a] = -M[a * Mdim + b + 1] * cutinv;
           // columns 1-2 (2D)  or 1-3 (3D)
 
           //Second derivatives
-          if (kernel_style == RK2)
+          if (kernel_style == RK2) C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * cutsqinv;
-            C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * cutsqinv;
+          // columns 3-4 (2D) or 4-6 (3D)
-            // columns 3-4 (2D) or 4-6 (3D)
         }
       }
     }
@@ -780,7 +771,6 @@ void ComputeRHEOKernel::compute_peratom()
   comm->forward_comm(this);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 void ComputeRHEOKernel::compute_coordination()
@@ -819,8 +809,7 @@ void ComputeRHEOKernel::compute_coordination()
       dx[2] = ztmp - x[j][2];
       rsq = lensq3(dx);
 
-      if (rsq < cutsq)
+      if (rsq < cutsq) coordination[i] += 1;
-        coordination[i] += 1;
     }
   }
 
@@ -847,8 +836,8 @@ void ComputeRHEOKernel::grow_arrays(int nmax)
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
+int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                         int * /*pbc*/)
 {
   int m = 0;
   for (int i = 0; i < n; i++) {
@@ -860,8 +849,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
         buf[m++] = C0[j];
       } else {
         for (int a = 0; a < ncor; a++)
-          for (int b = 0; b < Mdim; b++)
+          for (int b = 0; b < Mdim; b++) buf[m++] = C[j][a][b];
-            buf[m++] = C[j][a][b];
       }
     }
   }
@@ -882,8 +870,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
         C0[i] = buf[m++];
       } else {
         for (int a = 0; a < ncor; a++)
-          for (int b = 0; b < Mdim; b++)
+          for (int b = 0; b < Mdim; b++) C[i][a][b] = buf[m++];
-            C[i][a][b] = buf[m++];
       }
     }
   }

src/RHEO/compute_rheo_kernel.h

@@ -36,13 +36,13 @@ class ComputeRHEOKernel : public Compute {
   void unpack_forward_comm(int, int, double *) override;
   double memory_usage() override;
   void compute_coordination();
-  double calc_w_self(int,int);
+  double calc_w_self();
-  double calc_w(int,int,double,double,double,double);
+  double calc_w(int, int, double, double, double, double);
-  double calc_dw(int,int,double,double,double,double);
+  double calc_dw(int, int, double, double, double, double);
-  double calc_w_quintic(int,int,double,double,double,double);
+  double calc_w_quintic(double);
-  double calc_dw_quintic(int,int,double,double,double,double,double *,double *);
+  double calc_dw_quintic(double, double, double, double, double *, double *);
-  double calc_w_wendlandc4(int,int,double,double,double,double);
+  double calc_w_wendlandc4(double);
-  double calc_dw_wendlandc4(int,int,double,double,double,double,double *,double *);
+  double calc_dw_wendlandc4(double, double, double, double, double *, double *);
   void grow_arrays(int);
 
   double dWij[3], dWji[3], Wij, Wji;
@@ -68,14 +68,12 @@ class ComputeRHEOKernel : public Compute {
 
   int check_corrections(int);
 
-  double calc_w_rk0(int,int,double,double,double,double);
+  double calc_w_rk0(int, int, double);
-  double calc_w_rk1(int,int,double,double,double,double);
+  double calc_w_rk1(int, int, double, double, double, double);
-  double calc_w_rk2(int,int,double,double,double,double);
+  double calc_w_rk2(int, int, double, double, double, double);
-  void calc_dw_rk1(int,int,double,double,double,double,double *);
+  void calc_dw_rk1(int, double, double, double, double, double *);
-  void calc_dw_rk2(int,int,double,double,double,double,double *);
+  void calc_dw_rk2(int, double, double, double, double, double *);
 };
-
 }    // namespace LAMMPS_NS
-
 #endif
 #endif

src/RHEO/compute_rheo_property_atom.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -21,11 +20,11 @@
 
 #include "atom.h"
 #include "atom_vec.h"
+#include "compute_rheo_grad.h"
 #include "compute_rheo_interface.h"
 #include "compute_rheo_kernel.h"
 #include "compute_rheo_surface.h"
 #include "compute_rheo_vshift.h"
-#include "compute_rheo_grad.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
@@ -46,12 +45,12 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), avec_index(nullptr), col_index(nullptr), col_t_index(nullptr),
+    Compute(lmp, narg, arg), avec_index(nullptr), col_index(nullptr), col_t_index(nullptr),
-  buf(nullptr), pack_choice(nullptr), fix_rheo(nullptr), fix_pressure(nullptr),
+    buf(nullptr), pack_choice(nullptr), fix_rheo(nullptr), fix_pressure(nullptr),
-  fix_thermal(nullptr),  compute_interface(nullptr), compute_kernel(nullptr),
+    fix_thermal(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
-  compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr)
+    compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr)
 {
-  if (narg < 4)  utils::missing_cmd_args(FLERR, "compute property/atom", error);
+  if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error);
 
   peratom_flag = 1;
   int dim = domain->dimension;
@@ -74,8 +73,10 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
     }
   }
 
-  if (nvalues == 1) size_peratom_cols = 0;
+  if (nvalues == 1)
-  else size_peratom_cols = nvalues;
+    size_peratom_cols = 0;
+  else
+    size_peratom_cols = nvalues;
 
   pressure_flag = thermal_flag = interface_flag = 0;
   surface_flag = shift_flag = shell_flag = 0;
@@ -173,15 +174,19 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
 void ComputeRHEOPropertyAtom::init()
 {
   auto fixes = modify->get_fix_by_style("^rheo$");
-  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/property/atom");
+  if (fixes.size() == 0)
+    error->all(FLERR, "Need to define fix rheo to use compute rheo/property/atom");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
   if (interface_flag && !(fix_rheo->interface_flag))
-    error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo");
+    error->all(FLERR,
+               "Cannot request interfacial property without corresponding option in fix rheo");
   if (surface_flag && !(fix_rheo->surface_flag))
     error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo");
   if (shift_flag && !(fix_rheo->shift_flag))
-    error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo");
+    error->all(
+        FLERR,
+        "Cannot request velocity shifting property without corresponding option in fix rheo");
   if (thermal_flag && !(fix_rheo->thermal_flag))
     error->all(FLERR, "Cannot request thermal property without fix rheo/thermal");
 
@@ -237,10 +242,11 @@ void ComputeRHEOPropertyAtom::compute_peratom()
     buf = vector_atom;
     (this->*pack_choice[0])(0);
   } else {
-    if (nmax) buf = &array_atom[0][0];
+    if (nmax)
-    else buf = nullptr;
+      buf = &array_atom[0][0];
-    for (int n = 0; n < nvalues; n++)
+    else
-      (this->*pack_choice[n])(n);
+      buf = nullptr;
+    for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n);
   }
 }
 
@@ -250,7 +256,7 @@ void ComputeRHEOPropertyAtom::compute_peratom()
 
 double ComputeRHEOPropertyAtom::memory_usage()
 {
-  double bytes = (double)nmax * nvalues * sizeof(double);
+  double bytes = (double) nmax * nvalues * sizeof(double);
   return bytes;
 }
 
@@ -269,8 +275,10 @@ void ComputeRHEOPropertyAtom::pack_phase(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = (status[i] & PHASECHECK);
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = (status[i] & PHASECHECK);
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -284,8 +292,10 @@ void ComputeRHEOPropertyAtom::pack_status(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = status[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = status[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -299,8 +309,10 @@ void ComputeRHEOPropertyAtom::pack_chi(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = chi[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = chi[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -334,8 +346,10 @@ void ComputeRHEOPropertyAtom::pack_surface_r(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = rsurface[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = rsurface[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -349,8 +363,10 @@ void ComputeRHEOPropertyAtom::pack_surface_divr(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = divr[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = divr[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -365,8 +381,10 @@ void ComputeRHEOPropertyAtom::pack_surface_n(int n)
   int index = col_index[n];
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = nsurface[i][index];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = nsurface[i][index];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -380,8 +398,10 @@ void ComputeRHEOPropertyAtom::pack_coordination(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = coordination[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = coordination[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -395,8 +415,10 @@ void ComputeRHEOPropertyAtom::pack_cv(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i, type[i]);
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = fix_thermal->calc_cv(type[i]);
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -411,8 +433,10 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = fix_pressure->calc_pressure(rho[i], type[i]);
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = fix_pressure->calc_pressure(rho[i], type[i]);
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -432,8 +456,10 @@ void ComputeRHEOPropertyAtom::pack_viscous_stress(int n)
   int index_transpose = b * dim + a;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]);
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]);
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -462,7 +488,8 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n)
       else
         p = 0.0;
       buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p;
-    } else buf[n] = 0.0;
+    } else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -476,8 +503,10 @@ void ComputeRHEOPropertyAtom::pack_nbond_shell(int n)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = nbond[i];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = nbond[i];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -492,8 +521,10 @@ void ComputeRHEOPropertyAtom::pack_shift_v(int n)
   int index = col_index[n];
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = vshift[i][index];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = vshift[i][index];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -508,8 +539,10 @@ void ComputeRHEOPropertyAtom::pack_gradv(int n)
   int index = col_index[n];
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) buf[n] = gradv[i][index];
+    if (mask[i] & groupbit)
-    else buf[n] = 0.0;
+      buf[n] = gradv[i][index];
+    else
+      buf[n] = 0.0;
     n += nvalues;
   }
 }
@@ -523,7 +556,7 @@ void ComputeRHEOPropertyAtom::pack_atom_style(int n)
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack pack_function)
+int ComputeRHEOPropertyAtom::add_tensor_component(char *option, int i, FnPtrPack pack_function)
 {
   int shift;
   int dim = domain->dimension;
@@ -547,11 +580,15 @@ int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack
       index = 2;
       if (dim == 2) dim_error = 1;
     } else if (utils::strmatch(option, "yx$")) {
-      if (dim == 2) index = 2;
+      if (dim == 2)
-      else index = 3;
+        index = 2;
+      else
+        index = 3;
     } else if (utils::strmatch(option, "yy$")) {
-      if (dim == 2) index = 3;
+      if (dim == 2)
-      else index = 4;
+        index = 3;
+      else
+        index = 4;
     } else if (utils::strmatch(option, "yz$")) {
       index = 5;
       if (dim == 2) dim_error = 1;
@@ -581,7 +618,7 @@ int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOPropertyAtom::add_vector_component(char* option, int i, FnPtrPack pack_function)
+int ComputeRHEOPropertyAtom::add_vector_component(char *option, int i, FnPtrPack pack_function)
 {
   int shift;
   int dim = domain->dimension;

src/RHEO/compute_rheo_property_atom.h

@@ -61,8 +61,8 @@ class ComputeRHEOPropertyAtom : public Compute {
   void pack_nbond_shell(int);
   void pack_atom_style(int);
 
-  int add_vector_component(char*, int, FnPtrPack);
+  int add_vector_component(char *, int, FnPtrPack);
-  int add_tensor_component(char*, int, FnPtrPack);
+  int add_tensor_component(char *, int, FnPtrPack);
 
   class FixRHEO *fix_rheo;
   class FixRHEOPressure *fix_pressure;
@@ -73,11 +73,9 @@ class ComputeRHEOPropertyAtom : public Compute {
   class ComputeRHEOSurface *compute_surface;
   class ComputeRHEOVShift *compute_vshift;
   class ComputeRHEOGrad *compute_grad;
-
 };
 
 }    // namespace LAMMPS_NS
 
 #endif
 #endif
-

src/RHEO/compute_rheo_rho_sum.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -25,18 +24,18 @@
 #include "error.h"
 #include "fix_rheo.h"
 #include "force.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEORhoSum::ComputeRHEORhoSum(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr)
+    Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr)
 {
-  if (narg != 4) error->all(FLERR,"Illegal compute RHEO/rho command");
+  if (narg != 4) error->all(FLERR, "Illegal compute RHEO/rho command");
 
   self_mass_flag = utils::bnumeric(FLERR, arg[3], false, lmp);
 
@@ -69,7 +68,6 @@ void ComputeRHEORhoSum::init_list(int /*id*/, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
-
 void ComputeRHEORhoSum::compute_peratom()
 {
   int i, j, ii, jj, inum, jnum;
@@ -94,9 +92,11 @@ void ComputeRHEORhoSum::compute_peratom()
 
   // initialize arrays, local with quintic self-contribution, ghosts are zeroed
   for (i = 0; i < nlocal; i++) {
-    w = compute_kernel->calc_w_self(i, i);
+    w = compute_kernel->calc_w_self();
-    if (rmass) rho[i] = w * rmass[i];
+    if (rmass)
-    else rho[i] = w * mass[type[i]];
+      rho[i] = w * rmass[i];
+    else
+      rho[i] = w * mass[type[i]];
   }
 
   for (i = nlocal; i < nall; i++) rho[i] = 0.0;
@@ -123,22 +123,18 @@ void ComputeRHEORhoSum::compute_peratom()
         if (rmass) {
           if (self_mass_flag) {
             rho[i] += w * rmass[i];
-            if (newton || j < nlocal)
+            if (newton || j < nlocal) rho[j] += w * rmass[j];
-              rho[j] += w * rmass[j];
           } else {
             rho[i] += w * rmass[j];
-            if (newton || j < nlocal)
+            if (newton || j < nlocal) rho[j] += w * rmass[i];
-              rho[j] += w * rmass[i];
           }
         } else {
           if (self_mass_flag) {
             rho[i] += w * mass[type[i]];
-            if (newton || j < nlocal)
+            if (newton || j < nlocal) rho[j] += w * mass[type[j]];
-              rho[j] += w * mass[type[j]];
           } else {
             rho[i] += w * mass[type[j]];
-            if (newton || j < nlocal)
+            if (newton || j < nlocal) rho[j] += w * mass[type[i]];
-              rho[j] += w * mass[type[i]];
           }
         }
       }
@@ -151,8 +147,8 @@ void ComputeRHEORhoSum::compute_peratom()
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf,
+int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                         int * /*pbc*/)
 {
   int i, j, m;
   double *rho = atom->rho;
@@ -173,9 +169,7 @@ void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf)
 
   m = 0;
   last = first + n;
-  for (i = first; i < last; i++) {
+  for (i = first; i < last; i++) { rho[i] = buf[m++]; }
-    rho[i] = buf[m++];
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -187,9 +181,7 @@ int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf)
 
   m = 0;
   last = first + n;
-  for (i = first; i < last; i++) {
+  for (i = first; i < last; i++) { buf[m++] = rho[i]; }
-    buf[m++] = rho[i];
-  }
   return m;
 }
 

src/RHEO/compute_rheo_surface.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -29,9 +28,9 @@
 #include "force.h"
 #include "math_extra.h"
 #include "memory.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
@@ -43,11 +42,11 @@ static constexpr double EPSILON = 1e-10;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), nsurface(nullptr), rsurface(nullptr), divr(nullptr), fix_rheo(nullptr),
+    Compute(lmp, narg, arg), nsurface(nullptr), rsurface(nullptr), divr(nullptr), fix_rheo(nullptr),
-  rho0(nullptr), B(nullptr), gradC(nullptr), list(nullptr), compute_kernel(nullptr),
+    rho0(nullptr), B(nullptr), gradC(nullptr), list(nullptr), compute_kernel(nullptr),
-  compute_interface(nullptr)
+    compute_interface(nullptr)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute RHEO/SURFACE command");
+  if (narg != 3) error->all(FLERR, "Illegal compute RHEO/SURFACE command");
 
   int dim = domain->dimension;
   comm_forward = 2;
@@ -123,8 +122,7 @@ void ComputeRHEOSurface::compute_peratom()
   firstneigh = list->firstneigh;
 
   // Grow and zero arrays
-  if (nmax_store < atom->nmax)
+  if (nmax_store < atom->nmax) grow_arrays(atom->nmax);
-    grow_arrays(atom->nmax);
 
   size_t nbytes = nmax_store * sizeof(double);
   memset(&divr[0], 0, nbytes);
@@ -133,7 +131,6 @@ void ComputeRHEOSurface::compute_peratom()
   memset(&gradC[0][0], 0, dim * dim * nbytes);
   memset(&B[0][0], 0, dim * dim * nbytes);
 
-
   // loop over neighbors to calculate the average orientation of neighbors
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -165,9 +162,9 @@ void ComputeRHEOSurface::compute_peratom()
         // Add corrections for walls
         if (interface_flag) {
           if (fluidi && (!fluidj)) {
-            rhoj = compute_interface->correct_rho(j, i);
+            rhoj = compute_interface->correct_rho(j);
           } else if ((!fluidi) && fluidj) {
-            rhoi = compute_interface->correct_rho(i, j);
+            rhoi = compute_interface->correct_rho(i);
           } else if ((!fluidi) && (!fluidj)) {
             rhoi = rho0[itype];
             rhoj = rho0[jtype];
@@ -181,7 +178,7 @@ void ComputeRHEOSurface::compute_peratom()
           Voli = mass[itype] / rhoi;
           Volj = mass[jtype] / rhoj;
         }
-        compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
+        compute_kernel->calc_dw_quintic(dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
 
         for (a = 0; a < dim; a++) {
           divr[i] -= dWij[a] * dx[a] * Volj;
@@ -189,7 +186,7 @@ void ComputeRHEOSurface::compute_peratom()
         }
 
         if ((j < nlocal) || newton) {
-          for (a = 0; a < dim; a++){
+          for (a = 0; a < dim; a++) {
             divr[j] += dWji[a] * dx[a] * Voli;
             gradC[j][a] += dWji[a] * Voli;
           }
@@ -211,12 +208,10 @@ void ComputeRHEOSurface::compute_peratom()
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       maggC = 0.0;
-      for (a = 0;a < dim; a++)
+      for (a = 0; a < dim; a++) maggC += gradC[i][a] * gradC[i][a];
-        maggC += gradC[i][a] * gradC[i][a];
       maggC = sqrt(maggC) + EPSILON;
       maggC = 1.0 / maggC;
-      for (a = 0; a < dim; a++)
+      for (a = 0; a < dim; a++) nsurface[i][a] = -gradC[i][a] * maggC;
-        nsurface[i][a] = -gradC[i][a] * maggC;
     }
   }
 
@@ -233,8 +228,7 @@ void ComputeRHEOSurface::compute_peratom()
         test = coordination[i] < threshold_z;
 
       // Treat nonfluid particles as bulk
-      if (status[i] & PHASECHECK)
+      if (status[i] & PHASECHECK) test = 0;
-        test = 0;
 
       if (test) {
         if (coordination[i] < threshold_splash)
@@ -249,7 +243,6 @@ void ComputeRHEOSurface::compute_peratom()
     }
   }
 
-
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     xtmp = x[i][0];
@@ -277,18 +270,17 @@ void ComputeRHEOSurface::compute_peratom()
             status[i] |= STATUS_LAYER;
           }
 
-          if (status[j] & STATUS_SURFACE)
+          if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
-            rsurface[i] = MIN(rsurface[i], sqrt(rsq));
         }
 
         if (fluidj && (j < nlocal || newton)) {
-          if ((status[j] & STATUS_BULK) && (status[j] & PHASECHECK) && (status[i] & STATUS_SURFACE)) {
+          if ((status[j] & STATUS_BULK) && (status[j] & PHASECHECK) &&
+              (status[i] & STATUS_SURFACE)) {
             status[j] &= SURFACEMASK;
             status[j] |= STATUS_LAYER;
           }
 
-          if (status[i] & STATUS_SURFACE)
+          if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
-            rsurface[j] = MIN(rsurface[j], sqrt(rsq));
         }
       }
     }
@@ -299,8 +291,7 @@ void ComputeRHEOSurface::compute_peratom()
   for (i = 0; i < nall; i++) {
     if (mask[i] & groupbit) {
       if (!(status[i] & STATUS_SURFACE))
-        for (a = 0; a < dim; a++)
+        for (a = 0; a < dim; a++) nsurface[i][a] = 0.0;
-          nsurface[i][a] = 0.0;
     }
   }
 
@@ -323,9 +314,8 @@ int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
   for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
       buf[m++] = divr[i];
-      for (int a = 0; a < dim; a ++ )
+      for (int a = 0; a < dim; a++)
-        for (int b = 0; b < dim; b ++)
+        for (int b = 0; b < dim; b++) buf[m++] = gradC[i][a * dim + b];
-          buf[m++] = gradC[i][a * dim + b];
     } else if (comm_stage == 1) {
       buf[m++] = (double) status[i];
       buf[m++] = rsurface[i];
@@ -345,9 +335,8 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
     int j = list[i];
     if (comm_stage == 0) {
       divr[j] += buf[m++];
-      for (int a = 0; a < dim; a ++ )
+      for (int a = 0; a < dim; a++)
-        for (int b = 0; b < dim; b ++)
+        for (int b = 0; b < dim; b++) gradC[j][a * dim + b] += buf[m++];
-          gradC[j][a * dim + b] += buf[m++];
     } else if (comm_stage == 1) {
       auto tmp1 = (int) buf[m++];
       if ((status[j] & STATUS_BULK) && (tmp1 & STATUS_LAYER)) {
@@ -360,11 +349,10 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
   }
 }
 
-
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
+int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                          int * /*pbc*/)
 {
   int *status = atom->rheo_status;
   int m = 0;

src/RHEO/compute_rheo_vshift.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -29,9 +28,9 @@
 #include "fix_rheo.h"
 #include "force.h"
 #include "memory.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 
 #include "update.h"
 
@@ -41,10 +40,10 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr),
+    Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr),
-  compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr)
+    compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");
+  if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command");
 
   comm_reverse = 3;
   surface_flag = 0;
@@ -125,8 +124,7 @@ void ComputeRHEOVShift::compute_peratom()
   }
 
   for (i = 0; i < nall; i++)
-    for (a = 0; a < dim; a++)
+    for (a = 0; a < dim; a++) vshift[i][a] = 0.0;
-      vshift[i][a] = 0.0;
 
   for (a = 0; a < 3; a++) {
     vi[a] = 0.0;
@@ -141,8 +139,10 @@ void ComputeRHEOVShift::compute_peratom()
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    if (rmass) imass = rmass[i];
+    if (rmass)
-    else imass = mass[itype];
+      imass = rmass[i];
+    else
+      imass = mass[itype];
     fluidi = !(status[i] & PHASECHECK);
 
     for (jj = 0; jj < jnum; jj++) {
@@ -162,8 +162,10 @@ void ComputeRHEOVShift::compute_peratom()
 
       if (rsq < cutsq) {
         jtype = type[j];
-        if (rmass) jmass = rmass[j];
+        if (rmass)
-        else jmass = mass[jtype];
+          jmass = rmass[j];
+        else
+          jmass = mass[jtype];
 
         r = sqrt(rsq);
         rinv = 1 / r;
@@ -180,10 +182,10 @@ void ComputeRHEOVShift::compute_peratom()
         if (interface_flag) {
           if (fluidi && (!fluidj)) {
             compute_interface->correct_v(vj, vi, j, i);
-            rhoj = compute_interface->correct_rho(j, i);
+            rhoj = compute_interface->correct_rho(j);
           } else if ((!fluidi) && fluidj) {
             compute_interface->correct_v(vi, vj, i, j);
-            rhoi = compute_interface->correct_rho(i, j);
+            rhoi = compute_interface->correct_rho(i);
           } else if ((!fluidi) && (!fluidj)) {
             rhoi = rho0[itype];
             rhoj = rho0[jtype];
@@ -195,7 +197,7 @@ void ComputeRHEOVShift::compute_peratom()
 
         wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2], r);
         w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
-        w0 = compute_kernel->calc_w(i, j, 0, 0, 0, cutthird); // dx, dy, dz irrelevant
+        w0 = compute_kernel->calc_w(i, j, 0, 0, 0, cutthird);    // dx, dy, dz irrelevant
         w4 = w * w * w * w / (w0 * w0 * w0 * w0);
         dr = -2 * cutthird * (1 + 0.2 * w4) * wp * rinv;
 
@@ -225,7 +227,6 @@ void ComputeRHEOVShift::compute_peratom()
   if (newton_pair) comm->reverse_comm(this);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 void ComputeRHEOVShift::correct_surfaces()
@@ -300,7 +301,7 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,j,m;
+  int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {

src/RHEO/fix_rheo.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -23,9 +22,9 @@
 #include "citeme.h"
 #include "compute_rheo_grad.h"
 #include "compute_rheo_interface.h"
-#include "compute_rheo_surface.h"
 #include "compute_rheo_kernel.h"
 #include "compute_rheo_rho_sum.h"
+#include "compute_rheo_surface.h"
 #include "compute_rheo_vshift.h"
 #include "domain.h"
 #include "error.h"
@@ -39,20 +38,23 @@ using namespace LAMMPS_NS;
 using namespace RHEO_NS;
 using namespace FixConst;
 
+#if 0
+// publication was removed from documentation
 static const char cite_rheo[] =
-  "@article{PalermoInPrep,\n"
+    "@article{PalermoInPrep,\n"
-  " journal = {in prep},\n"
+    " journal = {in prep},\n"
-  " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n"
+    " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n"
-  " year = {2024},\n"
+    " year = {2024},\n"
-  " author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n"
+    " author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n"
-  "}\n\n";
+    "}\n\n";
+#endif
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), compute_kernel(nullptr),
+    Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr),
-  compute_interface(nullptr), compute_surface(nullptr), compute_rhosum(nullptr),
+    compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr),
-  compute_vshift(nullptr)
+    compute_rhosum(nullptr), compute_vshift(nullptr)
 {
   time_integrate = 1;
 
@@ -78,44 +80,42 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     csq[i] = 1.0;
   }
 
-  if (igroup != 0)
+  if (igroup != 0) error->all(FLERR, "fix rheo command requires group all");
-    error->all(FLERR, "fix rheo command requires group all");
 
   if (atom->pressure_flag != 1)
     error->all(FLERR, "fix rheo command requires atom_style with pressure");
-  if (atom->rho_flag != 1)
+  if (atom->rho_flag != 1) error->all(FLERR, "fix rheo command requires atom_style with density");
-    error->all(FLERR, "fix rheo command requires atom_style with density");
   if (atom->viscosity_flag != 1)
     error->all(FLERR, "fix rheo command requires atom_style with viscosity");
   if (atom->rheo_status_flag != 1)
     error->all(FLERR, "fix rheo command requires atom_style with status");
 
-  if (narg < 5)
+  if (narg < 5) utils::missing_cmd_args(FLERR, "fix rheo", error);
-    utils::missing_cmd_args(FLERR, "fix rheo", error);
 
   cut = utils::numeric(FLERR, arg[3], false, lmp);
   if (strcmp(arg[4], "quintic") == 0) {
-      kernel_style = QUINTIC;
+    kernel_style = QUINTIC;
   } else if (strcmp(arg[4], "wendland/c4") == 0) {
-      kernel_style = WENDLANDC4;
+    kernel_style = WENDLANDC4;
   } else if (strcmp(arg[4], "RK0") == 0) {
-      kernel_style = RK0;
+    kernel_style = RK0;
   } else if (strcmp(arg[4], "RK1") == 0) {
-      kernel_style = RK1;
+    kernel_style = RK1;
   } else if (strcmp(arg[4], "RK2") == 0) {
-      kernel_style = RK2;
+    kernel_style = RK2;
-  } else error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]);
+  } else
+    error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]);
   zmin_kernel = utils::numeric(FLERR, arg[5], false, lmp);
 
   int iarg = 6;
-  while (iarg < narg){
+  while (iarg < narg) {
     if (strcmp(arg[iarg], "shift") == 0) {
       shift_flag = 1;
     } else if (strcmp(arg[iarg], "thermal") == 0) {
       thermal_flag = 1;
     } else if (strcmp(arg[iarg], "surface/detection") == 0) {
       surface_flag = 1;
-      if(iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo surface/detection", error);
+      if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo surface/detection", error);
       if (strcmp(arg[iarg + 1], "coordination") == 0) {
         surface_style = COORDINATION;
         zmin_surface = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
@@ -137,8 +137,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       self_mass_flag = 1;
     } else if (strcmp(arg[iarg], "density") == 0) {
       if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error);
-      for (i = 1; i <= n; i++)
+      for (i = 1; i <= n; i++) rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
-        rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
       iarg += n;
     } else if (strcmp(arg[iarg], "speed/sound") == 0) {
       if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo speed/sound", error);
@@ -156,7 +155,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   if (self_mass_flag && (!rhosum_flag))
     error->all(FLERR, "Cannot use self/mass setting without rho/sum");
 
+#if 0
   if (lmp->citeme) lmp->citeme->add(cite_rheo);
+#endif
 }
 
 /* ---------------------------------------------------------------------- */
@@ -174,15 +175,14 @@ FixRHEO::~FixRHEO()
   memory->destroy(rho0);
 }
 
-
 /* ----------------------------------------------------------------------
   Create necessary internal computes
 ------------------------------------------------------------------------- */
 
 void FixRHEO::post_constructor()
 {
-  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(
+  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(
-    fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style)));
+      modify->add_compute(fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style)));
   compute_kernel->fix_rheo = this;
 
   std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity";
@@ -191,26 +191,26 @@ void FixRHEO::post_constructor()
   compute_grad->fix_rheo = this;
 
   if (rhosum_flag) {
-    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(
+    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(
-      fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag)));
+        modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag)));
     compute_rhosum->fix_rheo = this;
   }
 
   if (shift_flag) {
-    compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(
+    compute_vshift =
-      "rheo_vshift all RHEO/VSHIFT"));
+        dynamic_cast<ComputeRHEOVShift *>(modify->add_compute("rheo_vshift all RHEO/VSHIFT"));
     compute_vshift->fix_rheo = this;
   }
 
   if (interface_flag) {
-    compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(
+    compute_interface = dynamic_cast<ComputeRHEOInterface *>(
-      "rheo_interface all RHEO/INTERFACE"));
+        modify->add_compute("rheo_interface all RHEO/INTERFACE"));
     compute_interface->fix_rheo = this;
   }
 
   if (surface_flag) {
-    compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(
+    compute_surface =
-      "rheo_surface all RHEO/SURFACE"));
+        dynamic_cast<ComputeRHEOSurface *>(modify->add_compute("rheo_surface all RHEO/SURFACE"));
     compute_surface->fix_rheo = this;
   }
 }
@@ -237,23 +237,23 @@ void FixRHEO::init()
     error->all(FLERR, "Can only specify one instance of fix rheo");
 
   if (atom->rheo_status_flag != 1)
-    error->all(FLERR,"fix rheo command requires atom property status");
+    error->all(FLERR, "fix rheo command requires atom property status");
-  if (atom->rho_flag != 1)
+  if (atom->rho_flag != 1) error->all(FLERR, "fix rheo command requires atom property rho");
-    error->all(FLERR,"fix rheo command requires atom property rho");
   if (atom->pressure_flag != 1)
-    error->all(FLERR,"fix rheo command requires atom property pressure");
+    error->all(FLERR, "fix rheo command requires atom property pressure");
   if (atom->viscosity_flag != 1)
-    error->all(FLERR,"fix rheo command requires atom property viscosity");
+    error->all(FLERR, "fix rheo command requires atom property viscosity");
 
   if (thermal_flag) {
     if (atom->esph_flag != 1)
-      error->all(FLERR,"fix rheo command requires atom property esph with thermal setting");
+      error->all(FLERR, "fix rheo command requires atom property esph with thermal setting");
     if (atom->temperature_flag != 1)
-      error->all(FLERR,"fix rheo command requires atom property temperature with thermal setting");
+      error->all(FLERR, "fix rheo command requires atom property temperature with thermal setting");
     if (atom->heatflow_flag != 1)
-      error->all(FLERR,"fix rheo command requires atom property heatflow with thermal setting");
+      error->all(FLERR, "fix rheo command requires atom property heatflow with thermal setting");
     if (atom->conductivity_flag != 1)
-      error->all(FLERR,"fix rheo command requires atom property conductivity with thermal setting");
+      error->all(FLERR,
+                 "fix rheo command requires atom property conductivity with thermal setting");
   }
 }
 
@@ -281,14 +281,11 @@ void FixRHEO::setup(int /*vflag*/)
 {
   // Confirm all accessory fixes are defined
   // Note: fixes set this flag in setup_pre_force()
-  if (!viscosity_fix_defined)
+  if (!viscosity_fix_defined) error->all(FLERR, "Missing fix rheo/viscosity");
-    error->all(FLERR, "Missing fix rheo/viscosity");
 
-  if (!pressure_fix_defined)
+  if (!pressure_fix_defined) error->all(FLERR, "Missing fix rheo/pressure");
-    error->all(FLERR, "Missing fix rheo/pressure");
 
-  if(thermal_flag && !thermal_fix_defined)
+  if (thermal_flag && !thermal_fix_defined) error->all(FLERR, "Missing fix rheo/thermal");
-    error->all(FLERR, "Missing fix rheo/thermal");
 
   // Reset to zero for future runs
   thermal_fix_defined = 0;
@@ -296,8 +293,7 @@ void FixRHEO::setup(int /*vflag*/)
   pressure_fix_defined = 0;
   oxidation_fix_defined = 0;
 
-  if (rhosum_flag)
+  if (rhosum_flag) compute_rhosum->compute_peratom();
-    compute_rhosum->compute_peratom();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -321,15 +317,13 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
   double **vshift;
-  if (shift_flag)
+  if (shift_flag) vshift = compute_vshift->vshift;
-    vshift = compute_vshift->vshift;
 
   int nlocal = atom->nlocal;
   int rmass_flag = atom->rmass_flag;
   int dim = domain->dimension;
 
-  if (igroup == atom->firstgroup)
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
-    nlocal = atom->nfirst;
 
   //Density Half-step
   for (i = 0; i < nlocal; i++) {
@@ -349,7 +343,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   }
 
   // Update gradients and interpolate solid properties
-  compute_grad->forward_fields(); // also forwards v and rho for chi
+  compute_grad->forward_fields();    // also forwards v and rho for chi
   if (interface_flag) {
     // Need to save, wiped in exchange
     compute_interface->store_forces();
@@ -360,9 +354,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   // Position half-step
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-      for (a = 0; a < dim; a++) {
+      for (a = 0; a < dim; a++) { x[i][a] += dtv * v[i][a]; }
-        x[i][a] += dtv * v[i][a];
-      }
     }
   }
 
@@ -374,9 +366,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
         if (status[i] & PHASECHECK) continue;
 
         divu = 0;
-        for (a = 0; a < dim; a++) {
+        for (a = 0; a < dim; a++) { divu += gradv[i][a * (1 + dim)]; }
-          divu += gradv[i][a * (1 + dim)];
-        }
         rho[i] += dtf * (drho[i] - rho[i] * divu);
       }
     }
@@ -392,16 +382,12 @@ void FixRHEO::initial_integrate(int /*vflag*/)
       if (mask[i] & groupbit) {
         for (a = 0; a < dim; a++) {
           x[i][a] += dtv * vshift[i][a];
-          for (b = 0; b < dim; b++) {
+          for (b = 0; b < dim; b++) { v[i][a] += dtv * vshift[i][b] * gradv[i][a * dim + b]; }
-            v[i][a] += dtv * vshift[i][b] * gradv[i][a * dim + b];
-          }
         }
 
         if (!rhosum_flag) {
           if (status[i] & PHASECHECK) continue;
-          for (a = 0; a < dim; a++) {
+          for (a = 0; a < dim; a++) { rho[i] += dtv * vshift[i][a] * gradr[i][a]; }
-            rho[i] += dtv * vshift[i][a] * gradr[i][a];
-          }
         }
       }
     }
@@ -412,10 +398,9 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
 void FixRHEO::pre_force(int /*vflag*/)
 {
-  compute_kernel->compute_coordination(); // Needed for rho sum
+  compute_kernel->compute_coordination();    // Needed for rho sum
 
-  if (rhosum_flag)
+  if (rhosum_flag) compute_rhosum->compute_peratom();
-    compute_rhosum->compute_peratom();
 
   compute_kernel->compute_peratom();
 
@@ -428,22 +413,19 @@ void FixRHEO::pre_force(int /*vflag*/)
   compute_grad->compute_peratom();
   compute_grad->forward_gradients();
 
-  if (shift_flag)
+  if (shift_flag) compute_vshift->compute_peratom();
-    compute_vshift->compute_peratom();
 
   // Remove temporary options
   int *mask = atom->mask;
   int *status = atom->rheo_status;
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++)
-    if (mask[i] & groupbit)
+    if (mask[i] & groupbit) status[i] &= OPTIONSMASK;
-      status[i] &= OPTIONSMASK;
 
   // Calculate surfaces, update status
   if (surface_flag) {
     compute_surface->compute_peratom();
-    if (shift_flag)
+    if (shift_flag) compute_vshift->correct_surfaces();
-      compute_vshift->correct_surfaces();
   }
 }
 
@@ -452,8 +434,7 @@ void FixRHEO::pre_force(int /*vflag*/)
 void FixRHEO::final_integrate()
 {
   int nlocal = atom->nlocal;
-  if (igroup == atom->firstgroup)
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
-    nlocal = atom->nfirst;
 
   double dtfm, divu;
   int i, a;
@@ -482,9 +463,7 @@ void FixRHEO::final_integrate()
         dtfm = dtf / mass[type[i]];
       }
 
-      for (a = 0; a < dim; a++) {
+      for (a = 0; a < dim; a++) { v[i][a] += dtfm * f[i][a]; }
-        v[i][a] += dtfm * f[i][a];
-      }
     }
   }
 
@@ -496,9 +475,7 @@ void FixRHEO::final_integrate()
         if (status[i] & PHASECHECK) continue;
 
         divu = 0;
-        for (a = 0; a < dim; a++) {
+        for (a = 0; a < dim; a++) { divu += gradv[i][a * (1 + dim)]; }
-          divu += gradv[i][a * (1 + dim)];
-        }
         rho[i] += dtf * (drho[i] - rho[i] * divu);
       }
     }

src/RHEO/fix_rheo.h

@@ -72,11 +72,11 @@ class FixRHEO : public Fix {
 
 namespace RHEO_NS {
 
-  enum {QUINTIC, WENDLANDC4, RK0, RK1, RK2};
+  enum { QUINTIC, WENDLANDC4, RK0, RK1, RK2 };
-  enum {COORDINATION, DIVR};
+  enum { COORDINATION, DIVR };
 
   // Status variables
-  enum Status{
+  enum Status {
     // Phase status
     STATUS_SOLID = 1 << 0,
     // Gap for future phase: STATUS_ = 1 << 1,
@@ -95,12 +95,13 @@ namespace RHEO_NS {
   };
 
   // Masks and their inverses
-  #define PHASEMASK 0xFFFFFFFC     // 11111111111111111111111111111100
+  enum {
-  #define PHASECHECK 0x00000003    // 00000000000000000000000000000011
+    PHASECHECK = 0x00000003,      // 00000000000000000000000000000011
-  #define SURFACEMASK 0xFFFFFFC3   // 11111111111111111111111111000011
+    SURFACECHECK = 0x0000003C,    // 00000000000000000000000000111100
-  #define SURFACECHECK 0x0000003C  // 00000000000000000000000000111100
+    OPTIONSMASK = 0xFFFFFC3F,     // 11111111111111111111110000111111
-  #define OPTIONSMASK 0xFFFFFC3F   // 11111111111111111111110000111111
+    SURFACEMASK = 0xFFFFFFC3,     // 11111111111111111111111111000011
-
+    PHASEMASK = 0xFFFFFFFC        // 11111111111111111111111111111100
+  };
 }    // namespace RHEO_NS
 }    // namespace LAMMPS_NS
 

src/RHEO/fix_rheo_oxidation.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -28,34 +27,35 @@
 #include "fix_rheo.h"
 #include "force.h"
 #include "modify.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 #include "update.h"
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
 using namespace FixConst;
-enum {NONE, CONSTANT};
+enum { NONE, CONSTANT };
 
 static const char cite_rheo_oxide[] =
-  "@article{ApplMathModel.130.310,\n"
+    "@article{ApplMathModel.130.310,\n"
-  " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
+    " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
-  " journal = {Applied Mathematical Modelling},\n"
+    " journal = {Applied Mathematical Modelling},\n"
-  " volume = {130},\n"
+    " volume = {130},\n"
-  " pages = {310-326},\n"
+    " pages = {310-326},\n"
-  " year = {2024},\n"
+    " year = {2024},\n"
-  " issn = {0307-904X},\n"
+    " issn = {0307-904X},\n"
-  " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
+    " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
-  " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
+    " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and "
-  "}\n\n";
+    "Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
+    "}\n\n";
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), compute_surface(nullptr), fix_rheo(nullptr)
+    Fix(lmp, narg, arg), compute_surface(nullptr), fix_rheo(nullptr)
 {
-  if (narg != 6) error->all(FLERR,"Illegal fix command");
+  if (narg != 6) error->all(FLERR, "Illegal fix command");
 
   force_reneighbor = 1;
   next_reneighbor = -1;
@@ -65,7 +65,8 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) :
   if (cut <= 0.0) error->all(FLERR, "Illegal bond cutoff {} in fix rheo/oxidation", cut);
 
   btype = utils::inumeric(FLERR, arg[4], false, lmp);
-  if (btype < 1 || btype > atom->nbondtypes) error->all(FLERR, "Illegal value {} for bond type in fix rheo/oxidation", btype);
+  if (btype < 1 || btype > atom->nbondtypes)
+    error->all(FLERR, "Illegal value {} for bond type in fix rheo/oxidation", btype);
 
   rsurf = utils::numeric(FLERR, arg[5], false, lmp);
   if (rsurf <= 0.0) error->all(FLERR, "Illegal surface distance {} in fix rheo/oxidation", cut);
@@ -77,9 +78,7 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-FixRHEOOxidation::~FixRHEOOxidation()
+FixRHEOOxidation::~FixRHEOOxidation() {}
-{
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -100,19 +99,16 @@ void FixRHEOOxidation::init()
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/oxidation");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
-  if (cut > fix_rheo->cut)
+  if (cut > fix_rheo->cut) error->all(FLERR, "Bonding length exceeds kernel cutoff");
-    error->all(FLERR, "Bonding length exceeds kernel cutoff");
 
-  if (rsurf >= fix_rheo->cut)
+  if (rsurf >= fix_rheo->cut) error->all(FLERR, "Surface distance must be less than kernel cutoff");
-    error->all(FLERR, "Surface distance must be less than kernel cutoff");
 
   if (!force->bond) error->all(FLERR, "Must define a bond style with fix rheo/oxidation");
   if (!atom->avec->bonds_allow) error->all(FLERR, "Fix rheo/oxidation requires atom bonds");
 
   int tmp1, tmp2;
   index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
-  if (index_nb == -1)
+  if (index_nb == -1) error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation");
-    error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation");
   nbond = atom->ivector[index_nb];
 
   // need a half neighbor list
@@ -135,16 +131,14 @@ void FixRHEOOxidation::setup_pre_force(int /*vflag*/)
   // but enforce to be consistent with other RHEO fixes
   fix_rheo->oxidation_fix_defined = 1;
 
-  if (!fix_rheo->surface_flag) error->all(FLERR,
+  if (!fix_rheo->surface_flag)
-      "fix rheo/oxidation requires surface calculation in fix rheo");
+    error->all(FLERR, "fix rheo/oxidation requires surface calculation in fix rheo");
   compute_surface = fix_rheo->compute_surface;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixRHEOOxidation::pre_force(int /*vflag*/)
+void FixRHEOOxidation::pre_force(int /*vflag*/) {}
-{
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -235,7 +229,8 @@ void FixRHEOOxidation::post_integrate()
 
       if (!newton_bond || (tagi < tagj)) {
         if (num_bond[i] == atom->bond_per_atom)
-          error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/oxidation for atom {}", tagi);
+          error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/oxidation for atom {}",
+                     tagi);
         bond_type[i][num_bond[i]] = btype;
         bond_atom[i][num_bond[i]] = tagj;
         num_bond[i]++;
@@ -246,8 +241,7 @@ void FixRHEOOxidation::post_integrate()
   int added_bonds_all;
   MPI_Allreduce(&added_bonds, &added_bonds_all, 1, MPI_INT, MPI_SUM, world);
 
-  if (added_bonds_all > 0)
+  if (added_bonds_all > 0) next_reneighbor = update->ntimestep;
-    next_reneighbor = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -257,14 +251,13 @@ void FixRHEOOxidation::post_force(int /*vflag*/)
   int *status = atom->rheo_status;
   int *num_bond = atom->num_bond;
   for (int i = 0; i < atom->nlocal; i++)
-    if (num_bond[i] != 0)
+    if (num_bond[i] != 0) status[i] |= STATUS_NO_SHIFT;
-      status[i] |= STATUS_NO_SHIFT;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf,
+int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                        int * /*pbc*/)
 {
   int i, j, m;
   double **x = atom->x;

src/RHEO/fix_rheo_pressure.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -32,24 +31,23 @@
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
-enum {NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL};
+enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL };
 
 static constexpr double SEVENTH = 1.0 / 7.0;
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr),
+    Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr),
-  rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr),
+    rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr),
-  fix_rheo(nullptr)
+    fix_rheo(nullptr)
 {
-  if (narg < 4) error->all(FLERR,"Illegal fix command");
+  if (narg < 4) error->all(FLERR, "Illegal fix command");
 
   comm_forward = 1;
 
   // Currently can only have one instance of fix rheo/pressure
-  if (igroup != 0)
+  if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all");
-    error->all(FLERR,"fix rheo/pressure command requires group all");
 
   int i, nlo, nhi;
   int n = atom->ntypes;
@@ -66,11 +64,9 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
     if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure", error);
 
     if (strcmp(arg[iarg + 1], "linear") == 0) {
-      for (i = nlo; i <= nhi; i++)
+      for (i = nlo; i <= nhi; i++) pressure_style[i] = LINEAR;
-        pressure_style[i] = LINEAR;
     } else if (strcmp(arg[iarg + 1], "tait/water") == 0) {
-      for (i = nlo; i <= nhi; i++)
+      for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER;
-        pressure_style[i] = TAITWATER;
     } else if (strcmp(arg[iarg + 1], "tait/general") == 0) {
       if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait", error);
 
@@ -94,14 +90,14 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
         c_cubic[i] = c_cubic_one;
       }
     } else {
-      error->all(FLERR,"Illegal fix command, {}", arg[iarg]);
+      error->all(FLERR, "Illegal fix command, {}", arg[iarg]);
     }
     iarg += 2;
   }
 
   for (i = 1; i <= n; i++)
     if (pressure_style[i] == NONE)
-      error->all(FLERR,"Must specify pressure for atom type {} in fix/rheo/pressure", i);
+      error->all(FLERR, "Must specify pressure for atom type {} in fix/rheo/pressure", i);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -170,16 +166,15 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit)
+    if (mask[i] & groupbit) pressure[i] = calc_pressure(rho[i], type[i]);
-      pressure[i] = calc_pressure(rho[i], type[i]);
 
   if (comm_forward) comm->forward_comm(this);
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
+int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                       int * /*pbc*/)
 {
   double *pressure = atom->pressure;
   int m = 0;
@@ -197,9 +192,7 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
   double *pressure = atom->pressure;
   int m = 0;
   int last = first + n;
-  for (int i = first; i < last; i++) {
+  for (int i = first; i < last; i++) { pressure[i] = buf[m++]; }
-    pressure[i] = buf[m++];
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -235,7 +228,8 @@ double FixRHEOPressure::calc_rho(double p, int type)
   if (pressure_style[type] == LINEAR) {
     rho = csqinv[type] * p + rho0[type];
   } else if (pressure_style[type] == CUBIC) {
-    error->one(FLERR, "Rho calculation from pressure not yet supported for cubic pressure equation");
+    error->one(FLERR,
+               "Rho calculation from pressure not yet supported for cubic pressure equation");
   } else if (pressure_style[type] == TAITWATER) {
     rho = pow(7.0 * p + csq[type] * rho0[type], SEVENTH);
     rho *= pow(rho0[type], 6.0 * SEVENTH);

src/RHEO/fix_rheo_thermal.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -27,44 +26,45 @@
 #include "compute_rheo_vshift.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_rheo.h"
 #include "fix_bond_history.h"
+#include "fix_rheo.h"
 #include "fix_update_special_bonds.h"
 #include "force.h"
 #include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "neighbor.h"
 #include "pair.h"
 #include "update.h"
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
 using namespace FixConst;
-enum {NONE, CONSTANT};
+enum { NONE, CONSTANT };
 
 static const char cite_rheo_oxide[] =
-  "@article{ApplMathModel.130.310,\n"
+    "@article{ApplMathModel.130.310,\n"
-  " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
+    " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
-  " journal = {Applied Mathematical Modelling},\n"
+    " journal = {Applied Mathematical Modelling},\n"
-  " volume = {130},\n"
+    " volume = {130},\n"
-  " pages = {310-326},\n"
+    " pages = {310-326},\n"
-  " year = {2024},\n"
+    " year = {2024},\n"
-  " issn = {0307-904X},\n"
+    " issn = {0307-904X},\n"
-  " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
+    " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
-  " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
+    " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and "
-  "}\n\n";
+    "Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
+    "}\n\n";
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), cv(nullptr), Tc(nullptr), kappa(nullptr), L(nullptr), cv_style(nullptr),
+    Fix(lmp, narg, arg), cv(nullptr), Tc(nullptr), kappa(nullptr), L(nullptr), cv_style(nullptr),
-  Tc_style(nullptr), kappa_style(nullptr),  L_style(nullptr), list(nullptr), fix_rheo(nullptr),
+    Tc_style(nullptr), kappa_style(nullptr), L_style(nullptr), list(nullptr), fix_rheo(nullptr),
-  compute_grad(nullptr), compute_vshift(nullptr), fix_update_special_bonds(nullptr)
+    compute_grad(nullptr), compute_vshift(nullptr), fix_update_special_bonds(nullptr)
 {
-  if (narg < 4) error->all(FLERR,"Illegal fix command");
+  if (narg < 4) error->all(FLERR, "Illegal fix command");
 
   force_reneighbor = 1;
   next_reneighbor = -1;
@@ -72,8 +72,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
   comm_forward = 0;
 
   // Currently can only have one instance of fix rheo/thermal
-  if (igroup != 0)
+  if (igroup != 0) error->all(FLERR, "fix rheo/thermal command requires group all");
-    error->all(FLERR,"fix rheo/thermal command requires group all");
 
   int i, nlo, nhi;
   int n = atom->ntypes;
@@ -97,13 +96,14 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
 
   int iarg = 3;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"conductivity") == 0) {
+    if (strcmp(arg[iarg], "conductivity") == 0) {
       if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity", error);
       utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
 
       // Conductivity arguments
       if (strcmp(arg[iarg + 2], "constant") == 0) {
-        if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error);
+        if (iarg + 3 >= narg)
+          utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error);
 
         double kappa_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         if (kappa_one < 0.0) error->all(FLERR, "The conductivity must be positive");
@@ -114,7 +114,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
           kappa[i] = kappa_one;
         }
       } else {
-        error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]);
+        error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]);
       }
 
       iarg += 2;
@@ -124,7 +124,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
 
       // Cv arguments
       if (strcmp(arg[iarg + 2], "constant") == 0) {
-        if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error);
+        if (iarg + 3 >= narg)
+          utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error);
 
         double cv_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         if (cv_one < 0.0) error->all(FLERR, "The specific heat must be positive");
@@ -136,7 +137,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
         }
 
       } else {
-        error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]);
+        error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]);
       }
 
       iarg += 2;
@@ -146,7 +147,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
 
       // T freeze arguments
       if (strcmp(arg[iarg + 2], "constant") == 0) {
-        if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error);
+        if (iarg + 3 >= narg)
+          utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error);
 
         double Tc_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         iarg += 2;
@@ -167,7 +169,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
 
       // Cv arguments
       if (strcmp(arg[iarg + 2], "constant") == 0) {
-        if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error);
+        if (iarg + 3 >= narg)
+          utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error);
 
         double L_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         if (L_one < 0.0) error->all(FLERR, "The latent heat must be positive");
@@ -179,7 +182,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
         }
 
       } else {
-        error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]);
+        error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]);
       }
 
       iarg += 2;
@@ -195,18 +198,20 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
       cutsq_bond = cut_bond * cut_bond;
       iarg += 3;
     } else {
-      error->all(FLERR,"Unknown fix rheo/thermal keyword: {}", arg[iarg]);
+      error->all(FLERR, "Unknown fix rheo/thermal keyword: {}", arg[iarg]);
     }
   }
 
-
   for (i = 1; i <= n; i++) {
     if (cv_style[i] == NONE)
       error->all(FLERR, "Must specify specific/heat for atom type {} in fix/rheo/thermal", i);
     if (kappa_style[i] == NONE)
       error->all(FLERR, "Must specify conductivity for atom type {} in fix/rheo/thermal", i);
     if (Tc_style[i] == NONE && L_style[i] != NONE)
-      error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i);
+      error->all(FLERR,
+                 "Must specify critical temperature for atom type {} to use latent heat in fix "
+                 "rheo/thermal",
+                 i);
   }
 
   if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide);
@@ -253,11 +258,9 @@ void FixRHEOThermal::init()
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
   cut_kernel = fix_rheo->cut;
 
-  if (cut_bond > cut_kernel)
+  if (cut_bond > cut_kernel) error->all(FLERR, "Bonding length exceeds kernel cutoff");
-    error->all(FLERR, "Bonding length exceeds kernel cutoff");
 
-  if (!fix_rheo->thermal_flag)
+  if (!fix_rheo->thermal_flag) error->all(FLERR, "Need to define thermal setting in fix rheo");
-    error->all(FLERR, "Need to define thermal setting in fix rheo");
   compute_grad = fix_rheo->compute_grad;
   compute_vshift = fix_rheo->compute_vshift;
 
@@ -274,19 +277,25 @@ void FixRHEOThermal::init()
     error->all(FLERR, "fix rheo/thermal command requires atom property conductivity");
 
   if (cut_bond > 0.0) {
-    if (!force->bond) error->all(FLERR, "Must define a bond style to use reactive bond generation with fix rheo/thermal");
+    if (!force->bond)
-    if (!atom->avec->bonds_allow) error->all(FLERR, "Reactive bond generation in fix rheo/thermal requires atom bonds");
+      error->all(FLERR,
+                 "Must define a bond style to use reactive bond generation with fix rheo/thermal");
+    if (!atom->avec->bonds_allow)
+      error->all(FLERR, "Reactive bond generation in fix rheo/thermal requires atom bonds");
 
     // all special weights must be 1.0 (no special neighbors) or there must be an instance of fix update/special/bonds
-    if (force->special_lj[0] != 1.0 || force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) {
+    if (force->special_lj[0] != 1.0 || force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 ||
+        force->special_lj[3] != 1.0) {
       auto fixes = modify->get_fix_by_style("UPDATE_SPECIAL_BONDS");
-      if (fixes.size() == 0) error->all(FLERR, "Without fix update/special/bonds, reactive bond generation in fix rheo/thermal requires special weights of 1.0");
+      if (fixes.size() == 0)
+        error->all(FLERR,
+                   "Without fix update/special/bonds, reactive bond generation in fix rheo/thermal "
+                   "requires special weights of 1.0");
       fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(fixes[0]);
     }
 
     // must have newton off so both processors will search nlist to build bonds
-    if (force->newton_pair)
+    if (force->newton_pair) error->all(FLERR, "Need Newton off for reactive bond generation");
-      error->all(FLERR, "Need Newton off for reactive bond generation");
 
     // need a half neighbor list, built only when particles freeze
     auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
@@ -334,13 +343,11 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
   int nlocal = atom->nlocal;
   int dim = domain->dimension;
 
-  if (igroup == atom->firstgroup)
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
-    nlocal = atom->nfirst;
 
   for (i = 0; i < nlocal; i++) {
     if (status[i] & STATUS_NO_SHIFT) continue;
-    for (a = 0; a < dim; a++)
+    for (a = 0; a < dim; a++) energy[i] += dt * vshift[i][a] * grade[i][a];
-      energy[i] += dt * vshift[i][a] * grade[i][a];
   }
 }
 
@@ -365,16 +372,16 @@ void FixRHEOThermal::post_integrate()
     if (status[i] & STATUS_NO_INTEGRATION) continue;
 
     itype = type[i];
-    cvi = calc_cv(i, itype);
+    cvi = calc_cv(itype);
     energy[i] += dth * heatflow[i];
     temperature[i] = energy[i] / cvi;
 
     if (Tc_style[itype] != NONE) {
       Ti = temperature[i];
-      Tci = calc_Tc(i, itype);
+      Tci = calc_Tc(itype);
 
       if (L_style[itype] != NONE) {
-        Li = calc_L(i, itype);
+        Li = calc_L(itype);
         if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi);
         temperature[i] = Ti;
       }
@@ -448,8 +455,7 @@ void FixRHEOThermal::post_neighbor()
 
   for (i = 0; i < nall; i++) {
     itype = type[i];
-    if (kappa_style[itype] == CONSTANT)
+    if (kappa_style[itype] == CONSTANT) conductivity[i] = kappa[itype];
-      conductivity[i] = kappa[itype];
   }
 }
 
@@ -470,15 +476,15 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
   // Calculate temperature
   for (int i = 0; i < nall; i++) {
     int itype = type[i];
-    cvi = calc_cv(i, itype);
+    cvi = calc_cv(itype);
     temperature[i] = energy[i] / cvi;
 
     if (Tc_style[itype] != NONE) {
       Ti = temperature[i];
-      Tci = calc_Tc(i, itype);
+      Tci = calc_Tc(itype);
 
       if (L_style[itype] != NONE) {
-        Li = calc_L(i, itype);
+        Li = calc_L(itype);
         if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi);
         temperature[i] = Ti;
       }
@@ -539,14 +545,14 @@ void FixRHEOThermal::break_bonds()
       nmax = num_bond[i] - 1;
       if (m == nmax) {
         if (n_histories > 0)
-          for (auto &ihistory: histories)
+          for (auto &ihistory : histories)
             dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i, m);
       } else {
         bond_type[i][m] = bond_type[i][nmax];
         bond_atom[i][m] = bond_atom[i][nmax];
         if (n_histories > 0) {
-          for (auto &ihistory: histories) {
+          for (auto &ihistory : histories) {
-            auto fix_bond_history = dynamic_cast<FixBondHistory *>  (ihistory);
+            auto fix_bond_history = dynamic_cast<FixBondHistory *>(ihistory);
             fix_bond_history->shift_history(i, m, nmax);
             fix_bond_history->delete_history(i, nmax);
           }
@@ -559,9 +565,7 @@ void FixRHEOThermal::break_bonds()
       //  only update for atom with lower tag
       if (fix_update_special_bonds) {
         if ((i < nlocal) && (j < nlocal) && melti && meltj) {
-          if (tag[i] < tag[j]) {
+          if (tag[i] < tag[j]) { fix_update_special_bonds->add_broken_bond(i, j); }
-            fix_update_special_bonds->add_broken_bond(i, j);
-          }
         } else {
           fix_update_special_bonds->add_broken_bond(i, j);
         }
@@ -592,8 +596,8 @@ void FixRHEOThermal::break_bonds()
           bond_type[i][m] = bond_type[i][nmax];
           bond_atom[i][m] = bond_atom[i][nmax];
           if (n_histories > 0)
-            for (auto &ihistory: histories) {
+            for (auto &ihistory : histories) {
-              auto fix_bond_history = dynamic_cast<FixBondHistory *> (ihistory);
+              auto fix_bond_history = dynamic_cast<FixBondHistory *>(ihistory);
               fix_bond_history->shift_history(i, m, nmax);
               fix_bond_history->delete_history(i, nmax);
             }
@@ -611,8 +615,8 @@ void FixRHEOThermal::break_bonds()
           bond_type[j][m] = bond_type[j][nmax];
           bond_atom[j][m] = bond_atom[j][nmax];
           if (n_histories > 0)
-            for (auto &ihistory: histories) {
+            for (auto &ihistory : histories) {
-              auto fix_bond_history = dynamic_cast<FixBondHistory *>  (ihistory);
+              auto fix_bond_history = dynamic_cast<FixBondHistory *>(ihistory);
               fix_bond_history->shift_history(j, m, nmax);
               fix_bond_history->delete_history(j, nmax);
             }
@@ -691,7 +695,7 @@ void FixRHEOThermal::create_bonds()
       // If newton bond off, add to both, otherwise add to whichever has a smaller tag
       if ((i < nlocal) && (!newton_bond || (tag[i] < tag[j]))) {
         if (num_bond[i] == atom->bond_per_atom)
-          error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal");
+          error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/thermal");
         bond_type[i][num_bond[i]] = btype;
         bond_atom[i][num_bond[i]] = tag[j];
         num_bond[i]++;
@@ -699,7 +703,7 @@ void FixRHEOThermal::create_bonds()
 
       if ((j < nlocal) && (!newton_bond || (tag[j] < tag[i]))) {
         if (num_bond[j] == atom->bond_per_atom)
-          error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal");
+          error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/thermal");
         bond_type[j][num_bond[j]] = btype;
         bond_atom[j][num_bond[j]] = tag[i];
         num_bond[j]++;
@@ -712,41 +716,35 @@ void FixRHEOThermal::create_bonds()
 
 /* ---------------------------------------------------------------------- */
 
-double FixRHEOThermal::calc_cv(int i, int itype)
+double FixRHEOThermal::calc_cv(int itype)
 {
-  if (cv_style[itype] == CONSTANT) {
+  if (cv_style[itype] == CONSTANT) { return cv[itype]; }
-    return cv[itype];
-  }
 
   return 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-double FixRHEOThermal::calc_Tc(int i, int itype)
+double FixRHEOThermal::calc_Tc(int itype)
 {
-  if (Tc_style[itype] == CONSTANT) {
+  if (Tc_style[itype] == CONSTANT) { return Tc[itype]; }
-    return Tc[itype];
-  }
 
   return 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-double FixRHEOThermal::calc_L(int i, int itype)
+double FixRHEOThermal::calc_L(int itype)
 {
-  if (L_style[itype] == CONSTANT) {
+  if (L_style[itype] == CONSTANT) { return L[itype]; }
-    return L[itype];
-  }
 
   return 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf,
+int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                      int * /*pbc*/)
 {
   int *status = atom->rheo_status;
   double **x = atom->x;

src/RHEO/fix_rheo_thermal.h

@@ -42,9 +42,9 @@ class FixRHEOThermal : public Fix {
   int pack_forward_comm(int, int *, double *, int, int *) override;
   void unpack_forward_comm(int, int, double *) override;
   void reset_dt() override;
-  double calc_cv(int, int);
+  double calc_cv(int);
-  double calc_Tc(int, int);
+  double calc_Tc(int);
-  double calc_L(int, int);
+  double calc_L(int);
 
  private:
   double *cv, *Tc, *kappa, *L;

src/RHEO/fix_rheo_viscosity.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -33,13 +32,13 @@
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
-enum {NONE, CONSTANT, POWER};
+enum { NONE, CONSTANT, POWER };
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), eta(nullptr), npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr),
+    Fix(lmp, narg, arg), eta(nullptr), npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr),
-  viscosity_style(nullptr), fix_rheo(nullptr), compute_grad(nullptr)
+    viscosity_style(nullptr), fix_rheo(nullptr), compute_grad(nullptr)
 {
   if (narg < 4) error->all(FLERR, "Illegal fix command");
 
@@ -48,8 +47,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
   evolve_flag = 0;
 
   // Currently can only have one instance of fix rheo/viscosity
-  if (igroup != 0)
+  if (igroup != 0) error->all(FLERR, "fix rheo/viscosity command requires group all");
-    error->all(FLERR,"fix rheo/viscosity command requires group all");
 
   int i, nlo, nhi;
   int n = atom->ntypes;
@@ -107,7 +105,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
 
   for (i = 1; i <= n; i++)
     if (viscosity_style[i] == NONE)
-      error->all(FLERR,"Must specify viscosity for atom type {} in fix/rheo/viscosity", i);
+      error->all(FLERR, "Must specify viscosity for atom type {} in fix/rheo/viscosity", i);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -172,8 +170,7 @@ void FixRHEOViscosity::post_neighbor()
 
   for (i = 0; i < nall; i++) {
     itype = type[i];
-    if (viscosity_style[itype])
+    if (viscosity_style[itype]) viscosity[i] = eta[itype];
-      viscosity[i] = eta[itype];
   }
 }
 
@@ -195,7 +192,6 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
   int nlocal = atom->nlocal;
   int dim = domain->dimension;
 
-
   for (i = 0; i < nlocal; i++) {
     itype = type[i];
     if (viscosity_style[itype] == POWER) {
@@ -222,8 +218,8 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
 
 /* ---------------------------------------------------------------------- */
 
-int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
+int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                        int * /*pbc*/)
 {
   double *viscosity = atom->viscosity;
   int m = 0;
@@ -241,7 +237,5 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
   double *viscosity = atom->viscosity;
   int m = 0;
   int last = first + n;
-  for (int i = first; i < last; i++) {
+  for (int i = first; i < last; i++) { viscosity[i] = buf[m++]; }
-    viscosity[i] = buf[m++];
-  }
 }

src/RHEO/pair_rheo.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -21,9 +20,9 @@
 
 #include "atom.h"
 #include "comm.h"
-#include "compute_rheo_kernel.h"
 #include "compute_rheo_grad.h"
 #include "compute_rheo_interface.h"
+#include "compute_rheo_kernel.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
@@ -32,8 +31,8 @@
 #include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
-#include "neighbor.h"
 #include "neigh_list.h"
+#include "neighbor.h"
 #include "utils.h"
 
 #include <cmath>
@@ -47,9 +46,8 @@ static constexpr double EPSILON = 1e-2;
 /* ---------------------------------------------------------------------- */
 
 PairRHEO::PairRHEO(LAMMPS *lmp) :
-  Pair(lmp), csq(nullptr), rho0(nullptr), cs(nullptr), compute_kernel(nullptr),
+    Pair(lmp), csq(nullptr), rho0(nullptr), cs(nullptr), compute_kernel(nullptr),
-  compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr),
+    compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr), fix_pressure(nullptr)
-  fix_pressure(nullptr)
 {
   restartinfo = 0;
   single_enable = 0;
@@ -82,7 +80,8 @@ void PairRHEO::compute(int eflag, int vflag)
   int pair_force_flag, pair_rho_flag, pair_avisc_flag;
   int fluidi, fluidj;
   double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave;
-  double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, kappa_ave, dT_prefactor;
+  double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave,
+      kappa_ave, dT_prefactor;
   double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
   double *dWij, *dWji;
   double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3];
@@ -150,8 +149,10 @@ void PairRHEO::compute(int eflag, int vflag)
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    if (rmass) imass = rmass[i];
+    if (rmass)
-    else imass = mass[itype];
+      imass = rmass[i];
+    else
+      imass = mass[itype];
     etai = viscosity[i];
     fluidi = !(status[i] & PHASECHECK);
     if (thermal_flag) {
@@ -173,8 +174,10 @@ void PairRHEO::compute(int eflag, int vflag)
         r = sqrt(rsq);
         rinv = 1 / r;
 
-        if (rmass) jmass = rmass[i];
+        if (rmass)
-        else jmass = mass[jtype];
+          jmass = rmass[i];
+        else
+          jmass = mass[jtype];
         etaj = viscosity[j];
         fluidj = !(status[j] & PHASECHECK);
         if (thermal_flag) {
@@ -217,7 +220,7 @@ void PairRHEO::compute(int eflag, int vflag)
         if (interface_flag) {
           if (fluidi && (!fluidj)) {
             compute_interface->correct_v(vj, vi, j, i);
-            rhoj = compute_interface->correct_rho(j, i);
+            rhoj = compute_interface->correct_rho(j);
             Pj = fix_pressure->calc_pressure(rhoj, jtype);
 
             if ((chi[j] > 0.9) && (r < (cutk * 0.5)))
@@ -225,7 +228,7 @@ void PairRHEO::compute(int eflag, int vflag)
 
           } else if ((!fluidi) && fluidj) {
             compute_interface->correct_v(vi, vj, i, j);
-            rhoi = compute_interface->correct_rho(i, j);
+            rhoi = compute_interface->correct_rho(i);
             Pi = fix_pressure->calc_pressure(rhoi, itype);
 
             if (chi[i] > 0.9 && r < (cutk * 0.5))
@@ -251,7 +254,8 @@ void PairRHEO::compute(int eflag, int vflag)
           } else {
             kappa_ave = 0.5 * (kappai + kappaj);
           }
-          dT_prefactor = 2.0 * kappa_ave * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj);
+          dT_prefactor =
+              2.0 * kappa_ave * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj);
 
           dT = dot3(dx, dWij);
           heatflow[i] += dT * dT_prefactor;
@@ -295,8 +299,7 @@ void PairRHEO::compute(int eflag, int vflag)
           // Now compute viscous eta*Lap[v] terms
           for (a = 0; a < dim; a++) {
             fv[a] = 0.0;
-            for (b = 0; b < dim; b++)
+            for (b = 0; b < dim; b++) fv[a] += dv[a] * dx[b] * dWij[b];
-              fv[a] += dv[a] * dx[b] * dWij[b];
             fv[a] *= 2.0 * eta_ave * voli * volj * rinv * rinv;
           }
 
@@ -309,13 +312,13 @@ void PairRHEO::compute(int eflag, int vflag)
 
           // Note the virial's definition is hazy, e.g. viscous contributions will depend on rotation
           if (evflag)
-            ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]);
+            ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1],
+                         dx[2]);
 
           if (newton_pair || j < nlocal) {
-            for (a = 0; a < dim; a ++) {
+            for (a = 0; a < dim; a++) {
               fv[a] = 0.0;
-              for (b = 0; b < dim; b++)
+              for (b = 0; b < dim; b++) fv[a] += (vi[a] - vj[a]) * dx[b] * dWji[b];
-                fv[a] += (vi[a] - vj[a]) * dx[b] * dWji[b];
               fv[a] *= -2.0 * eta_ave * voli * volj * rinv * rinv;
               // flip sign here b/c -= at accummulator
             }
@@ -329,7 +332,8 @@ void PairRHEO::compute(int eflag, int vflag)
             f[j][2] -= ft[2];
 
             if (evflag)
-              ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], -dx[0], -dx[1], -dx[2]);
+              ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], -dx[0], -dx[1],
+                           -dx[2]);
           }
 
           if (compute_interface) {
@@ -352,8 +356,7 @@ void PairRHEO::compute(int eflag, int vflag)
           if (laplacian_order >= 1) {
             psi_ij = rhoj - rhoi;
             Fij = -rinv * rinv * dot3(dx, dWij);
-            for (a = 0; a < dim; a++)
+            for (a = 0; a < dim; a++) psi_ij += 0.5 * (gradr[i][a] + gradr[j][a]) * dx[a];
-              psi_ij += 0.5 * (gradr[i][a] + gradr[j][a]) * dx[a];
             drho[i] += 2 * rho_damp * psi_ij * Fij * volj;
           } else {
             drho_damp = 2 * rho_damp * ((rhoj - rho0[jtype]) - (rhoi - rho0[itype])) * rinv * wp;
@@ -393,8 +396,7 @@ void PairRHEO::allocate()
 
   memory->create(setflag, n + 1, n + 1, "pair:setflag");
   for (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++)
+    for (int j = i; j <= n; j++) setflag[i][j] = 0;
-      setflag[i][j] = 0;
 
   memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
 }
@@ -423,7 +425,8 @@ void PairRHEO::settings(int narg, char **arg)
       iarg++;
     } else if (strcmp(arg[iarg], "harmonic/means") == 0) {
       harmonic_means_flag = 1;
-    } else error->all(FLERR, "Illegal pair_style command, {}", arg[iarg]);
+    } else
+      error->all(FLERR, "Illegal pair_style command, {}", arg[iarg]);
     iarg++;
   }
 }
@@ -434,10 +437,8 @@ void PairRHEO::settings(int narg, char **arg)
 
 void PairRHEO::coeff(int narg, char **arg)
 {
-  if (narg != 2)
+  if (narg != 2) error->all(FLERR, "Incorrect number of args for pair_style rheo coefficients");
-    error->all(FLERR, "Incorrect number of args for pair_style rheo coefficients");
+  if (!allocated) allocate();
-  if (!allocated)
-    allocate();
 
   int ilo, ihi, jlo, jhi;
   utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
@@ -451,8 +452,7 @@ void PairRHEO::coeff(int narg, char **arg)
     }
   }
 
-  if (count == 0)
+  if (count == 0) error->all(FLERR, "Incorrect args for pair rheo coefficients");
-    error->all(FLERR, "Incorrect args for pair rheo coefficients");
 }
 
 /* ----------------------------------------------------------------------
@@ -481,7 +481,8 @@ void PairRHEO::setup()
   rho0 = fix_rheo->rho0;
 
   if (cutk != fix_rheo->cut)
-    error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, fix_rheo->cut);
+    error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk,
+               fix_rheo->cut);
 
   cutksq = cutk * cutk;
   cutkinv = 1.0 / cutk;
@@ -490,11 +491,9 @@ void PairRHEO::setup()
 
   int n = atom->ntypes;
   memory->create(cs, n + 1, "rheo:cs");
-  for (int i = 1; i <= n; i++)
+  for (int i = 1; i <= n; i++) cs[i] = sqrt(csq[i]);
-    cs[i] = sqrt(csq[i]);
 
-  if (comm->ghost_velocity == 0)
+  if (comm->ghost_velocity == 0) error->all(FLERR, "Pair RHEO requires ghost atoms store velocity");
-    error->all(FLERR, "Pair RHEO requires ghost atoms store velocity");
 
   if (laplacian_order == -1) {
     if (fix_rheo->kernel_style == RK2)
@@ -512,8 +511,7 @@ void PairRHEO::setup()
 
 double PairRHEO::init_one(int i, int j)
 {
-  if (setflag[i][j] == 0)
+  if (setflag[i][j] == 0) error->all(FLERR, "All pair rheo coeffs are not set");
-    error->all(FLERR, "All pair rheo coeffs are not set");
 
   return cutk;
 }

src/RHEO/pair_rheo.h

@@ -37,7 +37,7 @@ class PairRHEO : public Pair {
   void unpack_reverse_comm(int, int *, double *) override;
 
  protected:
-  double cutk, *csq, *rho0;        // From fix RHEO
+  double cutk, *csq, *rho0;    // From fix RHEO
   double *cs, cutksq, cutkinv, cutkinv3, av, rho_damp;
 
   int laplacian_order;

src/RHEO/pair_rheo_solid.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -221,15 +220,13 @@ void PairRHEOSolid::coeff(int narg, char **arg)
 void PairRHEOSolid::init_style()
 {
   if (comm->ghost_velocity == 0)
-    error->all(FLERR,"Pair rheo/solid requires ghost atoms store velocity");
+    error->all(FLERR, "Pair rheo/solid requires ghost atoms store velocity");
 
-  if (!atom->rheo_status_flag)
+  if (!atom->rheo_status_flag) error->all(FLERR, "Pair rheo/solid requires atom_style rheo");
-    error->all(FLERR,"Pair rheo/solid requires atom_style rheo");
 
   neighbor->add_request(this);
 }
 
-
 /* ----------------------------------------------------------------------
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */