Cleaning up math, fixing tension bug, patching bond creation

src/RHEO/compute_rheo_interface.cpp

@@ -26,6 +26,7 @@
 #include "force.h"
 #include "fix_rheo.h"
 #include "fix_rheo_pressure.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
 #include "neighbor.h"
@@ -36,6 +37,7 @@
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
+using namespace MathExtra;
 
 static constexpr double EPSILON = 1e-1;
 
@@ -107,8 +109,8 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr)
 
 void ComputeRHEOInterface::compute_peratom()
 {
-  int i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match;
+  int a, i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match;
-  double xtmp, ytmp, ztmp, delx, dely, delz, rsq, w, dot;
+  double xtmp, ytmp, ztmp, rsq, w, dot, dx[3];
 
   int inum, *ilist, *jlist, *numneigh, **firstneigh;
   int nlocal = atom->nlocal;
@@ -153,15 +155,15 @@ void ComputeRHEOInterface::compute_peratom()
       j = jlist[jj];
       j &= NEIGHMASK;
 
-      delx = xtmp - x[j][0];
+      dx[0] = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
+      dx[1] = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
+      dx[2] = ztmp - x[j][2];
-      rsq = delx * delx + dely * dely + delz * delz;
+      rsq = lensq3(dx);
 
       if (rsq < cutsq) {
         jtype = type[j];
         fluidj = !(status[j] & PHASECHECK);
-        w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq));
+        w = compute_kernel->calc_w_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq));
 
         status_match = 0;
         norm[i] += w;
@@ -172,9 +174,9 @@ void ComputeRHEOInterface::compute_peratom()
           chi[i] += w;
         } else {
           if (!fluidi) {
-            dot = (-fp_store[j][0] + fp_store[i][0]) * delx;
+            dot = 0;
-            dot += (-fp_store[j][1] + fp_store[i][1]) * dely;
+            for (a = 0; a < 3; a++)
-            dot += (-fp_store[j][2] + fp_store[i][2]) * delz;
+              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;
@@ -187,9 +189,9 @@ void ComputeRHEOInterface::compute_peratom()
             chi[j] += w;
           } else {
             if (!fluidj) {
-              dot = (-fp_store[i][0] + fp_store[j][0]) * delx;
+              dot = 0;
-              dot += (-fp_store[i][1] + fp_store[j][1]) * dely;
+              for (a = 0; a < 3; a++)
-              dot += (-fp_store[i][2] + fp_store[j][2]) * delz;
+                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;
@@ -225,7 +227,7 @@ void ComputeRHEOInterface::compute_peratom()
 
 int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i,j,k,m;
+  int i, j, k, m;
   m = 0;
   double *rho = atom->rho;
 
@@ -267,7 +269,7 @@ void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
 
 int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,k,m,last;
+  int i, k, m, last;
   double *rho = atom->rho;
 
   m = 0;
@@ -285,7 +287,7 @@ int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,k,j,m;
+  int i, k, j, m;
   double *rho = atom->rho;
   int *status = atom->status;
   m = 0;
@@ -350,23 +352,19 @@ void ComputeRHEOInterface::store_forces()
     for (const auto &fix : fixlist) {
       for (int i = 0; i < atom->nlocal; i++) {
         minv = 1.0 / mass[type[i]];
-        if (mask[i] & fix->groupbit) {
+        if (mask[i] & fix->groupbit)
-          fp_store[i][0] = f[i][0] * minv;
+          for (int a = 0; a < 3; a++)
-          fp_store[i][1] = f[i][1] * minv;
+            fp_store[i][a] = f[i][a] * minv;
-          fp_store[i][2] = f[i][2] * minv;
+        else
-        } else {
+          for (int a = 0; a < 3; a++)
-          fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv;
+            fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
-          fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv;
-          fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv;
-        }
       }
     }
   } else {
     for (int i = 0; i < atom->nlocal; i++) {
       minv = 1.0 / mass[type[i]];
-      fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv;
+      for (int a = 0; a < 3; a++)
-      fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv;
+        fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
-      fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv;
     }
   }
 

src/RHEO/fix_rheo.cpp

@@ -422,12 +422,6 @@ void FixRHEO::pre_force(int /*vflag*/)
     if (mask[i] & groupbit)
       status[i] &= OPTIONSMASK;
 
-  // Reinstate temporary options
-  for (int i = 0; i < nall; i++)
-    if (mask[i] & groupbit)
-      if (status[i] & STATUS_SOLID)
-        status[i] |= STATUS_NO_SHIFT;
-
   // Calculate surfaces, update status
   if (surface_flag) {
     compute_surface->compute_peratom();

src/RHEO/fix_rheo_tension.cpp

@@ -467,7 +467,7 @@ void FixRHEOTension::pre_force(int vflag)
 
     if (wsame[i] < wmin) continue;
 
-    weight = MAX(1.0, wsame[i] * wmin_inv);
+    weight = MIN(1.0, wsame[i] * wmin_inv); //MAX -> MIN 2/14/24
     itype = type[i];
 
     if (norm[i] != 0)

src/RHEO/fix_rheo_thermal.cpp

@@ -423,6 +423,7 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
   double *energy = atom->esph;
   double *temperature = atom->temperature;
   int *type = atom->type;
+  int *status = atom->status;
 
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
@@ -444,6 +445,11 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
       }
     }
   }
+
+  // Add temporary options, wiped by preceding fix rheo preforce
+  for (int i = 0; i < nall; i++)
+    if (status[i] & STATUS_SOLID)
+      status[i] |= STATUS_NO_SHIFT;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -481,39 +487,79 @@ void FixRHEOThermal::break_bonds()
   tagint **bond_atom = atom->bond_atom;
   int *num_bond = atom->num_bond;
 
-  int nlocal = atom->nlocal;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
 
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+
+  // Rapidly delete all bonds for local atoms that melt (no shifting)
   for (int i = 0; i < nlocal; i++) {
+    if (!(status[i] & STATUS_MELTING)) continue;
     for (m = 0; m < num_bond[i]; m++) {
       j = atom->map(bond_atom[i][m]);
-      if (!(status[i] & STATUS_MELTING) && !(status[j] & STATUS_MELTING)) continue;
+      bond_type[i][m] = 0;
 
       if (n_histories > 0)
         for (auto &ihistory: histories)
-          dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i, num_bond[i] - 1);
+          dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i, m);
 
-      if (fix_update_special_bonds) fix_update_special_bonds->add_broken_bond(i, j);
+      if (fix_update_special_bonds)
+        fix_update_special_bonds->add_broken_bond(i, j);
+    }
+    num_bond[i] = 0;
+  }
 
-      // For non-melting neighbors, selectively delete bond if necessary
+  // Update bond list and break solid-melted bonds
-      if (j >= nlocal || (status[j] & STATUS_MELTING)) continue;
+  for (n = 0; n < nbondlist; n++) {
-      for (n = 0; n < num_bond[j]; n++) {
+
-        if (bond_atom[j][n] == tag[i]) {
+    // skip bond if already broken
-          bond_type[j][n] = 0;
+    if (bondlist[n][2] <= 0) continue;
-          nmax = num_bond[j] - 1;
+    i = bondlist[n][0];
-          bond_type[j][n] = bond_type[j][nmax];
+    j = bondlist[n][1];
-          bond_atom[j][n] = bond_atom[j][nmax];
+
-          if (n_histories > 0) {
+    if (!(status[i] & STATUS_MELTING) && !(status[j] & STATUS_MELTING)) continue;
+
+    bondlist[n][2] = 0;
+
+    // Delete bonds for non-melted local atoms (shifting)
+    if (i < nlocal) {
+      for (m = 0; m < num_bond[i]; m++) {
+        if (bond_atom[i][m] == tag[j]) {
+          nmax = num_bond[i] - 1;
+          bond_type[i][m] = bond_type[i][nmax];
+          bond_atom[i][m] = bond_atom[i][nmax];
+          if (n_histories > 0)
             for (auto &ihistory: histories) {
-              dynamic_cast<FixBondHistory *>(ihistory)->shift_history(j, n, nmax);
+              auto fix_bond_history = dynamic_cast<FixBondHistory *>  (ihistory);
-              dynamic_cast<FixBondHistory *>(ihistory)->delete_history(j, nmax);
+              fix_bond_history->shift_history(i, m, nmax);
+              fix_bond_history->delete_history(i, nmax);
             }
-          }
+          bond_type[i][nmax] = 0;
+          num_bond[i]--;
+          break;
+        }
+      }
+    }
+
+    if (j < nlocal) {
+      for (m = 0; m < num_bond[j]; m++) {
+        if (bond_atom[j][m] == tag[i]) {
+          nmax = num_bond[j] - 1;
+          bond_type[j][m] = bond_type[j][nmax];
+          bond_atom[j][m] = bond_atom[j][nmax];
+          if (n_histories > 0)
+            for (auto &ihistory: histories) {
+              auto fix_bond_history = dynamic_cast<FixBondHistory *>  (ihistory);
+              fix_bond_history->shift_history(j, m, nmax);
+              fix_bond_history->delete_history(j, nmax);
+            }
+          bond_type[j][nmax] = 0;
           num_bond[j]--;
           break;
         }
       }
     }
-    num_bond[i] = 0;
   }
 }