Merge pull request #4633 from jtclemm/rheo

Fixes to RHEO package

doc/src/fix_deform.rst

@@ -298,7 +298,8 @@ Note that it should return the "change" in box length, not the
 absolute box length.  This means it should evaluate to 0.0 when
 invoked on the initial timestep of the run following the definition of
 fix deform.  It should evaluate to a value > 0.0 to dilate the box at
-future times, or a value < 0.0 to compress the box.
+future times, or a value < 0.0 to compress the box. The exception
+would be if the run command uses the *pre no* option.
 
 The variable *name2* must also be an :doc:`equal-style variable
 <variable>` and should calculate the rate of box length change, in

src/RHEO/bond_rheo_shell.cpp

@@ -147,7 +147,8 @@ void BondRHEOShell::store_data()
 
       // map to find index n
       j = atom->map(atom->bond_atom[i][m]);
-      if (j == -1) error->one(FLERR, "Atom missing in BPM bond");
+      if (j == -1) error->one(FLERR, Error::NOLASTLINE,
+                              "Atom {} missing in BPM bond", atom->bond_atom[i][m]);
 
       fix_bond_history->update_atom_value(i, m, 0, 0.0);
       fix_bond_history->update_atom_value(i, m, 1, 0.0);
@@ -217,12 +218,6 @@ void BondRHEOShell::compute(int eflag, int vflag)
       i2 = itmp;
     }
 
-    delx = x[i1][0] - x[i2][0];
-    dely = x[i1][1] - x[i2][1];
-    delz = x[i1][2] - x[i2][2];
-    rsq = delx * delx + dely * dely + delz * delz;
-    r = sqrt(rsq);
-
     // If bond hasn't been set - zero data
     if (t < EPSILON || std::isnan(t)) t = store_bond(n, i1, i2);
 
@@ -358,19 +353,20 @@ void BondRHEOShell::init_style()
   BondBPM::init_style();
 
   if (comm->ghost_velocity == 0)
-    error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity");
+    error->all(FLERR, Error::NOLASTLINE, "Bond rheo/shell requires ghost atoms store velocity");
 
   auto fixes = modify->get_fix_by_style("^rheo$");
-  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell");
+  if (fixes.size() == 0)
+    error->all(FLERR, Error::NOLASTLINE, "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, "Bond rheo/shell requires surface calculation in fix rheo");
+    error->all(FLERR, Error::NOLASTLINE, "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");
+    error->all(FLERR, Error::NOLASTLINE, "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;
@@ -389,14 +385,13 @@ void BondRHEOShell::settings(int narg, char **arg)
     if (strcmp(arg[iarg], "t/form") == 0) {
       if (iarg + 1 > narg) utils::missing_cmd_args(FLERR, "bond rheo/shell t/form", error);
       tform = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      if (tform < 0.0)
+        error->all(FLERR, iarg + 1, "Illegal bond rheo/shell command, must specify positive formation time");
       i += 1;
     } else {
-      error->all(FLERR, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]);
+      error->all(FLERR, iarg, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]);
     }
   }
-
-  if (tform < 0.0)
-    error->all(FLERR, "Illegal bond rheo/shell command, must specify positive formation time");
 }
 
 /* ----------------------------------------------------------------------
@@ -472,7 +467,7 @@ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf)
   m = 0;
   last = first + n;
 
-  for (i = first; i < last; i++) { buf[m++] = dbond[i]; }
+  for (i = first; i < last; i++) { buf[m++] = ubuf(dbond[i]).d; }
   return m;
 }
 
@@ -485,7 +480,7 @@ void BondRHEOShell::unpack_reverse_comm(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    dbond[j] += buf[m++];
+    dbond[j] += (int) ubuf(buf[m++]).i;
   }
 }
 
@@ -495,14 +490,20 @@ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce)
 {
   if (type <= 0) return 0.0;
 
-  double r0, t;
+  double r0 = -1;
+  double t = -1;
+  const auto *tag = atom->tag;
   for (int n = 0; n < atom->num_bond[i]; n++) {
-    if (atom->bond_atom[i][n] == atom->tag[j]) {
+    if (atom->bond_atom[i][n] == tag[j]) {
       r0 = fix_bond_history->get_atom_value(i, n, 0);
       t = fix_bond_history->get_atom_value(i, n, 1);
     }
   }
 
+  if (r0 == -1)
+    error->one(FLERR, Error::NOLASTLINE,
+               "Could not find bond between atoms {} and {}", tag[i], tag[j]);
+
   svector[1] = t;
   if (t < tform) return 0.0;
 

src/RHEO/compute_rheo_kernel.cpp

@@ -251,13 +251,14 @@ 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) { w = 0.0; }
-
   if (s <= 3.0) {
     tmp3 = 3.0 - s;
     tmp3sq = tmp3 * tmp3;
     w = tmp3sq * tmp3sq * tmp3;
+  } else {
+    w = 0.0;
   }
+
   if (s <= 2.0) {
     tmp2 = 2.0 - s;
     tmp2sq = tmp2 * tmp2;
@@ -285,12 +286,14 @@ double ComputeRHEOKernel::calc_dw_scalar_quintic(double r)
 
   s = r * 3.0 * cutinv;
 
-  if (s > 3.0) { wp = 0.0; }
   if (s <= 3.0) {
     tmp3 = 3.0 - s;
     tmp3sq = tmp3 * tmp3;
     wp = -5.0 * tmp3sq * tmp3sq;
+  } else {
+    wp = 0.0;
   }
+
   if (s <= 2.0) {
     tmp2 = 2.0 - s;
     tmp2sq = tmp2 * tmp2;
@@ -855,7 +858,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb
   for (int i = 0; i < n; i++) {
     int j = list[i];
     if (comm_stage == 0) {
-      buf[m++] = coordination[j];
+      buf[m++] = ubuf(coordination[j]).d;
     } else {
       if (kernel_style == RK0) {
         buf[m++] = C0[j];
@@ -878,7 +881,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
 
   for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
-      coordination[i] = buf[m++];
+      coordination[i] = (int) ubuf(buf[m++]).i;
     } else {
       if (kernel_style == RK0) {
         C0[i] = buf[m++];

src/RHEO/compute_rheo_vshift.cpp

@@ -271,6 +271,8 @@ void ComputeRHEOVShift::correct_surfaces()
   int dim = domain->dimension;
 
   double nx, ny, nz, vx, vy, vz, dot;
+  nz = 0.0;
+  vz = 0.0;
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
 

src/RHEO/fix_rheo.cpp

@@ -110,7 +110,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     kernel_style = RK2;
   } else
     error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]);
-  zmin_kernel = utils::numeric(FLERR, arg[5], false, lmp);
+  zmin_kernel = utils::inumeric(FLERR, arg[5], false, lmp);
 
   int iarg = 6;
   while (iarg < narg) {
@@ -342,8 +342,6 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   double *rmass = atom->rmass;
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
-  double **vshift;
-  if (shift_flag) vshift = compute_vshift->vshift;
 
   int nlocal = atom->nlocal;
   int rmass_flag = atom->rmass_flag;
@@ -400,6 +398,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   // Shifting atoms
   if (shift_flag) {
+    double **vshift = compute_vshift->vshift;
     for (i = 0; i < nlocal; i++) {
       if (status[i] & STATUS_NO_SHIFT) continue;
       if (status[i] & PHASECHECK) continue;

src/RHEO/fix_rheo_thermal.cpp

@@ -189,7 +189,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg], "react") == 0) {
       if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal react", error);
       cut_bond = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      btype = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+      btype = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
       comm_forward = 4;
       if (cut_bond <= 0.0) error->all(FLERR, "Illegal max bond length must be greater than zero");
       if ((btype < 1) || (btype > atom->nbondtypes))

src/RHEO/pair_rheo.cpp

@@ -80,7 +80,7 @@ void PairRHEO::compute(int eflag, int vflag)
   int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
   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 xtmp, ytmp, ztmp, wp, Ti, Tj, dT, cs_ave;
   double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj;
   double eta_ave, kappa_ave, dT_prefactor;
   double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
@@ -189,11 +189,10 @@ void PairRHEO::compute(int eflag, int vflag)
         if (!variable_csq) {
           cs_ave = 0.5 * (cs[itype] + cs[jtype]);
         } else {
-          csqi = fix_pressure->calc_csq(rhoi, i);
+          csqi = fix_pressure->calc_csq(rho[i], i);
-          csqj = fix_pressure->calc_csq(rhoj, j);
+          csqj = fix_pressure->calc_csq(rho[j], j);
           cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
         }
-        csq_ave = cs_ave * cs_ave;
 
         pair_rho_flag = 0;
         pair_force_flag = 0;
@@ -231,7 +230,7 @@ void PairRHEO::compute(int eflag, int vflag)
             Pj = fix_pressure->calc_pressure(rhoj, j);
 
             if ((chi[j] > 0.9) && (r < (cutk * 0.5)))
-              fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv;
+              fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * cs_ave * cs_ave * cutk * rinv;
 
           } else if ((!fluidi) && fluidj) {
             compute_interface->correct_v(vi, vj, i, j);
@@ -239,7 +238,7 @@ void PairRHEO::compute(int eflag, int vflag)
             Pi = fix_pressure->calc_pressure(rhoi, i);
 
             if (chi[i] > 0.9 && r < (cutk * 0.5))
-              fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv;
+              fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * cs_ave * cs_ave * cutk * rinv;
 
           } else if ((!fluidi) && (!fluidj)) {
             rhoi = rho0i;
@@ -251,7 +250,6 @@ void PairRHEO::compute(int eflag, int vflag)
             csqi = fix_pressure->calc_csq(rhoi, i);
             csqj = fix_pressure->calc_csq(rhoj, j);
             cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
-            csq_ave = cs_ave * cs_ave;
           }
         }
 

src/set.cpp

@@ -48,7 +48,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
 
 enum{ANGLE,ANGMOM,APIP_LAMBDA,BOND,CC,CHARGE,DENSITY,DIAMETER,DIHEDRAL,DIPOLE,
   DIPOLE_RANDOM,DPD_THETA,EDPD_CV,EDPD_TEMP,EPSILON,IMAGE,IMPROPER,LENGTH,
-  MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,SHAPE,
+  MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,RHEO_STATUS,SHAPE,
   SMD_CONTACT_RADIUS,SMD_MASS_DENSITY,SPH_CV,SPH_E,SPH_RHO,
   SPIN_ATOM,SPIN_ATOM_RANDOM,SPIN_ELECTRON,TEMPERATURE,THETA,THETA_RANDOM,
   TRI,TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,VOLUME,VX,VY,VZ,X,Y,Z,
@@ -296,6 +296,14 @@ void Set::process_args(int caller_flag, int narg, char **arg)
       action->keyword = RADIUS_ELECTRON;
       process_radius_election(iarg,narg,arg,action);
       invoke_choice[naction++] = &Set::invoke_radius_election;
+    } else if (strcmp(arg[iarg],"rheo/rho") == 0) {
+      action->keyword = SPH_RHO;
+      process_sph_rho(iarg,narg,arg,action);
+      invoke_choice[naction++] = &Set::invoke_sph_rho;
+    } else if (strcmp(arg[iarg],"rheo/status") == 0) {
+      action->keyword = RHEO_STATUS;
+      process_rheo_status(iarg,narg,arg,action);
+      invoke_choice[naction++] = &Set::invoke_rheo_status;
     } else if (strcmp(arg[iarg],"shape") == 0) {
       action->keyword = SHAPE;
       process_shape(iarg,narg,arg,action);
@@ -2029,6 +2037,44 @@ void Set::invoke_radius_election(Action *action)
 
 /* ---------------------------------------------------------------------- */
 
+void Set::process_rheo_status(int &iarg, int narg, char **arg, Action *action)
+{
+  if (!atom->rheo_status_flag)
+    error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
+  if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error);
+
+  if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action);
+  else {
+    action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+    if (action->ivalue1 != 0 && action->ivalue1 != 1)
+      error->one(FLERR,"Invalid rheo/status {} in set command", action->ivalue1);
+  }
+
+  iarg += 2;
+}
+
+void Set::invoke_rheo_status(Action *action)
+{
+  int nlocal = atom->nlocal;
+  int *status = atom->rheo_status;
+
+  int varflag = action->varflag;
+  int rheo_status;
+  if (!action->varflag1) rheo_status = action->ivalue1;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!select[i]) continue;
+    if (varflag) {
+      rheo_status = static_cast<int> (vec1[i]);
+      if (rheo_status != 0 && rheo_status != 1)
+        error->one(FLERR,"Invalid rheo/status in set command");
+    }
+    status[i] = rheo_status;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Set::process_shape(int &iarg, int narg, char **arg, Action *action)
 {
   if (!atom->ellipsoid_flag)

src/set.h

@@ -111,6 +111,7 @@ class Set : public Command {
   void process_quat(int &, int, char **, Action *);
   void process_quat_random(int &, int, char **, Action *);
   void process_radius_election(int &, int, char **, Action *);
+  void process_rheo_status(int &, int, char **, Action *);
   void process_shape(int &, int, char **, Action *);
   void process_smd_contact_radius(int &, int, char **, Action *);
   void process_smd_mass_density(int &, int, char **, Action *);
@@ -164,6 +165,7 @@ class Set : public Command {
   void invoke_quat(Action *);
   void invoke_quat_random(Action *);
   void invoke_radius_election(Action *);
+  void invoke_rheo_status(Action *);
   void invoke_shape(Action *);
   void invoke_smd_contact_radius(Action *);
   void invoke_smd_mass_density(Action *);