Improvements & bug fixes to fix def/press

doc/src/fix_deform_pressure.rst

@@ -29,10 +29,12 @@ Syntax
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *xy*, *xz*, *yz* args = style value
-         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure*
+         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure* or *erate/rescale*
            *pressure* values = target gain
              target = target pressure (pressure units)
              gain = proportional gain constant (1/(time * pressure) or 1/time units)
+           *erate/rescale* value = R
+             R = engineering shear strain rate (1/time units)
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *box* = style value
@@ -159,6 +161,21 @@ details of a simulation and testing different values is
 recommended. One can also apply a maximum limit to the magnitude of
 the applied strain using the :ref:`max/rate <deform_max_rate>` option.
 
+The *erate/rescale* style operates similarly to the *erate* style with
+a specified strain rate in units of 1/time. The difference is that
+the change in the tilt factor will depend on the current length of
+the box perpendicular to the shear direction, L, instead of the
+original length, L0. The tilt factor T as a function of time will
+change as
+
+.. parsed-literal::
+
+   T(t) = T(t-1) + L\*erate\* \Delta t
+
+where T(t-1) is the tilt factor on the previous timestep and :math:`\Delta t`
+is the timestep size. This option may be useful in scenarios where
+L changes in time.
+
 ----------
 
 The *box* parameter provides an additional control over the *x*, *y*,

src/EXTRA-FIX/fix_deform_pressure.cpp

@@ -110,6 +110,12 @@ FixDeformPressure::FixDeformPressure(LAMMPS *lmp, int narg, char **arg) :
         }
         set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         i += 4;
+      } else if (strcmp(arg[iarg + 1], "erate/rescale") == 0) {
+        if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure erate/rescale", error);
+        set[index].style = ERATERS;
+        set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+        iarg += 3;
+        i += 3;
       } else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg + 1]);
 
     } else if (strcmp(arg[iarg], "box") == 0) {
@@ -424,16 +430,31 @@ void FixDeformPressure::init()
     if (!pressure)
       error->all(FLERR, "Pressure ID {} for fix deform/pressure does not exist", id_press);
   }
+
+  // if yz [3] changes and will cause box flip, then xy [5] cannot be changing
+  // this is b/c the flips would induce continuous changes in xz
+  //   in order to keep the edge vectors of the flipped shape matrix
+  //   an integer combination of the edge vectors of the unflipped shape matrix
+  // error if style PRESSURE/ERATEER for yz, can't calculate if box flip occurs
+
+  if (set[3].style && set[5].style) {
+    int flag = 0;
+    double lo,hi;
+    if (flipflag && set[3].style == PRESSURE)
+      error->all(FLERR, "Fix {} cannot use yz pressure with xy", style);
+    if (flipflag && set[3].style == ERATERS)
+      error->all(FLERR, "Fix {} cannot use yz erate/rescale with xy", style);
+  }
 }
 
 /* ----------------------------------------------------------------------
-   compute T,P if needed before integrator starts
+   compute T,P before integrator starts
 ------------------------------------------------------------------------- */
 
 void FixDeformPressure::setup(int /*vflag*/)
 {
-  // trigger virial computation on next timestep
+  // trigger virial computation, if needed, on next timestep
-  if (pressure_flag) pressure->addstep(update->ntimestep+1);
+  if (pressure_flag) pressure->addstep(update->ntimestep + 1);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -446,7 +467,20 @@ void FixDeformPressure::end_of_step()
 
   // set new box size for strain-based dims
 
-  if (strain_flag) FixDeform::apply_strain();
+  if (strain_flag) {
+    FixDeform::apply_strain();
+
+    for (int i = 3; i < 6; i++) {
+      if (set[i].style == ERATERS) {
+        double L = domain->zprd;
+        if (i == 5) L = domain->yprd;
+
+        h_rate[i] = set[i].rate * L;
+        set_extra[i].cumulative_shift += update->dt * h_rate[i];
+        set[i].tilt_target = set[i].tilt_start + set_extra[i].cumulative_shift;
+      }
+    }
+  }
 
   // set new box size for pressure-based dims
 
@@ -479,12 +513,33 @@ void FixDeformPressure::end_of_step()
     for (int i = 0; i < 3; i++) {
       set_extra[i].prior_pressure = pressure->vector[i];
       set_extra[i].prior_rate = ((set[i].hi_target - set[i].lo_target) /
-                           (domain->boxhi[i] - domain->boxlo[i]) - 1.0)  / update->dt;
+                           domain->prd[i] - 1.0)  / update->dt;
     }
   }
 
   if (varflag) modify->addstep_compute(update->ntimestep + nevery);
 
+  // If tilting while evolving linear dimension, sum remapping effects
+  // otherwise, update_domain() will inaccurately use the current
+  // linear dimension to apply prior remappings
+
+  for (int i = 3; i < 6; i++) {
+    int idenom = 0;
+    if (i == 3) idenom = 1;
+    if (set[i].style && (set_box.style || set[idenom].style)) {
+      // Add prior remappings. If the box remaps this timestep, don't
+      // add it yet so update_domain() will first detect the remapping
+      set[i].tilt_target += set_extra[i].cumulative_remap;
+
+      // Update remapping for next timestep
+      double prd = set[idenom].hi_target - set[idenom].lo_target;
+      double prdinv = 1.0 / prd;
+      if (set[i].tilt_target * prdinv < -0.5)
+        set_extra[i].cumulative_remap += prd;
+      if (set[i].tilt_target * prdinv > 0.5)
+        set_extra[i].cumulative_remap -= prd;
+    }
+  }
 
   FixDeform::update_domain();
 
@@ -556,26 +611,24 @@ void FixDeformPressure::apply_pressure()
 
     h_ratelo[i] = -0.5 * h_rate[i];
 
-    double offset = 0.5 * (domain->boxhi[i] - domain->boxlo[i]) * (1.0 + update->dt * h_rate[i]);
+    double shift = domain->prd[i] * update->dt * h_rate[i];
-    set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - offset;
+    set_extra[i].cumulative_shift += shift;
-    set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + offset;
+    set[i].lo_target = set[i].lo_start - 0.5 * set_extra[i].cumulative_shift;
+    set[i].hi_target = set[i].hi_start + 0.5 * set_extra[i].cumulative_shift;
   }
 
   for (int i = 3; i < 6; i++) {
     if (set[i].style != PRESSURE) continue;
 
-    double L, tilt, pcurrent;
+    double L, pcurrent;
     if (i == 3) {
       L = domain->zprd;
-      tilt = domain->yz;
       pcurrent = tensor[5];
     } else if (i == 4) {
       L = domain->zprd;
-      tilt = domain->xz + update->dt;
       pcurrent = tensor[4];
     } else {
       L = domain->yprd;
-      tilt = domain->xy;
       pcurrent = tensor[3];
     }
 
@@ -592,7 +645,8 @@ void FixDeformPressure::apply_pressure()
       if (fabs(h_rate[i]) > max_h_rate)
         h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
 
-    set[i].tilt_target = tilt + update->dt * h_rate[i];
+    set_extra[i].cumulative_shift += update->dt * h_rate[i];
+    set[i].tilt_target = set[i].tilt_start + set_extra[i].cumulative_shift;
   }
 }
 
@@ -629,9 +683,9 @@ void FixDeformPressure::apply_volume()
         double dt = update->dt;
         double e1i = set_extra[i].prior_rate;
         double e2i = set_extra[fixed].prior_rate;
-        double L1i = domain->boxhi[i] - domain->boxlo[i];
+        double L1i = domain->prd[i];
-        double L2i = domain->boxhi[fixed] - domain->boxlo[fixed];
+        double L2i = domain->prd[fixed];
-        double L3i = domain->boxhi[dynamic1] - domain->boxlo[dynamic1];
+        double L3i = domain->prd[dynamic1];
         double L3 = (set[dynamic1].hi_target - set[dynamic1].lo_target);
         double Vi = L1i * L2i * L3i;
         double V = L3 * L1i * L2i;
@@ -680,7 +734,7 @@ void FixDeformPressure::apply_volume()
       }
     }
 
-    h_rate[i] = (2.0 * shift / (domain->boxhi[i] - domain->boxlo[i]) - 1.0) / update->dt;
+    h_rate[i] = (2.0 * shift / domain->prd[i] - 1.0) / update->dt;
     h_ratelo[i] = -0.5 * h_rate[i];
 
     set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
@@ -742,7 +796,7 @@ void FixDeformPressure::apply_box()
       set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
 
       // Recalculate h_rate
-      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / domain->prd[i] - 1.0;
       h_rate[i] /= update->dt;
       h_ratelo[i] = -0.5 * h_rate[i];
     }
@@ -767,14 +821,14 @@ void FixDeformPressure::apply_box()
       if (fabs(v_rate) > max_h_rate)
         v_rate = max_h_rate * v_rate / fabs(v_rate);
 
-    scale = (1.0 + update->dt * v_rate);
     for (i = 0; i < 3; i++) {
-      shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
+      shift = (set[i].hi_target - set[i].lo_target) * update->dt * v_rate;
-      set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
+      set_extra[6].cumulative_vshift[i] += shift;
-      set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
+      set[i].lo_target -= 0.5 * set_extra[6].cumulative_vshift[i];
+      set[i].hi_target += 0.5 * set_extra[6].cumulative_vshift[i];
 
       // Recalculate h_rate
-      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / domain->prd[i] - 1.0;
       h_rate[i] /= update->dt;
       h_ratelo[i] = -0.5 * h_rate[i];
     }

src/EXTRA-FIX/fix_deform_pressure.h

@@ -51,6 +51,9 @@ class FixDeformPressure : public FixDeform {
   struct SetExtra {
     double ptarget, pgain;
     double prior_pressure, prior_rate;
+    double cumulative_shift;
+    double cumulative_vshift[3];
+    double cumulative_remap;
     int saved;
     char *pstr;
     int pvar, pvar_flag;

src/fix_deform.cpp

@@ -63,7 +63,7 @@ irregular(nullptr), set(nullptr)
   int nskip;
   if (utils::strmatch(style, "^deform/pressure")) {
     child_parameters.insert("box");
-    child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"volume", 2}});
+    child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"erate/rescale", 3}, {"volume", 2}});
   }
 
   // set defaults

src/fix_deform.h

@@ -29,7 +29,7 @@ class FixDeform : public Fix {
   int remapflag;     // whether x,v are remapped across PBC
   int dimflag[6];    // which dims are deformed
 
-  enum { NONE, FINAL, DELTA, SCALE, VEL, ERATE, TRATE, VOLUME, WIGGLE, VARIABLE, PRESSURE, PMEAN };
+  enum { NONE, FINAL, DELTA, SCALE, VEL, ERATE, TRATE, VOLUME, WIGGLE, VARIABLE, PRESSURE, PMEAN, ERATERS };
   enum { ONE_FROM_ONE, ONE_FROM_TWO, TWO_FROM_ONE };
 
   FixDeform(class LAMMPS *, int, char **);