diff --git a/src/VTK/dump_vtk.cpp b/src/VTK/dump_vtk.cpp
index a7aabe3f51..54dc0d4b0e 100644
--- a/src/VTK/dump_vtk.cpp
+++ b/src/VTK/dump_vtk.cpp
@@ -31,6 +31,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
+#include "fix_store.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -87,9 +88,9 @@ enum{X,Y,Z, // required for vtk, must come first
      Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
      OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
      TQX,TQY,TQZ,
-     VARIABLE,COMPUTE,FIX,INAME,DNAME,
+     COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY,
      ATTRIBUTES}; // must come last
-enum{LT,LE,GT,GE,EQ,NEQ};
+enum{LT,LE,GT,GE,EQ,NEQ,XOR};
 enum{VTK,VTP,VTU,PVTP,PVTU}; // file formats
 
 #define ONEFIELD 32
@@ -119,11 +120,10 @@ DumpVTK::DumpVTK(LAMMPS *lmp, int narg, char **arg) :
   // ioptional = start of additional optional args
   // only dump image and dump movie styles process optional args
 
-  ioptional = parse_fields(narg,arg);
+  ioptional = parse_fields(nargnew,earg);
 
-  if (ioptional < narg &&
-      strcmp(style,"image") != 0 && strcmp(style,"movie") != 0)
-    error->all(FLERR,"Invalid attribute in dump vtk command");
+  if (ioptional < nargnew)
+    error->all(FLERR,"Invalid attribute {} in dump vtk command", earg[ioptional]);
   size_one = pack_choice.size();
   current_pack_choice_key = -1;
 
@@ -210,38 +210,40 @@ void DumpVTK::init_style()
   else
     write_choice = &DumpVTK::write_vtk;
 
-  // find current ptr for each compute,fix,variable
+  // find current ptr for each compute,fix,variable and custom atom property
   // check that fix frequency is acceptable
 
-  int icompute;
   for (int i = 0; i < ncompute; i++) {
-    icompute = modify->find_compute(id_compute[i]);
+    int icompute = modify->find_compute(id_compute[i]);
     if (icompute < 0) error->all(FLERR,"Could not find dump vtk compute ID");
     compute[i] = modify->compute[icompute];
   }
 
-  int ifix;
   for (int i = 0; i < nfix; i++) {
-    ifix = modify->find_fix(id_fix[i]);
+    int ifix = modify->find_fix(id_fix[i]);
     if (ifix < 0) error->all(FLERR,"Could not find dump vtk fix ID");
     fix[i] = modify->fix[ifix];
     if (nevery % modify->fix[ifix]->peratom_freq)
       error->all(FLERR,"Dump vtk and fix not computed at compatible times");
   }
 
-  int ivariable;
   for (int i = 0; i < nvariable; i++) {
-    ivariable = input->variable->find(id_variable[i]);
+    int ivariable = input->variable->find(id_variable[i]);
     if (ivariable < 0)
       error->all(FLERR,"Could not find dump vtk variable name");
     variable[i] = ivariable;
   }
 
-  int icustom;
+  int icustom,flag,cols;
   for (int i = 0; i < ncustom; i++) {
-    icustom = atom->find_custom(id_custom[i],flag_custom[i]);
+    icustom = atom->find_custom(id_custom[i],flag,cols);
     if (icustom < 0)
       error->all(FLERR,"Could not find custom per-atom property ID");
+    custom[i] = icustom;
+    if (!flag && !cols) custom_flag[i] = IVEC;
+    else if (flag && !cols) custom_flag[i] = DVEC;
+    else if (!flag && cols) custom_flag[i] = IARRAY;
+    else if (flag && cols) custom_flag[i] = DARRAY;
   }
 
   // set index and check validity of region
@@ -275,7 +277,7 @@ int DumpVTK::count()
 
   // grow choose and variable vbuf arrays if needed
 
-  int nlocal = atom->nlocal;
+  const int nlocal = atom->nlocal;
   if (atom->nmax > maxlocal) {
     maxlocal = atom->nmax;
 
@@ -345,10 +347,10 @@ int DumpVTK::count()
   // un-choose if any threshold criterion isn't met
 
   if (nthresh) {
-    double *ptr;
+    double *ptr,*ptrhold;
+    double *values;
     double value;
-    int nstride;
-    int nlocal = atom->nlocal;
+    int nstride,lastflag;
 
     for (int ithresh = 0; ithresh < nthresh; ithresh++) {
 
@@ -635,26 +637,22 @@ int DumpVTK::count()
         nstride = 1;
       } else if (thresh_array[ithresh] == MUX) {
         if (!atom->mu_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
+          error->all(FLERR,"Threshold for an atom property that isn't allocated");
         ptr = &atom->mu[0][0];
         nstride = 4;
       } else if (thresh_array[ithresh] == MUY) {
         if (!atom->mu_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
+          error->all(FLERR,"Threshold for an atom property that isn't allocated");
         ptr = &atom->mu[0][1];
         nstride = 4;
       } else if (thresh_array[ithresh] == MUZ) {
         if (!atom->mu_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
+          error->all(FLERR,"Threshold for an atom property that isn't allocated");
         ptr = &atom->mu[0][2];
         nstride = 4;
       } else if (thresh_array[ithresh] == MU) {
         if (!atom->mu_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
+          error->all(FLERR,"Threshold for an atom property that isn't allocated");
         ptr = &atom->mu[0][3];
         nstride = 4;
 
@@ -753,9 +751,8 @@ int DumpVTK::count()
         nstride = 1;
 
       } else if (thresh_array[ithresh] == IVEC) {
-        int iwhich,flag,cols
         i = ATTRIBUTES + nfield + ithresh;
-        iwhich = atom->find_custom(id_custom[field2index[i]],flag,cols);
+        int iwhich = custom[field2index[i]];
         int *ivector = atom->ivector[iwhich];
         for (i = 0; i < nlocal; i++)
           dchoose[i] = ivector[i];
@@ -763,16 +760,14 @@ int DumpVTK::count()
         nstride = 1;
 
       } else if (thresh_array[ithresh] == DVEC) {
-        int iwhich,flag,cols;
         i = ATTRIBUTES + nfield + ithresh;
-        iwhich = atom->find_custom(id_custom[field2index[i]],flag,cols);
+        int iwhich = custom[field2index[i]];
         ptr = atom->dvector[iwhich];
         nstride = 1;
 
       } else if (thresh_array[ithresh] == IARRAY) {
-        int iwhich,flag,cols;
         i = ATTRIBUTES + nfield + ithresh;
-        iwhich = atom->find_custom(id_custom[field2index[i]],flag,cols);
+        int iwhich = custom[field2index[i]];
         int **iarray = atom->iarray[iwhich];
         int icol = argindex[i] - 1;
         for (i = 0; i < nlocal; i++)
@@ -781,43 +776,99 @@ int DumpVTK::count()
         nstride = 1;
 
       } else if (thresh_array[ithresh] == DARRAY) {
-        int iwhich,flag,cols;
         i = ATTRIBUTES + nfield + ithresh;
-        iwhich = atom->find_custom(id_custom[field2index[i]],flag,cols)
+        int iwhich = custom[field2index[i]];
         double **darray = atom->darray[iwhich];
         ptr = &darray[0][argindex[i]-1];
         nstride = atom->dcols[iwhich];
       }
 
       // unselect atoms that don't meet threshold criterion
+      // compare to single value or values stored in threshfix
+      // copy ptr attribute into thresh_fix if this is first comparison
 
-      value = thresh_value[ithresh];
+      if (thresh_last[ithresh] < 0) {
+        lastflag = 0;
+        value = thresh_value[ithresh];
+      } else {
+        lastflag = 1;
+        int ilast = thresh_last[ithresh];
+        values = thresh_fix[ilast]->vstore;
+        ptrhold = ptr;
+        if (thresh_first[ilast]) {
+          thresh_first[ilast] = 0;
+          for (i = 0; i < nlocal; i++, ptr += nstride) values[i] = *ptr;
+          ptr = ptrhold;
+        }
+      }
 
-      switch (thresh_op[ithresh]) {
-      case LT:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr >= value) choose[i] = 0;
-        break;
-      case LE:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr > value) choose[i] = 0;
-        break;
-      case GT:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr <= value) choose[i] = 0;
-        break;
-      case GE:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr < value) choose[i] = 0;
-        break;
-      case EQ:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr != value) choose[i] = 0;
-        break;
-      case NEQ:
-        for (i = 0; i < nlocal; i++, ptr += nstride)
-          if (choose[i] && *ptr == value) choose[i] = 0;
-        break;
+      if (thresh_op[ithresh] == LT) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr >= values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr >= value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == LE) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr > values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr > value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == GT) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr <= values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr <= value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == GE) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr < values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr < value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == EQ) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr != values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr != value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == NEQ) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr == values[i]) choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if (choose[i] && *ptr == value) choose[i] = 0;
+        }
+      } else if (thresh_op[ithresh] == XOR) {
+        if (lastflag) {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if ((choose[i] && *ptr == 0.0 && values[i] == 0.0) ||
+                (*ptr != 0.0 && values[i] != 0.0))
+              choose[i] = 0;
+        } else {
+          for (i = 0; i < nlocal; i++, ptr += nstride)
+            if ((choose[i] && *ptr == 0.0 && value == 0.0) ||
+                (*ptr != 0.0 && value != 0.0))
+              choose[i] = 0;
+        }
+      }
+
+      // update values stored in threshfix
+
+      if (lastflag) {
+        ptr = ptrhold;
+        for (i = 0; i < nlocal; i++, ptr += nstride) values[i] = *ptr;
       }
     }
   }
@@ -1754,15 +1805,16 @@ int DumpVTK::parse_fields(int narg, char **arg)
 
     } else {
 
-      int n,tmp;
+      int n,flag,cols;
       ArgInfo argi(arg[iarg],ArgInfo::COMPUTE|ArgInfo::FIX|ArgInfo::VARIABLE
-                   |ArgInfo::DVEC|ArgInfo::IVEC);
+                   |ArgInfo::DNAME|ArgInfo::INAME);
       argindex[ATTRIBUTES+i] = argi.get_index1();
+      auto aname = argi.get_name();
 
       switch (argi.get_type()) {
 
       case ArgInfo::UNKNOWN:
-        error->all(FLERR,"Invalid attribute in dump vtk command");
+        error->all(FLERR,"Invalid attribute in dump vtk command: {}",arg[iarg]);
         break;
 
       // compute value = c_ID
@@ -1772,21 +1824,19 @@ int DumpVTK::parse_fields(int narg, char **arg)
         pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_compute;
         vtype[ATTRIBUTES+i] = Dump::DOUBLE;
 
-        n = modify->find_compute(argi.get_name());
-        if (n < 0) error->all(FLERR,"Could not find dump vtk compute ID");
+        n = modify->find_compute(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump vtk compute ID: {}",aname);
         if (modify->compute[n]->peratom_flag == 0)
-          error->all(FLERR,"Dump vtk compute does not compute per-atom info");
+          error->all(FLERR,"Dump vtk compute {} does not compute per-atom info",aname);
         if (argi.get_dim() == 0 && modify->compute[n]->size_peratom_cols > 0)
-          error->all(FLERR,
-                     "Dump vtk compute does not calculate per-atom vector");
+          error->all(FLERR,"Dump vtk compute {} does not calculate per-atom vector",aname);
         if (argi.get_dim() > 0 && modify->compute[n]->size_peratom_cols == 0)
-          error->all(FLERR,
-                     "Dump vtk compute does not calculate per-atom array");
+          error->all(FLERR,"Dump vtk compute {} does not calculate per-atom array",aname);
         if (argi.get_dim() > 0 &&
             argi.get_index1() > modify->compute[n]->size_peratom_cols)
-          error->all(FLERR,"Dump vtk compute vector is accessed out-of-range");
+          error->all(FLERR,"Dump vtk compute {} vector is accessed out-of-range",aname);
 
-        field2index[ATTRIBUTES+i] = add_compute(argi.get_name());
+        field2index[ATTRIBUTES+i] = add_compute(aname);
         name[ATTRIBUTES+i] = arg[iarg];
         break;
 
@@ -1797,19 +1847,19 @@ int DumpVTK::parse_fields(int narg, char **arg)
         pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_fix;
         vtype[ATTRIBUTES+i] = Dump::DOUBLE;
 
-        n = modify->find_fix(argi.get_name());
-        if (n < 0) error->all(FLERR,"Could not find dump vtk fix ID");
+        n = modify->find_fix(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump vtk fix ID: {}",aname);
         if (modify->fix[n]->peratom_flag == 0)
-          error->all(FLERR,"Dump vtk fix does not compute per-atom info");
+          error->all(FLERR,"Dump vtk fix {} does not compute per-atom info",aname);
         if (argi.get_dim() == 0 && modify->fix[n]->size_peratom_cols > 0)
-          error->all(FLERR,"Dump vtk fix does not compute per-atom vector");
+          error->all(FLERR,"Dump vtk fix {} does not compute per-atom vector",aname);
         if (argi.get_dim() > 0 && modify->fix[n]->size_peratom_cols == 0)
-          error->all(FLERR,"Dump vtk fix does not compute per-atom array");
+          error->all(FLERR,"Dump vtk fix {} does not compute per-atom array",aname);
         if (argi.get_dim() > 0 &&
             argi.get_index1() > modify->fix[n]->size_peratom_cols)
-          error->all(FLERR,"Dump vtk fix vector is accessed out-of-range");
+          error->all(FLERR,"Dump vtk fix {} vector is accessed out-of-range",aname);
 
-        field2index[ATTRIBUTES+i] = add_fix(argi.get_name());
+        field2index[ATTRIBUTES+i] = add_fix(aname);
         name[ATTRIBUTES+i] = arg[iarg];
         break;
 
@@ -1819,61 +1869,62 @@ int DumpVTK::parse_fields(int narg, char **arg)
         pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_variable;
         vtype[ATTRIBUTES+i] = Dump::DOUBLE;
 
-        n = input->variable->find(argi.get_name());
-        if (n < 0) error->all(FLERR,"Could not find dump vtk variable name");
+        n = input->variable->find(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump vtk variable name {}",aname);
         if (input->variable->atomstyle(n) == 0)
-          error->all(FLERR,"Dump vtk variable is not atom-style variable");
+          error->all(FLERR,"Dump vtk variable {} is not atom-style variable",aname);
 
-        field2index[ATTRIBUTES+i] = add_variable(argi.get_name());
+        field2index[ATTRIBUTES+i] = add_variable(aname);
         name[ATTRIBUTES+i] = arg[iarg];
         break;
 
-      // custom per-atom integer vector = i_ID
-
-      case ArgInfo::INAME:
-        pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_custom;
-        vtype[ATTRIBUTES+i] = Dump::INT;
-
-        tmp = -1;
-        n = atom->find_custom(argi.get_name(),tmp);
-        if (n < 0)
-          error->all(FLERR,"Could not find custom per-atom property ID");
-
-        if (tmp != 0)
-          error->all(FLERR,"Custom per-atom property ID is not integer");
-
-        field2index[ATTRIBUTES+i] = add_custom(argi.get_name(),0);
-        name[ATTRIBUTES+i] = arg[iarg];
-        break;
-
-      // custom per-atom floating point vector = d_ID
+      // custom per-atom floating point vector or array = d_ID d2_ID
 
       case ArgInfo::DNAME:
         pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_custom;
         vtype[ATTRIBUTES+i] = Dump::DOUBLE;
 
-        tmp = -1;
-        n = atom->find_custom(argi.get_name(),tmp);
+        n = atom->find_custom(aname,flag,cols);
+
         if (n < 0)
-          error->all(FLERR,"Could not find custom per-atom property ID");
-
-        if (tmp != 1)
-          error->all(FLERR,"Custom per-atom property ID is not floating point");
-
-        field2index[ATTRIBUTES+i] = add_custom(argi.get_name(),1);
+          error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
+        if (argindex[ATTRIBUTES+i] == 0) {
+          if (!flag || cols)
+            error->all(FLERR,"Property double vector {} for dump vtk does not exist",aname);
+        } else {
+          if (!flag || !cols)
+            error->all(FLERR,"Property double array {} for dump vtk does not exist",aname);
+          if (argindex[ATTRIBUTES+i] > atom->dcols[n])
+            error->all(FLERR,"Dump vtk property array {} is accessed out-of-range",aname);
+        }
+        field2index[ATTRIBUTES+i] = add_custom(aname,1);
         name[ATTRIBUTES+i] = arg[iarg];
         break;
 
-      // NEWSTYLE
-      // custom per-atom integer array = i2_ID
+      // custom per-atom integer vector or array = i_ID or i2_ID
 
-      case ArgInfo::IARRAY:
-        return iarg;
+      case ArgInfo::INAME:
+        pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_custom;
+        vtype[ATTRIBUTES+i] = Dump::INT;
 
-      // custom per-atom floating point array = d2_ID
+        n = atom->find_custom(aname,flag,cols);
 
-      case ArgInfo::DARRAY:
-        return iarg;
+        if (n < 0)
+          error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
+        if (argindex[ATTRIBUTES+i] == 0) {
+          if (flag || cols)
+            error->all(FLERR,"Property integer vector {} for dump vtk does not exist",aname);
+        } else {
+          if (flag || !cols)
+            error->all(FLERR,"Property integer array {} for dump vtk does not exist",aname);
+          if (argindex[ATTRIBUTES+i] > atom->icols[n])
+            error->all(FLERR,"Dump vtk property array {} is accessed out-of-range",aname);
+        }
+        field2index[ATTRIBUTES+i] = add_custom(aname,0);
+        name[ATTRIBUTES+i] = arg[iarg];
+        break;
+
+      // no match
 
       default:
         return iarg;
@@ -1948,12 +1999,10 @@ int DumpVTK::add_compute(const char *id)
 
   id_compute = (char **)
     memory->srealloc(id_compute,(ncompute+1)*sizeof(char *),"dump:id_compute");
-  delete [] compute;
+  delete[] compute;
   compute = new Compute*[ncompute+1];
 
-  int n = strlen(id) + 1;
-  id_compute[ncompute] = new char[n];
-  strcpy(id_compute[ncompute],id);
+  id_compute[ncompute] = utils::strdup(id);
   ncompute++;
   return ncompute-1;
 }
@@ -1973,12 +2022,10 @@ int DumpVTK::add_fix(const char *id)
 
   id_fix = (char **)
     memory->srealloc(id_fix,(nfix+1)*sizeof(char *),"dump:id_fix");
-  delete [] fix;
+  delete[] fix;
   fix = new Fix*[nfix+1];
 
-  int n = strlen(id) + 1;
-  id_fix[nfix] = new char[n];
-  strcpy(id_fix[nfix],id);
+  id_fix[nfix] = utils::strdup(id);
   nfix++;
   return nfix-1;
 }
@@ -1999,22 +2046,20 @@ int DumpVTK::add_variable(const char *id)
   id_variable = (char **)
     memory->srealloc(id_variable,(nvariable+1)*sizeof(char *),
                      "dump:id_variable");
-  delete [] variable;
+  delete[] variable;
   variable = new int[nvariable+1];
-  delete [] vbuf;
+  delete[] vbuf;
   vbuf = new double*[nvariable+1];
   for (int i = 0; i <= nvariable; i++) vbuf[i] = nullptr;
 
-  int n = strlen(id) + 1;
-  id_variable[nvariable] = new char[n];
-  strcpy(id_variable[nvariable],id);
+  id_variable[nvariable] = utils::strdup(id);
   nvariable++;
   return nvariable-1;
 }
 
 /* ----------------------------------------------------------------------
    add custom atom property to list used by dump
-   return index of where this property is in list
+   return index of where this property is in Atom class custom lists
    if already in list, do not add, just return index, else add to list
 ------------------------------------------------------------------------- */
 
@@ -2022,21 +2067,17 @@ int DumpVTK::add_custom(const char *id, int flag)
 {
   int icustom;
   for (icustom = 0; icustom < ncustom; icustom++)
-    if ((strcmp(id,id_custom[icustom]) == 0)
-        && (flag == flag_custom[icustom])) break;
+    if (strcmp(id,id_custom[icustom]) == 0) break;
   if (icustom < ncustom) return icustom;
 
-  id_custom = (char **)
-    memory->srealloc(id_custom,(ncustom+1)*sizeof(char *),"dump:id_custom");
-  flag_custom = (int *)
-    memory->srealloc(flag_custom,(ncustom+1)*sizeof(int),"dump:flag_custom");
-
-  int n = strlen(id) + 1;
-  id_custom[ncustom] = new char[n];
-  strcpy(id_custom[ncustom],id);
-  flag_custom[ncustom] = flag;
+  id_custom = (char **) memory->srealloc(id_custom,(ncustom+1)*sizeof(char *),"dump:id_custom");
+  custom = (int *) memory->srealloc(custom,(ncustom+1)*sizeof(int),"dump:custom");
+  custom_flag = (int *) memory->srealloc(custom_flag,(ncustom+1)*sizeof(int),"dump:custom_flag");
 
+  id_custom[ncustom] = utils::strdup(id);
+  custom_flag[ncustom] = flag;
   ncustom++;
+
   return ncustom-1;
 }
 
@@ -2050,21 +2091,17 @@ int DumpVTK::modify_param(int narg, char **arg)
     else {
       iregion = domain->find_region(arg[1]);
       if (iregion == -1)
-        error->all(FLERR,"Dump_modify region ID does not exist");
-      delete [] idregion;
-      int n = strlen(arg[1]) + 1;
-      idregion = new char[n];
-      strcpy(idregion,arg[1]);
+        error->all(FLERR,"Dump_modify region ID {} does not exist",arg[1]);
+      delete[] idregion;
+      idregion = utils::strdup(arg[1]);
     }
     return 2;
   }
 
   if (strcmp(arg[0],"label") == 0) {
      if (narg < 2) error->all(FLERR,"Illegal dump_modify command [label]");
-     delete [] label;
-     int n = strlen(arg[1]) + 1;
-     label = new char[n];
-     strcpy(label,arg[1]);
+     delete[] label;
+     label = utils::strdup(arg[1]);
      return 2;
    }
 
@@ -2076,23 +2113,29 @@ int DumpVTK::modify_param(int narg, char **arg)
 
   if (strcmp(arg[0],"element") == 0) {
     if (narg < ntypes+1)
-      error->all(FLERR,"Dump modify: number of element names do not match atom types");
-
-    if (typenames) {
-      for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
-      delete [] typenames;
-      typenames = nullptr;
-    }
+      error->all(FLERR,"Number of dump_modify element names does not match number of atom types");
 
+    for (int i = 1; i <= ntypes; i++) delete[] typenames[i];
+    delete[] typenames;
     typenames = new char*[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) {
-      int n = strlen(arg[itype]) + 1;
-      typenames[itype] = new char[n];
-      strcpy(typenames[itype],arg[itype]);
+      typenames[itype] = utils::strdup(arg[itype]);
     }
     return ntypes+1;
   }
 
+  if (strcmp(arg[0],"refresh") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
+    ArgInfo argi(arg[1],ArgInfo::COMPUTE);
+    if ((argi.get_type() != ArgInfo::COMPUTE) || (argi.get_dim() != 0))
+      error->all(FLERR,"Illegal dump_modify command");
+    if (refreshflag) error->all(FLERR,"Dump_modify can only have one refresh");
+
+    refreshflag = 1;
+    refresh = argi.copy_name();
+    return 2;
+  }
+
   if (strcmp(arg[0],"thresh") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     if (strcmp(arg[1],"none") == 0) {
@@ -2103,8 +2146,16 @@ int DumpVTK::modify_param(int narg, char **arg)
         thresh_array = nullptr;
         thresh_op = nullptr;
         thresh_value = nullptr;
+        thresh_last = nullptr;
+        for (int i = 0; i < nthreshlast; i++) {
+          modify->delete_fix(thresh_fixID[i]);
+          delete[] thresh_fixID[i];
+        }
+        thresh_fix = nullptr;
+        thresh_fixID = nullptr;
+        thresh_first = nullptr;
       }
-      nthresh = 0;
+      nthresh = nthreshlast = 0;
       return 2;
     }
 
@@ -2115,6 +2166,7 @@ int DumpVTK::modify_param(int narg, char **arg)
     memory->grow(thresh_array,nthresh+1,"dump:thresh_array");
     memory->grow(thresh_op,(nthresh+1),"dump:thresh_op");
     memory->grow(thresh_value,(nthresh+1),"dump:thresh_value");
+    memory->grow(thresh_last,(nthresh+1),"dump:thresh_last");
 
     // set attribute type of threshold
     // customize by adding to if statement
@@ -2197,98 +2249,125 @@ int DumpVTK::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY;
     else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ;
 
-    // compute value = c_ID
-    // if no trailing [], then arg is set to 0, else arg is between []
+    // compute or fix or variable or custom vector/array
 
-    else if (strncmp(arg[1],"c_",2) == 0) {
-      thresh_array[nthresh] = COMPUTE;
-      int n = strlen(arg[1]);
-      char *suffix = new char[n];
-      strcpy(suffix,&arg[1][2]);
+    else {
+      int n,flag,cols;
+      ArgInfo argi(arg[1],ArgInfo::COMPUTE|ArgInfo::FIX|ArgInfo::VARIABLE
+                   |ArgInfo::DNAME|ArgInfo::INAME);
+      argindex[ATTRIBUTES+nfield+nthresh] = argi.get_index1();
+      auto aname = argi.get_name();
 
-      char *ptr = strchr(suffix,'[');
-      if (ptr) {
-        if (suffix[strlen(suffix)-1] != ']')
-          error->all(FLERR,"Invalid attribute in dump modify command");
-        argindex[ATTRIBUTES+nfield+nthresh] = atoi(ptr+1);
-        *ptr = '\0';
-      } else argindex[ATTRIBUTES+nfield+nthresh] = 0;
+      switch (argi.get_type()) {
 
-      n = modify->find_compute(suffix);
-      if (n < 0) error->all(FLERR,"Could not find dump modify compute ID");
+      case ArgInfo::UNKNOWN:
+        error->all(FLERR,"Invalid attribute in dump modify command");
+        break;
 
-      if (modify->compute[n]->peratom_flag == 0)
-        error->all(FLERR,
-                   "Dump modify compute ID does not compute per-atom info");
-      if (argindex[ATTRIBUTES+nfield+nthresh] == 0 &&
-          modify->compute[n]->size_peratom_cols > 0)
-        error->all(FLERR,
-                   "Dump modify compute ID does not compute per-atom vector");
-      if (argindex[ATTRIBUTES+nfield+nthresh] > 0 &&
-          modify->compute[n]->size_peratom_cols == 0)
-        error->all(FLERR,
-                   "Dump modify compute ID does not compute per-atom array");
-      if (argindex[ATTRIBUTES+nfield+nthresh] > 0 &&
-          argindex[ATTRIBUTES+nfield+nthresh] > modify->compute[n]->size_peratom_cols)
-        error->all(FLERR,"Dump modify compute ID vector is not large enough");
+      // compute value = c_ID
+      // if no trailing [], then arg is set to 0, else arg is between []
 
-      field2index[ATTRIBUTES+nfield+nthresh] = add_compute(suffix);
-      delete [] suffix;
+      case ArgInfo::COMPUTE:
+        thresh_array[nthresh] = COMPUTE;
+        n = modify->find_compute(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump modify compute ID: {}",aname);
 
-    // fix value = f_ID
-    // if no trailing [], then arg is set to 0, else arg is between []
+        if (modify->compute[n]->peratom_flag == 0)
+          error->all(FLERR,"Dump modify compute ID {} does not compute per-atom info",aname);
+        if (argi.get_dim() == 0 && modify->compute[n]->size_peratom_cols > 0)
+          error->all(FLERR,"Dump modify compute ID {} does not compute per-atom vector",aname);
+        if (argi.get_index1() > 0 && modify->compute[n]->size_peratom_cols == 0)
+          error->all(FLERR,"Dump modify compute ID {} does not compute per-atom array",aname);
+        if (argi.get_index1() > 0 &&
+            argi.get_index1() > modify->compute[n]->size_peratom_cols)
+          error->all(FLERR,"Dump modify compute ID {} vector is not large enough",aname);
 
-    } else if (strncmp(arg[1],"f_",2) == 0) {
-      thresh_array[nthresh] = FIX;
-      int n = strlen(arg[1]);
-      char *suffix = new char[n];
-      strcpy(suffix,&arg[1][2]);
+        field2index[ATTRIBUTES+nfield+nthresh] = add_compute(aname);
+        break;
 
-      char *ptr = strchr(suffix,'[');
-      if (ptr) {
-        if (suffix[strlen(suffix)-1] != ']')
-          error->all(FLERR,"Invalid attribute in dump modify command");
-        argindex[ATTRIBUTES+nfield+nthresh] = atoi(ptr+1);
-        *ptr = '\0';
-      } else argindex[ATTRIBUTES+nfield+nthresh] = 0;
+      // fix value = f_ID
+      // if no trailing [], then arg is set to 0, else arg is between []
 
-      n = modify->find_fix(suffix);
-      if (n < 0) error->all(FLERR,"Could not find dump modify fix ID");
+      case ArgInfo::FIX:
+        thresh_array[nthresh] = FIX;
+        n = modify->find_fix(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump modify fix ID: {}",aname);
 
-      if (modify->fix[n]->peratom_flag == 0)
-        error->all(FLERR,"Dump modify fix ID does not compute per-atom info");
-      if (argindex[ATTRIBUTES+nfield+nthresh] == 0 &&
-          modify->fix[n]->size_peratom_cols > 0)
-        error->all(FLERR,"Dump modify fix ID does not compute per-atom vector");
-      if (argindex[ATTRIBUTES+nfield+nthresh] > 0 &&
-          modify->fix[n]->size_peratom_cols == 0)
-        error->all(FLERR,"Dump modify fix ID does not compute per-atom array");
-      if (argindex[ATTRIBUTES+nfield+nthresh] > 0 &&
-          argindex[ATTRIBUTES+nfield+nthresh] > modify->fix[n]->size_peratom_cols)
-        error->all(FLERR,"Dump modify fix ID vector is not large enough");
+        if (modify->fix[n]->peratom_flag == 0)
+          error->all(FLERR,"Dump modify fix ID {} does not compute per-atom info",aname);
+        if (argi.get_dim() == 0 && modify->fix[n]->size_peratom_cols > 0)
+          error->all(FLERR,"Dump modify fix ID {} does not compute per-atom vector",aname);
+        if (argi.get_index1() > 0 && modify->fix[n]->size_peratom_cols == 0)
+          error->all(FLERR,"Dump modify fix ID {} does not compute per-atom array",aname);
+        if (argi.get_index1() > 0 && argi.get_index1() > modify->fix[n]->size_peratom_cols)
+          error->all(FLERR,"Dump modify fix ID {} vector is not large enough",aname);
 
-      field2index[ATTRIBUTES+nfield+nthresh] = add_fix(suffix);
-      delete [] suffix;
+        field2index[ATTRIBUTES+nfield+nthresh] = add_fix(aname);
+        break;
 
-    // variable value = v_ID
+      // variable value = v_ID
 
-    } else if (strncmp(arg[1],"v_",2) == 0) {
-      thresh_array[nthresh] = VARIABLE;
-      int n = strlen(arg[1]);
-      char *suffix = new char[n];
-      strcpy(suffix,&arg[1][2]);
+      case ArgInfo::VARIABLE:
+        thresh_array[nthresh] = VARIABLE;
+        n = input->variable->find(aname);
+        if (n < 0) error->all(FLERR,"Could not find dump modify variable name: {}",aname);
+        if (input->variable->atomstyle(n) == 0)
+          error->all(FLERR,"Dump modify variable {} is not atom-style variable",aname);
 
-      argindex[ATTRIBUTES+nfield+nthresh] = 0;
+        field2index[ATTRIBUTES+nfield+nthresh] = add_variable(aname);
+        break;
 
-      n = input->variable->find(suffix);
-      if (n < 0) error->all(FLERR,"Could not find dump modify variable name");
-      if (input->variable->atomstyle(n) == 0)
-        error->all(FLERR,"Dump modify variable is not atom-style variable");
+      // custom per atom floating point vector or array
 
-      field2index[ATTRIBUTES+nfield+nthresh] = add_variable(suffix);
-      delete [] suffix;
+      case ArgInfo::DNAME:
+        n = atom->find_custom(aname,flag,cols);
 
-    } else error->all(FLERR,"Invalid dump_modify threshold operator");
+        if (n < 0)
+          error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
+        if (argindex[ATTRIBUTES+nfield+nthresh] == 0) {
+          if (!flag || cols)
+            error->all(FLERR,"Property double vector for dump custom does not exist");
+          thresh_array[nthresh] = DVEC;
+        } else {
+          if (!flag || !cols)
+            error->all(FLERR,"Property double array for dump custom does not exist");
+          if (argindex[ATTRIBUTES+nfield+nthresh] > atom->dcols[n])
+            error->all(FLERR,"Dump custom property array is accessed out-of-range");
+          thresh_array[nthresh] = DARRAY;
+        }
+
+        field2index[ATTRIBUTES+nfield+nthresh] = add_custom(aname,thresh_array[nthresh]);
+        break;
+
+      // custom per atom integer vector or array
+
+      case ArgInfo::INAME:
+        n = atom->find_custom(aname,flag,cols);
+
+        if (n < 0)
+          error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
+        if (argindex[ATTRIBUTES+nfield+nthresh] == 0) {
+          if (flag || cols)
+            error->all(FLERR,"Property integer vector for dump custom does not exist");
+          thresh_array[nthresh] = IVEC;
+        } else {
+          if (flag || !cols)
+            error->all(FLERR,"Property integer array for dump custom does not exist");
+          if (argindex[ATTRIBUTES+nfield+nthresh] > atom->icols[n])
+            error->all(FLERR,"Dump custom property array is accessed out-of-range");
+          thresh_array[nthresh] = IARRAY;
+        }
+
+        field2index[ATTRIBUTES+nfield+nthresh] = add_custom(aname,thresh_array[nthresh]);
+        break;
+
+      // no match
+
+      default:
+        error->all(FLERR,"Invalid dump_modify thresh attribute: {}",aname);
+        break;
+      }
+    }
 
     // set operation type of threshold
 
@@ -2298,11 +2377,32 @@ int DumpVTK::modify_param(int narg, char **arg)
     else if (strcmp(arg[2],">=") == 0) thresh_op[nthresh] = GE;
     else if (strcmp(arg[2],"==") == 0) thresh_op[nthresh] = EQ;
     else if (strcmp(arg[2],"!=") == 0) thresh_op[nthresh] = NEQ;
-    else error->all(FLERR,"Invalid dump_modify threshold operator");
+    else if (strcmp(arg[2],"|^") == 0) thresh_op[nthresh] = XOR;
+    else error->all(FLERR,"Invalid dump_modify thresh operator");
 
-    // set threshold value
+    // set threshold value as number or special LAST keyword
+    // create FixStore to hold LAST values, should work with restart
+    // id = dump-ID + nthreshlast + DUMP_STORE, fix group = dump group
 
-    thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp);
+    if (strcmp(arg[3],"LAST") != 0) {
+      thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp);
+      thresh_last[nthresh] = -1;
+    } else {
+      thresh_fix = (FixStore **)
+        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStore *),"dump:thresh_fix");
+      thresh_fixID = (char **)
+        memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),"dump:thresh_fixID");
+      memory->grow(thresh_first,(nthreshlast+1),"dump:thresh_first");
+
+      std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast);
+      thresh_fixID[nthreshlast] = utils::strdup(threshid);
+      threshid += fmt::format(" {} STORE peratom 1 1", group->names[igroup]);
+      thresh_fix[nthreshlast] = (FixStore *) modify->add_fix(threshid);
+
+      thresh_last[nthreshlast] = nthreshlast;
+      thresh_first[nthreshlast] = 1;
+      nthreshlast++;
+    }
 
     nthresh++;
     return 4;
@@ -2387,25 +2487,35 @@ void DumpVTK::pack_variable(int n)
 
 void DumpVTK::pack_custom(int n)
 {
-  int index = field2index[n];
-
-  if (flag_custom[index] == 0) { // integer
-    int iwhich,tmp;
-    iwhich = atom->find_custom(id_custom[index],tmp);
+  int flag = custom_flag[field2index[current_pack_choice_key]];
+  int iwhich = custom[field2index[current_pack_choice_key]];
+  int index = argindex[current_pack_choice_key];
 
+  if (flag == IVEC) {
     int *ivector = atom->ivector[iwhich];
     for (int i = 0; i < nchoose; i++) {
       buf[n] = ivector[clist[i]];
       n += size_one;
     }
-  } else if (flag_custom[index] == 1) { // double
-    int iwhich,tmp;
-    iwhich = atom->find_custom(id_custom[index],tmp);
-
+  } else if (flag == DVEC) {
     double *dvector = atom->dvector[iwhich];
     for (int i = 0; i < nchoose; i++) {
       buf[n] = dvector[clist[i]];
       n += size_one;
     }
+  } else if (flag == IARRAY) {
+    index--;
+    int **iarray = atom->iarray[iwhich];
+    for (int i = 0; i < nchoose; i++) {
+      buf[n] = iarray[clist[i]][index];
+      n += size_one;
+    }
+  } else if (flag == DARRAY) {
+    index--;
+    double **darray = atom->darray[iwhich];
+    for (int i = 0; i < nchoose; i++) {
+      buf[n] = darray[clist[i]][index];
+      n += size_one;
+    }
   }
 }