Debugging some compute property items, allow surface shifting inward

src/RHEO/compute_rheo_property_atom.cpp

@@ -70,6 +70,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
 
     if (strcmp(arg[iarg],"phase") == 0) {
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase;
+    } else if (strcmp(arg[iarg],"rho") == 0) {
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_rho;
     } else if (strcmp(arg[iarg],"chi") == 0) {
       interface_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi;
@@ -82,7 +84,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
     } else if (strcmp(arg[iarg],"surface/divr") == 0) {
       surface_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
-    } else if (strcmp(arg[iarg],"^surface/n") == 0) {
+    } else if (utils::strmatch(arg[iarg], "^surface/n")) {
       surface_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n;
       col_index[i] = get_vector_index(arg[iarg]);
@@ -91,11 +93,11 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
     } else if (strcmp(arg[iarg],"cv") == 0) {
       thermal_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv;
-    } else if (strcmp(arg[iarg],"^shift/v") == 0) {
+    } else if (utils::strmatch(arg[iarg], "^shift/v")) {
       shift_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_v;
       col_index[i] = get_vector_index(arg[iarg]);
-    } else if (utils::strmatch(arg[iarg],"^gradv")) {
+    } else if (utils::strmatch(arg[iarg], "^grad/v")) {
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_gradv;
       col_index[i] = get_tensor_index(arg[iarg]);
     } else {
@@ -167,10 +169,10 @@ void ComputeRHEOPropertyAtom::compute_peratom()
     nmax = atom->nmax;
     if (nvalues == 1) {
       memory->destroy(vector_atom);
-      memory->create(vector_atom,nmax,"rheo/property/atom:vector");
+      memory->create(vector_atom, nmax, "rheo/property/atom:vector");
     } else {
       memory->destroy(array_atom);
-      memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array");
+      memory->create(array_atom, nmax, nvalues, "rheo/property/atom:array");
     }
   }
 
@@ -220,6 +222,21 @@ void ComputeRHEOPropertyAtom::pack_phase(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void ComputeRHEOPropertyAtom::pack_rho(int n)
+{
+  double *rho = atom->rho;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = rho[i];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputeRHEOPropertyAtom::pack_chi(int n)
 {
   double *chi = compute_interface->chi;
@@ -370,29 +387,29 @@ int ComputeRHEOPropertyAtom::get_tensor_index(char* option)
   int dim = domain->dimension;
   int dim_error = 0;
 
-  if (utils::strmatch(option,"xx$")) {
+  if (utils::strmatch(option, "xx$")) {
     index = 0;
-  } else if (utils::strmatch(option,"xy$")) {
+  } else if (utils::strmatch(option, "xy$")) {
     index = 1;
-  } else if (utils::strmatch(option,"xz$")) {
+  } else if (utils::strmatch(option, "xz$")) {
     index = 2;
     if (dim == 2) dim_error = 1;
-  } else if (utils::strmatch(option,"yx$")) {
+  } else if (utils::strmatch(option, "yx$")) {
     if (dim == 2) index = 2;
     else index = 3;
-  } else if (utils::strmatch(option,"yy$")) {
+  } else if (utils::strmatch(option, "yy$")) {
     if (dim == 2) index = 3;
     else index = 4;
-  } else if (utils::strmatch(option,"yz$")) {
+  } else if (utils::strmatch(option, "yz$")) {
     index = 5;
     if (dim == 2) dim_error = 1;
-  } else if (utils::strmatch(option,"zx$")) {
+  } else if (utils::strmatch(option, "zx$")) {
     index = 6;
     if (dim == 2) dim_error = 1;
-  } else if (utils::strmatch(option,"zy$")) {
+  } else if (utils::strmatch(option, "zy$")) {
     index = 7;
     if (dim == 2) dim_error = 1;
-  } else if (utils::strmatch(option,"zz$")) {
+  } else if (utils::strmatch(option, "zz$")) {
     index = 8;
     if (dim == 2) dim_error = 1;
   } else {

src/RHEO/compute_rheo_property_atom.h

@@ -43,6 +43,7 @@ class ComputeRHEOPropertyAtom : public Compute {
   FnPtrPack *pack_choice;    // ptrs to pack functions
 
   void pack_phase(int);
+  void pack_rho(int);
   void pack_chi(int);
   void pack_surface(int);
   void pack_surface_r(int);

src/RHEO/compute_rheo_vshift.cpp

@@ -232,7 +232,7 @@ void ComputeRHEOVShift::correct_surfaces()
   int nlocal = atom->nlocal;
   int dim = domain->dimension;
 
-  double nx, ny, nz, vx, vy, vz;
+  double nx, ny, nz, vx, vy, vz, dot;
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
@@ -240,11 +240,20 @@ void ComputeRHEOVShift::correct_surfaces()
         ny = nsurface[i][1];
         vx = vshift[i][0];
         vy = vshift[i][1];
-        vz = vshift[i][2];
+
+        dot = nx * vx + ny * vy;
+        if (dim == 3) {
+          nz = nsurface[i][2];
+          vz = vshift[i][2];
+          dot += nz * vz;
+        }
+
+        // Allowing shifting into the bulk
+        if (dot < 0.0) continue;
+
         vshift[i][0] = (1 - nx * nx) * vx - nx * ny * vy;
         vshift[i][1] = (1 - ny * ny) * vy - nx * ny * vx;
-        if (dim > 2) {
+        if (dim == 3) {
-          nz = nsurface[i][2];
           vshift[i][0] -= nx * nz * vz;
           vshift[i][1] -= ny * nz * vz;
           vshift[i][2] = (1 - nz * nz) * vz - nz * ny * vy - nx * nz * vx;

src/RHEO/fix_rheo.cpp

@@ -351,7 +351,6 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   // Shifting atoms
   if (shift_flag) {
-    compute_vshift->correct_surfaces(); // Could this be moved to preforce after the surface fix runs?
     for (i = 0; i < nlocal; i++) {
 
       if (status[i] & STATUS_NO_SHIFT) continue;
@@ -407,6 +406,9 @@ void FixRHEO::pre_force(int /*vflag*/)
 
   // Calculate surfaces, update status
   if (surface_flag) compute_surface->compute_peratom();
+
+  if (shift_flag)
+    compute_vshift->correct_surfaces();
 }
 
 /* ---------------------------------------------------------------------- */