Fixing errors in oxide model

src/RHEO/compute_rheo_property_atom.cpp

@@ -103,7 +103,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
       if (avec_index[i] < 0)
         error->all(FLERR,
                    "Invalid keyword {} for atom style {} in compute rheo/property/atom command ",
-                   atom->get_style(), arg[iarg]);
+                   arg[iarg], atom->get_style());
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
 
       if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1;

src/RHEO/compute_rheo_surface.cpp

@@ -152,16 +152,12 @@ void ComputeRHEOSurface::compute_peratom()
     grow_arrays(atom->nmax);
 
   size_t nbytes = nmax_store * sizeof(double);
-  memset(&divr, 0, nbytes);
+  memset(&divr[0], 0, nbytes);
-  memset(&rsurface, 0, nbytes);
+  memset(&rsurface[0], 0, nbytes);
-  memset(&nsurface, 0, 3 * nbytes);
+  memset(&nsurface[0][0], 0, dim * nbytes);
-  memset(&gradC, 0, 3 * 3 * nbytes);
+  memset(&gradC[0][0], 0, dim * dim * nbytes);
-  memset(&B, 0, 3 * 3 * nbytes);
+  memset(&B[0][0], 0, dim * dim * nbytes);
 
-  // Remove surface settings
-  int nall = nlocal + atom->nghost;
-  for (i = 0; i < nall; i++)
-    status[i] &= SURFACEMASK;
 
   // loop over neighbors to calculate the average orientation of neighbors
   for (ii = 0; ii < inum; ii++) {
@@ -208,7 +204,7 @@ void ComputeRHEOSurface::compute_peratom()
 
         wp = compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
 
-        for (a = 0; a < dim; a++){
+        for (a = 0; a < dim; a++) {
           divr[i] -= dWij[a] * dx[a] * Volj;
           gradC[i][a] += dWij[a] * Volj;
         }
@@ -245,35 +241,32 @@ void ComputeRHEOSurface::compute_peratom()
     }
   }
 
-  // Find the free-surface
+  // Remove surface settings and assign new values
-  if (threshold_style == DIVR) {
+  int nall = nlocal + atom->nghost;
-    for (i = 0; i < nall; i++) {
+  int test;
-      if (mask[i] & groupbit) {
+
+  for (i = 0; i < nall; i++) {
+    status[i] &= SURFACEMASK;
+    if (mask[i] & groupbit) {
+      if (threshold_style == DIVR)
+        test = divr[i] < threshold_divr;
+      else
+        test = coordination[i] < threshold_z;
+
+      if (test) {
+        if (coordination[i] < threshold_splash)
+          status[i] |= STATUS_SPLASH;
+        else
+          status[i] |= STATUS_SURFACE;
+        rsurface[i] = 0.0;
+      } else {
         status[i] |= STATUS_BULK;
         rsurface[i] = cut;
-        if (divr[i] < threshold_divr) {
-          status[i] |= STATUS_SURFACE;
-          rsurface[i] = 0.0;
-          if (coordination[i] < threshold_splash)
-            status[i] |= STATUS_SPLASH;
-        }
-      }
-    }
-  } else {
-    for (i = 0; i < nall; i++) {
-      if (mask[i] & groupbit) {
-        status[i] |= STATUS_BULK;
-        rsurface[i] = cut;
-        if (coordination[i] < threshold_z) {
-          status[i] |= STATUS_SURFACE;
-          rsurface[i] = 0.0;
-          if (coordination[i] < threshold_splash)
-            status[i] |= STATUS_SPLASH;
-        }
       }
     }
   }
 
+
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     xtmp = x[i][0];
@@ -297,7 +290,8 @@ void ComputeRHEOSurface::compute_peratom()
           status[i] |= STATUS_LAYER;
         }
 
-        if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
+        if (status[j] & STATUS_SURFACE)
+          rsurface[i] = MIN(rsurface[i], sqrt(rsq));
 
 
         if (j < nlocal || newton) {
@@ -306,7 +300,8 @@ void ComputeRHEOSurface::compute_peratom()
             status[j] |= STATUS_LAYER;
           }
 
-          if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
+          if (status[i] & STATUS_SURFACE)
+            rsurface[j] = MIN(rsurface[j], sqrt(rsq));
         }
       }
     }
@@ -351,7 +346,8 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
   int i,a,b,k,j,m;
   int dim = domain->dimension;
   int *status = atom->status;
-  int temp;
+  int tmp1;
+  double tmp2;
 
   m = 0;
   for (i = 0; i < n; i++) {
@@ -362,12 +358,13 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
         for (b = 0; b < dim; b ++)
           gradC[j][a * dim + b] += buf[m++];
     } else if (comm_stage == 1) {
-
+      tmp1 = (int) buf[m++];
-      temp = (int) buf[m++];
+      if ((status[j] & STATUS_BULK) && (tmp1 & STATUS_LAYER)) {
-      if ((status[j] & STATUS_BULK) && (temp & STATUS_LAYER))
+        status[j] &= SURFACEMASK;
-        status[j] = temp;
+        status[j] |= STATUS_LAYER;
-
+      }
-      rsurface[j] = MIN(rsurface[j], buf[m++]);
+      tmp2 = buf[m++];
+      rsurface[j] = MIN(rsurface[j], tmp2);
     }
   }
 }

src/RHEO/fix_rheo.cpp

@@ -133,7 +133,7 @@ FixRHEO::~FixRHEO()
   if (compute_kernel) modify->delete_compute("rheo_kernel");
   if (compute_grad) modify->delete_compute("rheo_grad");
   if (compute_interface) modify->delete_compute("rheo_interface");
-  if (compute_surface) modify->delete_compute("compute_surface");
+  if (compute_surface) modify->delete_compute("rheo_surface");
   if (compute_rhosum) modify->delete_compute("rheo_rhosum");
   if (compute_vshift) modify->delete_compute("rheo_vshift");
 }
@@ -211,7 +211,10 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
     error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
 
   // Calculate surfaces
-  if (surface_flag) compute_surface->compute_peratom();
+  if (surface_flag) {
+    compute_kernel->compute_coordination();
+    compute_surface->compute_peratom();
+  }
 
   pre_force(0);
 }

src/RHEO/pair_rheo_react.cpp

@@ -61,13 +61,13 @@ PairRHEOReact::PairRHEOReact(LAMMPS *lmp) : Pair(lmp),
   // between timesteps (fix property atom will handle callbacks)
 
   int tmp1, tmp2;
-  int index = atom->find_custom("react_nbond", tmp1, tmp2);
+  index_nb = atom->find_custom("react_nbond", tmp1, tmp2);
-  if (index == -1) {
+  if (index_nb == -1) {
     id_fix = utils::strdup("pair_rheo_react_fix_property_atom");
     modify->add_fix(fmt::format("{} all property/atom i_react_nbond", id_fix));
-    index = atom->find_custom("nbond", tmp1, tmp2);
+    index_nb = atom->find_custom("react_nbond", tmp1, tmp2);
   }
-  nbond = atom->ivector[index];
+  nbond = atom->ivector[index_nb];
 
   //Store non-persistent per atom quantities, intermediate
 
@@ -79,9 +79,10 @@ PairRHEOReact::PairRHEOReact(LAMMPS *lmp) : Pair(lmp),
 
 PairRHEOReact::~PairRHEOReact()
 {
-  if (modify->nfix && fix_history) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT");
+  if (modify->nfix && fix_history) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT" + std::to_string(instance_me));
-  if (modify->nfix && fix_dummy) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY");
+  if (modify->nfix && fix_dummy) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me));
-  if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_RHEO_REACT");
+  if (modify->nfix) modify->delete_fix(id_fix);
+  delete[] id_fix;
 
   if (allocated) {
     memory->destroy(setflag);
@@ -146,7 +147,7 @@ void PairRHEOReact::compute(int eflag, int vflag)
   }
 
   size_t nbytes = nmax_store * sizeof(int);
-  memset(&dbond, 0, nbytes);
+  memset(&dbond[0], 0, nbytes);
 
   // loop over neighbors of my atoms
   for (ii = 0; ii < inum; ii++) {
@@ -368,11 +369,11 @@ void PairRHEOReact::coeff(int narg, char **arg)
 
 void PairRHEOReact::init_style()
 {
-  int irequest = neighbor->request(this,instance_me);
+  int irequest = neighbor->request(this, instance_me);
   //neighbor->requests[irequest]->history = 1;
 
   if (fix_history == nullptr) {
-    auto cmd = fmt::format("NEIGH_HISTORY_RHEO_REACT {} all NEIGH_HISTORY {}", instance_me, size_history);
+    auto cmd = fmt::format("NEIGH_HISTORY_RHEO_REACT{} all NEIGH_HISTORY {}", instance_me, size_history);
     fix_history = dynamic_cast<FixNeighHistory *>(
         modify->replace_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me), cmd, 1));
     fix_history->pair = this;