From fcaaf75764f752ad9bcc824b02431543b40697c9 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 15:00:31 +0000
Subject: [PATCH 01/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12727
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/variable.cpp | 895 ++++++++++++++++++++++++++---------------------
 src/variable.h   |   7 +-
 2 files changed, 503 insertions(+), 399 deletions(-)

diff --git a/src/variable.cpp b/src/variable.cpp
index 9c8c04d0dc..551acabd1f 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -45,6 +45,7 @@ using namespace MathConst;
 #define MAXLINE 256
 #define CHUNK 1024
 #define VALUELENGTH 64
+#define MAXFUNCARG 6
 
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
@@ -59,7 +60,7 @@ enum{ARG,OP};
 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
      NOT,EQ,NE,LT,LE,GT,GE,AND,OR,
      SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
-     RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE,
+     RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE,STRIDE2,
      VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
      VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY};
 
@@ -973,7 +974,8 @@ double Variable::evaluate(char *str, Tree **tree)
         Tree *newtree = new Tree();
         newtree->type = VALUE;
         newtree->value = atof(number);
-        newtree->left = newtree->middle = newtree->right = NULL;
+        newtree->first = newtree->second = NULL;
+        newtree->extra = 0;
         treestack[ntreestack++] = newtree;
       } else argstack[nargstack++] = atof(number);
 
@@ -1057,7 +1059,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1082,7 +1085,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1110,7 +1114,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1179,7 +1184,8 @@ double Variable::evaluate(char *str, Tree **tree)
           newtree->array = compute->vector_atom;
           newtree->nstride = 1;
           newtree->selfalloc = 0;
-          newtree->left = newtree->middle = newtree->right = NULL;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
         // c_ID[i] = vector from per-atom array
@@ -1210,7 +1216,8 @@ double Variable::evaluate(char *str, Tree **tree)
             newtree->array = NULL;
           newtree->nstride = compute->size_peratom_cols;
           newtree->selfalloc = 0;
-          newtree->left = newtree->middle = newtree->right = NULL;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
         } else error->all(FLERR,"Mismatched compute in variable formula");
@@ -1265,7 +1272,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1284,7 +1292,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1306,7 +1315,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1360,7 +1370,8 @@ double Variable::evaluate(char *str, Tree **tree)
           newtree->array = fix->vector_atom;
           newtree->nstride = 1;
           newtree->selfalloc = 0;
-          newtree->left = newtree->middle = newtree->right = NULL;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
         // f_ID[i] = vector from per-atom array
@@ -1385,7 +1396,8 @@ double Variable::evaluate(char *str, Tree **tree)
             newtree->array = NULL;
           newtree->nstride = fix->size_peratom_cols;
           newtree->selfalloc = 0;
-          newtree->left = newtree->middle = newtree->right = NULL;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
         } else error->all(FLERR,"Mismatched fix in variable formula");
@@ -1430,7 +1442,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = atof(var);
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = atof(var);
 
@@ -1458,7 +1471,8 @@ double Variable::evaluate(char *str, Tree **tree)
           newtree->array = reader[ivar]->fix->vstore;
           newtree->nstride = 1;
           newtree->selfalloc = 0;
-          newtree->left = newtree->middle = newtree->right = NULL;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
         // v_name[N] = scalar from atom-style variable
@@ -1548,7 +1562,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
@@ -1568,7 +1583,8 @@ double Variable::evaluate(char *str, Tree **tree)
             Tree *newtree = new Tree();
             newtree->type = VALUE;
             newtree->value = value1;
-            newtree->left = newtree->middle = newtree->right = NULL;
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
         }
@@ -1645,12 +1661,13 @@ double Variable::evaluate(char *str, Tree **tree)
           Tree *newtree = new Tree();
           newtree->type = opprevious;
           if (opprevious == UNARY) {
-            newtree->left = treestack[--ntreestack];
-            newtree->middle = newtree->right = NULL;
+            newtree->first = treestack[--ntreestack];
+            newtree->second = NULL;
+            newtree->nextra = 0;
           } else {
-            newtree->right = treestack[--ntreestack];
-            newtree->middle = NULL;
-            newtree->left = treestack[--ntreestack];
+            newtree->second = treestack[--ntreestack];
+            newtree->first = treestack[--ntreestack];
+            newtree->nextra = 0;
           }
           treestack[ntreestack++] = newtree;
 
@@ -1762,36 +1779,36 @@ double Variable::collapse_tree(Tree *tree)
   if (tree->type == BIGINTARRAY) return 0.0;
 
   if (tree->type == ADD) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = arg1 + arg2;
     return tree->value;
   }
 
   if (tree->type == SUBTRACT) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = arg1 - arg2;
     return tree->value;
   }
 
   if (tree->type == MULTIPLY) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = arg1 * arg2;
     return tree->value;
   }
 
   if (tree->type == DIVIDE) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg2 == 0.0) error->one(FLERR,"Divide by 0 in variable formula");
     tree->value = arg1 / arg2;
@@ -1799,9 +1816,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == MODULO) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg2 == 0.0) error->one(FLERR,"Modulo 0 in variable formula");
     tree->value = fmod(arg1,arg2);
@@ -1809,9 +1826,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == CARAT) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg2 == 0.0) error->one(FLERR,"Power by 0 in variable formula");
     tree->value = pow(arg1,arg2);
@@ -1819,16 +1836,16 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == UNARY) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = -arg1;
     return tree->value;
   }
 
   if (tree->type == NOT) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 == 0.0) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1836,9 +1853,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == EQ) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 == arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1846,9 +1863,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == NE) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 != arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1856,9 +1873,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == LT) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 < arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1866,9 +1883,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == LE) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 <= arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1876,9 +1893,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == GT) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 > arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1886,9 +1903,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == GE) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 >= arg2) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1896,9 +1913,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == AND) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 != 0.0 && arg2 != 0.0) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1906,9 +1923,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == OR) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 != 0.0 || arg2 != 0.0) tree->value = 1.0;
     else tree->value = 0.0;
@@ -1916,8 +1933,8 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == SQRT) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 < 0.0)
       error->one(FLERR,"Sqrt of negative value in variable formula");
@@ -1926,16 +1943,16 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == EXP) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = exp(arg1);
     return tree->value;
   }
 
   if (tree->type == LN) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 <= 0.0)
       error->one(FLERR,"Log of zero/negative value in variable formula");
@@ -1944,8 +1961,8 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == LOG) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 <= 0.0)
       error->one(FLERR,"Log of zero/negative value in variable formula");
@@ -1954,40 +1971,40 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == ABS) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = fabs(arg1);
     return tree->value;
   }
 
   if (tree->type == SIN) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = sin(arg1);
     return tree->value;
   }
 
   if (tree->type == COS) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = cos(arg1);
     return tree->value;
   }
 
   if (tree->type == TAN) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = tan(arg1);
     return tree->value;
   }
 
   if (tree->type == ASIN) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 < -1.0 || arg1 > 1.0)
       error->one(FLERR,"Arcsin of invalid value in variable formula");
@@ -1996,8 +2013,8 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == ACOS) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg1 < -1.0 || arg1 > 1.0)
       error->one(FLERR,"Arccos of invalid value in variable formula");
@@ -2006,17 +2023,17 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == ATAN) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = atan(arg1);
     return tree->value;
   }
 
   if (tree->type == ATAN2) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = atan2(arg1,arg2);
     return tree->value;
@@ -2025,10 +2042,10 @@ double Variable::collapse_tree(Tree *tree)
   // random() or normal() do not become a single collapsed value
 
   if (tree->type == RANDOM) {
-    collapse_tree(tree->left);
-    collapse_tree(tree->middle);
+    collapse_tree(tree->first);
+    collapse_tree(tree->second);
     if (randomatom == NULL) {
-      int seed = static_cast<int> (collapse_tree(tree->right));
+      int seed = static_cast<int> (collapse_tree(tree->extra[0]));
       if (seed <= 0)
         error->one(FLERR,"Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
@@ -2037,12 +2054,12 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == NORMAL) {
-    collapse_tree(tree->left);
-    double sigma = collapse_tree(tree->middle);
+    collapse_tree(tree->first);
+    double sigma = collapse_tree(tree->second);
     if (sigma < 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
     if (randomatom == NULL) {
-      int seed = static_cast<int> (collapse_tree(tree->right));
+      int seed = static_cast<int> (collapse_tree(tree->extra[0]));
       if (seed <= 0)
         error->one(FLERR,"Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
@@ -2051,33 +2068,33 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == CEIL) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = ceil(arg1);
     return tree->value;
   }
 
   if (tree->type == FLOOR) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = floor(arg1);
     return tree->value;
   }
 
   if (tree->type == ROUND) {
-    arg1 = collapse_tree(tree->left);
-    if (tree->left->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    if (tree->first->type != VALUE) return 0.0;
     tree->type = VALUE;
     tree->value = MYROUND(arg1);
     return tree->value;
   }
 
   if (tree->type == RAMP) {
-    arg1 = collapse_tree(tree->left);
-    arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     double delta = update->ntimestep - update->beginstep;
     if (delta != 0.0) delta /= update->endstep - update->beginstep;
@@ -2086,9 +2103,9 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == STAGGER) {
-    int ivalue1 = static_cast<int> (collapse_tree(tree->left));
-    int ivalue2 = static_cast<int> (collapse_tree(tree->right));
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    int ivalue1 = static_cast<int> (collapse_tree(tree->first));
+    int ivalue2 = static_cast<int> (collapse_tree(tree->second));
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
       error->one(FLERR,"Invalid math function in variable formula");
@@ -2100,11 +2117,11 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == LOGFREQ) {
-    int ivalue1 = static_cast<int> (collapse_tree(tree->left));
-    int ivalue2 = static_cast<int> (collapse_tree(tree->middle));
-    int ivalue3 = static_cast<int> (collapse_tree(tree->right));
-    if (tree->left->type != VALUE || tree->middle->type != VALUE ||
-        tree->right->type != VALUE) return 0.0;
+    int ivalue1 = static_cast<int> (collapse_tree(tree->first));
+    int ivalue2 = static_cast<int> (collapse_tree(tree->second));
+    int ivalue3 = static_cast<int> (collapse_tree(tree->extra[0]));
+    if (tree->first->type != VALUE || tree->second->type != VALUE ||
+        tree->extra[0]->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
       error->one(FLERR,"Invalid math function in variable formula");
@@ -2120,11 +2137,11 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == STRIDE) {
-    int ivalue1 = static_cast<int> (collapse_tree(tree->left));
-    int ivalue2 = static_cast<int> (collapse_tree(tree->middle));
-    int ivalue3 = static_cast<int> (collapse_tree(tree->right));
-    if (tree->left->type != VALUE || tree->middle->type != VALUE ||
-        tree->right->type != VALUE) return 0.0;
+    int ivalue1 = static_cast<int> (collapse_tree(tree->first));
+    int ivalue2 = static_cast<int> (collapse_tree(tree->second));
+    int ivalue3 = static_cast<int> (collapse_tree(tree->extra[0]));
+    if (tree->first->type != VALUE || tree->second->type != VALUE ||
+        tree->extra[0]->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
       error->one(FLERR,"Invalid math function in variable formula");
@@ -2137,10 +2154,50 @@ double Variable::collapse_tree(Tree *tree)
     return tree->value;
   }
 
+  if (tree->type == STRIDE2) {
+    int ivalue1 = static_cast<int> (collapse_tree(tree->first));
+    int ivalue2 = static_cast<int> (collapse_tree(tree->second));
+    int ivalue3 = static_cast<int> (collapse_tree(tree->extra[0]));
+    int ivalue4 = static_cast<int> (collapse_tree(tree->extra[1]));
+    int ivalue5 = static_cast<int> (collapse_tree(tree->extra[2]));
+    int ivalue6 = static_cast<int> (collapse_tree(tree->extra[3]));
+    if (tree->first->type != VALUE || tree->second->type != VALUE ||
+        tree->extra[0]->type != VALUE || tree->extra[1]->type != VALUE || 
+	tree->extra[2]->type != VALUE || tree->extra[3]->type != VALUE) 
+      return 0.0;
+    tree->type = VALUE;
+    if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
+      error->one(FLERR,"Invalid math function in variable formula");
+    if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
+      error->one(FLERR,"Invalid math function in variable formula");
+    if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
+      error->one(FLERR,"Invalid math function in variable formula");
+    bigint istep;
+    if (update->ntimestep < ivalue1) istep = ivalue1;
+    else if (update->ntimestep < ivalue2) {
+      if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
+        int offset = update->ntimestep - ivalue1;
+        istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+        if (update->ntimestep < ivalue2 && istep > ivalue4) 
+          tree->value = ivalue4;
+      } else {
+        int offset = update->ntimestep - ivalue4;
+        istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
+        if (istep > ivalue5) {
+          int offset = ivalue5 - ivalue1;
+          istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+          if (istep > ivalue2) istep = 9.0e18;
+        }
+      }
+    } else istep = 9.0e18;
+    tree->value = istep;
+    return tree->value;
+  }
+
   if (tree->type == VDISPLACE) {
-    double arg1 = collapse_tree(tree->left);
-    double arg2 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+    double arg1 = collapse_tree(tree->first);
+    double arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
     tree->type = VALUE;
     double delta = update->ntimestep - update->beginstep;
     tree->value = arg1 + arg2*delta*update->dt;
@@ -2148,11 +2205,11 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == SWIGGLE) {
-    double arg1 = collapse_tree(tree->left);
-    double arg2 = collapse_tree(tree->middle);
-    double arg3 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->middle->type != VALUE ||
-        tree->right->type != VALUE) return 0.0;
+    double arg1 = collapse_tree(tree->first);
+    double arg2 = collapse_tree(tree->second);
+    double arg3 = collapse_tree(tree->extra[0]);
+    if (tree->first->type != VALUE || tree->second->type != VALUE ||
+        tree->extra[0]->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg3 == 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
@@ -2163,11 +2220,11 @@ double Variable::collapse_tree(Tree *tree)
   }
 
   if (tree->type == CWIGGLE) {
-    double arg1 = collapse_tree(tree->left);
-    double arg2 = collapse_tree(tree->middle);
-    double arg3 = collapse_tree(tree->right);
-    if (tree->left->type != VALUE || tree->middle->type != VALUE ||
-        tree->right->type != VALUE) return 0.0;
+    double arg1 = collapse_tree(tree->first);
+    double arg2 = collapse_tree(tree->second);
+    double arg3 = collapse_tree(tree->extra[0]);
+    if (tree->first->type != VALUE || tree->second->type != VALUE ||
+        tree->extra[0]->type != VALUE) return 0.0;
     tree->type = VALUE;
     if (arg3 == 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
@@ -2208,119 +2265,119 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == BIGINTARRAY) return (double) tree->barray[i*tree->nstride];
 
   if (tree->type == ADD)
-    return eval_tree(tree->left,i) + eval_tree(tree->right,i);
+    return eval_tree(tree->first,i) + eval_tree(tree->second,i);
   if (tree->type == SUBTRACT)
-    return eval_tree(tree->left,i) - eval_tree(tree->right,i);
+    return eval_tree(tree->first,i) - eval_tree(tree->second,i);
   if (tree->type == MULTIPLY)
-    return eval_tree(tree->left,i) * eval_tree(tree->right,i);
+    return eval_tree(tree->first,i) * eval_tree(tree->second,i);
   if (tree->type == DIVIDE) {
-    double denom = eval_tree(tree->right,i);
+    double denom = eval_tree(tree->second,i);
     if (denom == 0.0) error->one(FLERR,"Divide by 0 in variable formula");
-    return eval_tree(tree->left,i) / denom;
+    return eval_tree(tree->first,i) / denom;
   }
   if (tree->type == MODULO) {
-    double denom = eval_tree(tree->right,i);
+    double denom = eval_tree(tree->second,i);
     if (denom == 0.0) error->one(FLERR,"Modulo 0 in variable formula");
-    return fmod(eval_tree(tree->left,i),denom);
+    return fmod(eval_tree(tree->first,i),denom);
   }
   if (tree->type == CARAT) {
-    double exponent = eval_tree(tree->right,i);
+    double exponent = eval_tree(tree->second,i);
     if (exponent == 0.0) error->one(FLERR,"Power by 0 in variable formula");
-    return pow(eval_tree(tree->left,i),exponent);
+    return pow(eval_tree(tree->first,i),exponent);
   }
-  if (tree->type == UNARY) return -eval_tree(tree->left,i);
+  if (tree->type == UNARY) return -eval_tree(tree->first,i);
 
   if (tree->type == NOT) {
-    if (eval_tree(tree->left,i) == 0.0) return 1.0;
+    if (eval_tree(tree->first,i) == 0.0) return 1.0;
     else return 0.0;
   }
   if (tree->type == EQ) {
-    if (eval_tree(tree->left,i) == eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) == eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == NE) {
-    if (eval_tree(tree->left,i) != eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) != eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == LT) {
-    if (eval_tree(tree->left,i) < eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) < eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == LE) {
-    if (eval_tree(tree->left,i) <= eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) <= eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == GT) {
-    if (eval_tree(tree->left,i) > eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) > eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == GE) {
-    if (eval_tree(tree->left,i) >= eval_tree(tree->right,i)) return 1.0;
+    if (eval_tree(tree->first,i) >= eval_tree(tree->second,i)) return 1.0;
     else return 0.0;
   }
   if (tree->type == AND) {
-    if (eval_tree(tree->left,i) != 0.0 && eval_tree(tree->right,i) != 0.0)
+    if (eval_tree(tree->first,i) != 0.0 && eval_tree(tree->second,i) != 0.0)
       return 1.0;
     else return 0.0;
   }
   if (tree->type == OR) {
-    if (eval_tree(tree->left,i) != 0.0 || eval_tree(tree->right,i) != 0.0)
+    if (eval_tree(tree->first,i) != 0.0 || eval_tree(tree->second,i) != 0.0)
       return 1.0;
     else return 0.0;
   }
 
   if (tree->type == SQRT) {
-    arg1 = eval_tree(tree->left,i);
+    arg1 = eval_tree(tree->first,i);
     if (arg1 < 0.0)
       error->one(FLERR,"Sqrt of negative value in variable formula");
     return sqrt(arg1);
   }
   if (tree->type == EXP)
-    return exp(eval_tree(tree->left,i));
+    return exp(eval_tree(tree->first,i));
   if (tree->type == LN) {
-    arg1 = eval_tree(tree->left,i);
+    arg1 = eval_tree(tree->first,i);
     if (arg1 <= 0.0)
       error->one(FLERR,"Log of zero/negative value in variable formula");
     return log(arg1);
   }
   if (tree->type == LOG) {
-    arg1 = eval_tree(tree->left,i);
+    arg1 = eval_tree(tree->first,i);
     if (arg1 <= 0.0)
       error->one(FLERR,"Log of zero/negative value in variable formula");
     return log10(arg1);
   }
   if (tree->type == ABS)
-    return fabs(eval_tree(tree->left,i));
+    return fabs(eval_tree(tree->first,i));
 
   if (tree->type == SIN)
-    return sin(eval_tree(tree->left,i));
+    return sin(eval_tree(tree->first,i));
   if (tree->type == COS)
-    return cos(eval_tree(tree->left,i));
+    return cos(eval_tree(tree->first,i));
   if (tree->type == TAN)
-    return tan(eval_tree(tree->left,i));
+    return tan(eval_tree(tree->first,i));
 
   if (tree->type == ASIN) {
-    arg1 = eval_tree(tree->left,i);
+    arg1 = eval_tree(tree->first,i);
     if (arg1 < -1.0 || arg1 > 1.0)
       error->one(FLERR,"Arcsin of invalid value in variable formula");
     return asin(arg1);
   }
   if (tree->type == ACOS) {
-    arg1 = eval_tree(tree->left,i);
+    arg1 = eval_tree(tree->first,i);
     if (arg1 < -1.0 || arg1 > 1.0)
       error->one(FLERR,"Arccos of invalid value in variable formula");
     return acos(arg1);
   }
   if (tree->type == ATAN)
-    return atan(eval_tree(tree->left,i));
+    return atan(eval_tree(tree->first,i));
   if (tree->type == ATAN2)
-    return atan2(eval_tree(tree->left,i),eval_tree(tree->right,i));
+    return atan2(eval_tree(tree->first,i),eval_tree(tree->second,i));
 
   if (tree->type == RANDOM) {
-    double lower = eval_tree(tree->left,i);
-    double upper = eval_tree(tree->middle,i);
+    double lower = eval_tree(tree->first,i);
+    double upper = eval_tree(tree->second,i);
     if (randomatom == NULL) {
-      int seed = static_cast<int> (eval_tree(tree->right,i));
+      int seed = static_cast<int> (eval_tree(tree->extra[0],i));
       if (seed <= 0)
         error->one(FLERR,"Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
@@ -2328,12 +2385,12 @@ double Variable::eval_tree(Tree *tree, int i)
     return randomatom->uniform()*(upper-lower)+lower;
   }
   if (tree->type == NORMAL) {
-    double mu = eval_tree(tree->left,i);
-    double sigma = eval_tree(tree->middle,i);
+    double mu = eval_tree(tree->first,i);
+    double sigma = eval_tree(tree->second,i);
     if (sigma < 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
     if (randomatom == NULL) {
-      int seed = static_cast<int> (eval_tree(tree->right,i));
+      int seed = static_cast<int> (eval_tree(tree->extra[0],i));
       if (seed <= 0)
         error->one(FLERR,"Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
@@ -2342,15 +2399,15 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == CEIL)
-    return ceil(eval_tree(tree->left,i));
+    return ceil(eval_tree(tree->first,i));
   if (tree->type == FLOOR)
-    return floor(eval_tree(tree->left,i));
+    return floor(eval_tree(tree->first,i));
   if (tree->type == ROUND)
-    return MYROUND(eval_tree(tree->left,i));
+    return MYROUND(eval_tree(tree->first,i));
 
   if (tree->type == RAMP) {
-    arg1 = eval_tree(tree->left,i);
-    arg2 = eval_tree(tree->right,i);
+    arg1 = eval_tree(tree->first,i);
+    arg2 = eval_tree(tree->second,i);
     double delta = update->ntimestep - update->beginstep;
     if (delta != 0.0) delta /= update->endstep - update->beginstep;
     arg = arg1 + delta*(arg2-arg1);
@@ -2358,8 +2415,8 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == STAGGER) {
-    int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
-    int ivalue2 = static_cast<int> (eval_tree(tree->right,i));
+    int ivalue1 = static_cast<int> (eval_tree(tree->first,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->second,i));
     if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
       error->one(FLERR,"Invalid math function in variable formula");
     int lower = update->ntimestep/ivalue1 * ivalue1;
@@ -2370,9 +2427,9 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == LOGFREQ) {
-    int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
-    int ivalue2 = static_cast<int> (eval_tree(tree->middle,i));
-    int ivalue3 = static_cast<int> (eval_tree(tree->right,i));
+    int ivalue1 = static_cast<int> (eval_tree(tree->first,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->second,i));
+    int ivalue3 = static_cast<int> (eval_tree(tree->extra[0],i));
     if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
       error->one(FLERR,"Invalid math function in variable formula");
     if (update->ntimestep < ivalue1) arg = ivalue1;
@@ -2387,9 +2444,9 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == STRIDE) {
-    int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
-    int ivalue2 = static_cast<int> (eval_tree(tree->middle,i));
-    int ivalue3 = static_cast<int> (eval_tree(tree->right,i));
+    int ivalue1 = static_cast<int> (eval_tree(tree->first,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->second,i));
+    int ivalue3 = static_cast<int> (eval_tree(tree->extra[0],i));
     if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
       error->one(FLERR,"Invalid math function in variable formula");
     if (update->ntimestep < ivalue1) arg = ivalue1;
@@ -2401,18 +2458,53 @@ double Variable::eval_tree(Tree *tree, int i)
     return arg;
   }
 
+  if (tree->type == STRIDE2) {
+    int ivalue1 = static_cast<int> (eval_tree(tree->first,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->second,i));
+    int ivalue3 = static_cast<int> (eval_tree(tree->extra[0],i));
+    int ivalue4 = static_cast<int> (eval_tree(tree->extra[1],i));
+    int ivalue5 = static_cast<int> (eval_tree(tree->extra[2],i));
+    int ivalue6 = static_cast<int> (eval_tree(tree->extra[3],i));
+    if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
+      error->one(FLERR,"Invalid math function in variable formula");
+    if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
+      error->one(FLERR,"Invalid math function in variable formula");
+    if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
+      error->one(FLERR,"Invalid math function in variable formula");
+    bigint istep;
+    if (update->ntimestep < ivalue1) istep = ivalue1;
+    else if (update->ntimestep < ivalue2) {
+      if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
+        int offset = update->ntimestep - ivalue1;
+        istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+        if (update->ntimestep < ivalue2 && istep > ivalue4) 
+          tree->value = ivalue4;
+      } else {
+        int offset = update->ntimestep - ivalue4;
+        istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
+        if (istep > ivalue5) {
+          int offset = ivalue5 - ivalue1;
+          istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+          if (istep > ivalue2) istep = 9.0e18;
+        }
+      }
+    } else istep = 9.0e18;
+    arg = istep;
+    return arg;
+  }
+
   if (tree->type == VDISPLACE) {
-    arg1 = eval_tree(tree->left,i);
-    arg2 = eval_tree(tree->right,i);
+    arg1 = eval_tree(tree->first,i);
+    arg2 = eval_tree(tree->second,i);
     double delta = update->ntimestep - update->beginstep;
     arg = arg1 + arg2*delta*update->dt;
     return arg;
   }
 
   if (tree->type == SWIGGLE) {
-    arg1 = eval_tree(tree->left,i);
-    arg2 = eval_tree(tree->middle,i);
-    arg3 = eval_tree(tree->right,i);
+    arg1 = eval_tree(tree->first,i);
+    arg2 = eval_tree(tree->second,i);
+    arg3 = eval_tree(tree->extra[0],i);
     if (arg3 == 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
@@ -2422,9 +2514,9 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == CWIGGLE) {
-    arg1 = eval_tree(tree->left,i);
-    arg2 = eval_tree(tree->middle,i);
-    arg3 = eval_tree(tree->right,i);
+    arg1 = eval_tree(tree->first,i);
+    arg2 = eval_tree(tree->second,i);
+    arg3 = eval_tree(tree->extra[0],i);
     if (arg3 == 0.0)
       error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
@@ -2460,9 +2552,12 @@ double Variable::eval_tree(Tree *tree, int i)
 
 void Variable::free_tree(Tree *tree)
 {
-  if (tree->left) free_tree(tree->left);
-  if (tree->middle) free_tree(tree->middle);
-  if (tree->right) free_tree(tree->right);
+  if (tree->first) free_tree(tree->first);
+  if (tree->second) free_tree(tree->second);
+  if (tree->nextra) {
+    for (int i = 0; i < tree->nextra; i++) free_tree(tree->extra[i]);
+    delete [] tree->extra;
+  }
 
   if (tree->type == ATOMARRAY && tree->selfalloc)
     memory->destroy(tree->array);
@@ -2568,12 +2663,12 @@ int Variable::int_between_brackets(char *&ptr, int varallow)
    process a math function in formula
    push result onto tree or arg stack
    word = math function
-   contents = str between parentheses with one,two,three args
+   contents = str between parentheses with comma-separated args
    return 0 if not a match, 1 if successfully processed
    customize by adding a math function:
      sqrt(),exp(),ln(),log(),abs(),sin(),cos(),tan(),asin(),acos(),atan(),
      atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
-     ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z),
+     ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),
      vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z)
 ------------------------------------------------------------------------- */
 
@@ -2594,80 +2689,55 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       strcmp(word,"floor") && strcmp(word,"round") &&
       strcmp(word,"ramp") && strcmp(word,"stagger") &&
       strcmp(word,"logfreq") && strcmp(word,"stride") &&
-      strcmp(word,"vdisplace") &&
+      strcmp(word,"stride2") && strcmp(word,"vdisplace") &&
       strcmp(word,"swiggle") && strcmp(word,"cwiggle"))
     return 0;
 
-  // parse contents for arg1,arg2,arg3 separated by commas
-  // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
-
-  char *arg1,*arg2,*arg3;
-  char *ptr1,*ptr2;
-
-  ptr1 = find_next_comma(contents);
-  if (ptr1) {
-    *ptr1 = '\0';
-    ptr2 = find_next_comma(ptr1+1);
-    if (ptr2) *ptr2 = '\0';
-  } else ptr2 = NULL;
-
-  int n = strlen(contents) + 1;
-  arg1 = new char[n];
-  strcpy(arg1,contents);
-  int narg = 1;
-  if (ptr1) {
-    n = strlen(ptr1+1) + 1;
-    arg2 = new char[n];
-    strcpy(arg2,ptr1+1);
-    narg = 2;
-  } else arg2 = NULL;
-  if (ptr2) {
-    n = strlen(ptr2+1) + 1;
-    arg3 = new char[n];
-    strcpy(arg3,ptr2+1);
-    narg = 3;
-  } else arg3 = NULL;
-
-  // evaluate args
+  // parse contents for comma-separated args
+  // narg = number of args, args = strings between commas
+  
+  char *args[MAXFUNCARG];
+  int narg = parse_args(contents,args);
 
   Tree *newtree;
-  double value1,value2,value3;
+  double value1,value2;
+  double values[MAXFUNCARG-2];
 
   if (tree) {
     newtree = new Tree();
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     Tree *argtree;
-    if (narg == 1) {
-      evaluate(arg1,&argtree);
-      newtree->left = argtree;
-      newtree->middle = newtree->right = NULL;
-    } else if (narg == 2) {
-      evaluate(arg1,&argtree);
-      newtree->left = argtree;
-      newtree->middle = NULL;
-      evaluate(arg2,&argtree);
-      newtree->right = argtree;
-    } else if (narg == 3) {
-      evaluate(arg1,&argtree);
-      newtree->left = argtree;
-      evaluate(arg2,&argtree);
-      newtree->middle = argtree;
-      evaluate(arg3,&argtree);
-      newtree->right = argtree;
+    evaluate(args[0],&argtree);
+    newtree->first = argtree;
+    if (narg > 1) {
+      evaluate(args[1],&argtree);
+      newtree->second = argtree;
+      if (narg > 2) {
+        newtree->nextra = narg-2;
+        newtree->extra = new Tree*[narg-2];
+        for (int i = 2; i < narg; i++) {
+          evaluate(args[i],&argtree);
+          newtree->extra[i-2] = argtree;
+        }
+      }
     }
     treestack[ntreestack++] = newtree;
+
   } else {
-    if (narg == 1) {
-      value1 = evaluate(arg1,NULL);
-    } else if (narg == 2) {
-      value1 = evaluate(arg1,NULL);
-      value2 = evaluate(arg2,NULL);
-    } else if (narg == 3) {
-      value1 = evaluate(arg1,NULL);
-      value2 = evaluate(arg2,NULL);
-      value3 = evaluate(arg3,NULL);
+    value1 = evaluate(args[0],NULL);
+    if (narg > 1) {
+      value2 = evaluate(args[1],NULL);
+      if (narg > 2) {
+        for (int i = 2; i < narg; i++) 
+          values[i-2] = evaluate(args[i],NULL);
+      }
     }
   }
 
+  // individual math functions
+  // customize by adding a function
+
   if (strcmp(word,"sqrt") == 0) {
     if (narg != 1)
       error->all(FLERR,"Invalid math function in variable formula");
@@ -2758,7 +2828,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     if (tree) newtree->type = RANDOM;
     else {
       if (randomequal == NULL) {
-        int seed = static_cast<int> (value3);
+        int seed = static_cast<int> (values[0]);
         if (seed <= 0)
           error->all(FLERR,"Invalid math function in variable formula");
         randomequal = new RanMars(lmp,seed);
@@ -2773,7 +2843,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       if (value2 < 0.0)
         error->all(FLERR,"Invalid math function in variable formula");
       if (randomequal == NULL) {
-        int seed = static_cast<int> (value3);
+        int seed = static_cast<int> (values[0]);
         if (seed <= 0)
           error->all(FLERR,"Invalid math function in variable formula");
         randomequal = new RanMars(lmp,seed);
@@ -2836,7 +2906,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     else {
       int ivalue1 = static_cast<int> (value1);
       int ivalue2 = static_cast<int> (value2);
-      int ivalue3 = static_cast<int> (value3);
+      int ivalue3 = static_cast<int> (values[0]);
       if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
         error->all(FLERR,"Invalid math function in variable formula");
       double value;
@@ -2858,7 +2928,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     else {
       int ivalue1 = static_cast<int> (value1);
       int ivalue2 = static_cast<int> (value2);
-      int ivalue3 = static_cast<int> (value3);
+      int ivalue3 = static_cast<int> (values[0]);
       if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
         error->one(FLERR,"Invalid math function in variable formula");
       double value;
@@ -2871,6 +2941,44 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       argstack[nargstack++] = value;
     }
 
+  } else if (strcmp(word,"stride2") == 0) {
+    if (narg != 6)
+      error->all(FLERR,"Invalid math function in variable formula");
+    if (tree) newtree->type = STRIDE2;
+    else {
+      int ivalue1 = static_cast<int> (value1);
+      int ivalue2 = static_cast<int> (value2);
+      int ivalue3 = static_cast<int> (values[0]);
+      int ivalue4 = static_cast<int> (values[1]);
+      int ivalue5 = static_cast<int> (values[2]);
+      int ivalue6 = static_cast<int> (values[3]);
+      if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
+        error->one(FLERR,"Invalid math function in variable formula");
+      if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
+        error->one(FLERR,"Invalid math function in variable formula");
+      if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
+	error->one(FLERR,"Invalid math function in variable formula");
+      bigint istep;
+      if (update->ntimestep < ivalue1) istep = ivalue1;
+      else if (update->ntimestep < ivalue2) {
+        if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
+          int offset = update->ntimestep - ivalue1;
+          istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+          if (update->ntimestep < ivalue4 && istep > ivalue4) istep = ivalue4;
+        } else {
+          int offset = update->ntimestep - ivalue4;
+          istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
+          if (istep > ivalue5) {
+            int offset = ivalue5 - ivalue1;
+            istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
+            if (istep > ivalue2) istep = 9.0e18;
+          }
+        }
+      } else istep = 9.0e18;
+      double value = istep;
+      argstack[nargstack++] = value;
+    }
+
   } else if (strcmp(word,"vdisplace") == 0) {
     if (narg != 2)
       error->all(FLERR,"Invalid math function in variable formula");
@@ -2890,10 +2998,10 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       error->all(FLERR,"Cannot use swiggle in variable formula between runs");
     if (tree) newtree->type = CWIGGLE;
     else {
-      if (value3 == 0.0)
+      if (values[0] == 0.0)
         error->all(FLERR,"Invalid math function in variable formula");
       double delta = update->ntimestep - update->beginstep;
-      double omega = 2.0*MY_PI/value3;
+      double omega = 2.0*MY_PI/values[0];
       double value = value1 + value2*sin(omega*delta*update->dt);
       argstack[nargstack++] = value;
     }
@@ -2905,18 +3013,18 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       error->all(FLERR,"Cannot use cwiggle in variable formula between runs");
     if (tree) newtree->type = CWIGGLE;
     else {
-      if (value3 == 0.0)
+      if (values[0] == 0.0)
         error->all(FLERR,"Invalid math function in variable formula");
       double delta = update->ntimestep - update->beginstep;
-      double omega = 2.0*MY_PI/value3;
+      double omega = 2.0*MY_PI/values[0];
       double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
       argstack[nargstack++] = value;
     }
   }
 
-  delete [] arg1;
-  delete [] arg2;
-  delete [] arg3;
+  // delete stored args
+
+  for (int i = 0; i < narg; i++) delete [] args[i];
 
   return 1;
 }
@@ -2949,39 +3057,15 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       strcmp(word,"omega"))
     return 0;
 
-  // parse contents for arg1,arg2,arg3 separated by commas
-  // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
-
-  char *arg1,*arg2,*arg3;
-  char *ptr1,*ptr2;
-
-  ptr1 = find_next_comma(contents);
-  if (ptr1) {
-    *ptr1 = '\0';
-    ptr2 = find_next_comma(ptr1+1);
-    if (ptr2) *ptr2 = '\0';
-  } else ptr2 = NULL;
-
-  int n = strlen(contents) + 1;
-  arg1 = new char[n];
-  strcpy(arg1,contents);
-  int narg = 1;
-  if (ptr1) {
-    n = strlen(ptr1+1) + 1;
-    arg2 = new char[n];
-    strcpy(arg2,ptr1+1);
-    narg = 2;
-  } else arg2 = NULL;
-  if (ptr2) {
-    n = strlen(ptr2+1) + 1;
-    arg3 = new char[n];
-    strcpy(arg3,ptr2+1);
-    narg = 3;
-  } else arg3 = NULL;
+  // parse contents for comma-separated args
+  // narg = number of args, args = strings between commas
+  
+  char *args[MAXFUNCARG];
+  int narg = parse_args(contents,args);
 
   // group to operate on
 
-  int igroup = group->find(arg1);
+  int igroup = group->find(args[0]);
   if (igroup == -1)
     error->all(FLERR,"Group ID in variable formula does not exist");
 
@@ -2991,17 +3075,17 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
 
   if (strcmp(word,"count") == 0) {
     if (narg == 1) value = group->count(igroup);
-    else if (narg == 2) value = group->count(igroup,region_function(arg2));
+    else if (narg == 2) value = group->count(igroup,region_function(args[1]));
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"mass") == 0) {
     if (narg == 1) value = group->mass(igroup);
-    else if (narg == 2) value = group->mass(igroup,region_function(arg2));
+    else if (narg == 2) value = group->mass(igroup,region_function(args[1]));
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"charge") == 0) {
     if (narg == 1) value = group->charge(igroup);
-    else if (narg == 2) value = group->charge(igroup,region_function(arg2));
+    else if (narg == 2) value = group->charge(igroup,region_function(args[1]));
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"xcm") == 0) {
@@ -3011,13 +3095,13 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       double masstotal = group->mass(igroup);
       group->xcm(igroup,masstotal,xcm);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = xcm[0];
-    else if (strcmp(arg2,"y") == 0) value = xcm[1];
-    else if (strcmp(arg2,"z") == 0) value = xcm[2];
+    if (strcmp(args[1],"x") == 0) value = xcm[0];
+    else if (strcmp(args[1],"y") == 0) value = xcm[1];
+    else if (strcmp(args[1],"z") == 0) value = xcm[2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"vcm") == 0) {
@@ -3027,36 +3111,36 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       double masstotal = group->mass(igroup);
       group->vcm(igroup,masstotal,vcm);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->vcm(igroup,masstotal,vcm,iregion);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = vcm[0];
-    else if (strcmp(arg2,"y") == 0) value = vcm[1];
-    else if (strcmp(arg2,"z") == 0) value = vcm[2];
+    if (strcmp(args[1],"x") == 0) value = vcm[0];
+    else if (strcmp(args[1],"y") == 0) value = vcm[1];
+    else if (strcmp(args[1],"z") == 0) value = vcm[2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"fcm") == 0) {
     double fcm[3];
     if (narg == 2) group->fcm(igroup,fcm);
-    else if (narg == 3) group->fcm(igroup,fcm,region_function(arg3));
+    else if (narg == 3) group->fcm(igroup,fcm,region_function(args[2]));
     else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = fcm[0];
-    else if (strcmp(arg2,"y") == 0) value = fcm[1];
-    else if (strcmp(arg2,"z") == 0) value = fcm[2];
+    if (strcmp(args[1],"x") == 0) value = fcm[0];
+    else if (strcmp(args[1],"y") == 0) value = fcm[1];
+    else if (strcmp(args[1],"z") == 0) value = fcm[2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"bound") == 0) {
     double minmax[6];
     if (narg == 2) group->bounds(igroup,minmax);
-    else if (narg == 3) group->bounds(igroup,minmax,region_function(arg3));
+    else if (narg == 3) group->bounds(igroup,minmax,region_function(args[2]));
     else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"xmin") == 0) value = minmax[0];
-    else if (strcmp(arg2,"xmax") == 0) value = minmax[1];
-    else if (strcmp(arg2,"ymin") == 0) value = minmax[2];
-    else if (strcmp(arg2,"ymax") == 0) value = minmax[3];
-    else if (strcmp(arg2,"zmin") == 0) value = minmax[4];
-    else if (strcmp(arg2,"zmax") == 0) value = minmax[5];
+    if (strcmp(args[1],"xmin") == 0) value = minmax[0];
+    else if (strcmp(args[1],"xmax") == 0) value = minmax[1];
+    else if (strcmp(args[1],"ymin") == 0) value = minmax[2];
+    else if (strcmp(args[1],"ymax") == 0) value = minmax[3];
+    else if (strcmp(args[1],"zmin") == 0) value = minmax[4];
+    else if (strcmp(args[1],"zmax") == 0) value = minmax[5];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"gyration") == 0) {
@@ -3067,7 +3151,7 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       group->xcm(igroup,masstotal,xcm);
       value = group->gyration(igroup,masstotal,xcm);
     } else if (narg == 2) {
-      int iregion = region_function(arg2);
+      int iregion = region_function(args[1]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
       value = group->gyration(igroup,masstotal,xcm,iregion);
@@ -3075,7 +3159,7 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
 
   } else if (strcmp(word,"ke") == 0) {
     if (narg == 1) value = group->ke(igroup);
-    else if (narg == 2) value = group->ke(igroup,region_function(arg2));
+    else if (narg == 2) value = group->ke(igroup,region_function(args[1]));
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"angmom") == 0) {
@@ -3086,14 +3170,14 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       group->xcm(igroup,masstotal,xcm);
       group->angmom(igroup,xcm,lmom);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
       group->angmom(igroup,xcm,lmom,iregion);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = lmom[0];
-    else if (strcmp(arg2,"y") == 0) value = lmom[1];
-    else if (strcmp(arg2,"z") == 0) value = lmom[2];
+    if (strcmp(args[1],"x") == 0) value = lmom[0];
+    else if (strcmp(args[1],"y") == 0) value = lmom[1];
+    else if (strcmp(args[1],"z") == 0) value = lmom[2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"torque") == 0) {
@@ -3104,14 +3188,14 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       group->xcm(igroup,masstotal,xcm);
       group->torque(igroup,xcm,tq);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
       group->torque(igroup,xcm,tq,iregion);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = tq[0];
-    else if (strcmp(arg2,"y") == 0) value = tq[1];
-    else if (strcmp(arg2,"z") == 0) value = tq[2];
+    if (strcmp(args[1],"x") == 0) value = tq[0];
+    else if (strcmp(args[1],"y") == 0) value = tq[1];
+    else if (strcmp(args[1],"z") == 0) value = tq[2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"inertia") == 0) {
@@ -3122,17 +3206,17 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       group->xcm(igroup,masstotal,xcm);
       group->inertia(igroup,xcm,inertia);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
       group->inertia(igroup,xcm,inertia,iregion);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"xx") == 0) value = inertia[0][0];
-    else if (strcmp(arg2,"yy") == 0) value = inertia[1][1];
-    else if (strcmp(arg2,"zz") == 0) value = inertia[2][2];
-    else if (strcmp(arg2,"xy") == 0) value = inertia[0][1];
-    else if (strcmp(arg2,"yz") == 0) value = inertia[1][2];
-    else if (strcmp(arg2,"xz") == 0) value = inertia[0][2];
+    if (strcmp(args[1],"xx") == 0) value = inertia[0][0];
+    else if (strcmp(args[1],"yy") == 0) value = inertia[1][1];
+    else if (strcmp(args[1],"zz") == 0) value = inertia[2][2];
+    else if (strcmp(args[1],"xy") == 0) value = inertia[0][1];
+    else if (strcmp(args[1],"yz") == 0) value = inertia[1][2];
+    else if (strcmp(args[1],"xz") == 0) value = inertia[0][2];
     else error->all(FLERR,"Invalid group function in variable formula");
 
   } else if (strcmp(word,"omega") == 0) {
@@ -3145,22 +3229,22 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       group->inertia(igroup,xcm,inertia);
       group->omega(angmom,inertia,omega);
     } else if (narg == 3) {
-      int iregion = region_function(arg3);
+      int iregion = region_function(args[2]);
       double masstotal = group->mass(igroup,iregion);
       group->xcm(igroup,masstotal,xcm,iregion);
       group->angmom(igroup,xcm,angmom,iregion);
       group->inertia(igroup,xcm,inertia,iregion);
       group->omega(angmom,inertia,omega);
     } else error->all(FLERR,"Invalid group function in variable formula");
-    if (strcmp(arg2,"x") == 0) value = omega[0];
-    else if (strcmp(arg2,"y") == 0) value = omega[1];
-    else if (strcmp(arg2,"z") == 0) value = omega[2];
+    if (strcmp(args[1],"x") == 0) value = omega[0];
+    else if (strcmp(args[1],"y") == 0) value = omega[1];
+    else if (strcmp(args[1],"z") == 0) value = omega[2];
     else error->all(FLERR,"Invalid group function in variable formula");
   }
 
-  delete [] arg1;
-  delete [] arg2;
-  delete [] arg3;
+  // delete stored args
+
+  for (int i = 0; i < narg; i++) delete [] args[i];
 
   // save value in tree or on argstack
 
@@ -3168,7 +3252,8 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
     Tree *newtree = new Tree();
     newtree->type = VALUE;
     newtree->value = value;
-    newtree->left = newtree->middle = newtree->right = NULL;
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     treestack[ntreestack++] = newtree;
   } else argstack[nargstack++] = value;
 
@@ -3215,35 +3300,11 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       strcmp(word,"grmask") && strcmp(word,"next"))
     return 0;
 
-  // parse contents for arg1,arg2,arg3 separated by commas
-  // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
-
-  char *arg1,*arg2,*arg3;
-  char *ptr1,*ptr2;
-
-  ptr1 = find_next_comma(contents);
-  if (ptr1) {
-    *ptr1 = '\0';
-    ptr2 = find_next_comma(ptr1+1);
-    if (ptr2) *ptr2 = '\0';
-  } else ptr2 = NULL;
-
-  int n = strlen(contents) + 1;
-  arg1 = new char[n];
-  strcpy(arg1,contents);
-  int narg = 1;
-  if (ptr1) {
-    n = strlen(ptr1+1) + 1;
-    arg2 = new char[n];
-    strcpy(arg2,ptr1+1);
-    narg = 2;
-  } else arg2 = NULL;
-  if (ptr2) {
-    n = strlen(ptr2+1) + 1;
-    arg3 = new char[n];
-    strcpy(arg3,ptr2+1);
-    narg = 3;
-  } else arg3 = NULL;
+  // parse contents for comma-separated args
+  // narg = number of args, args = strings between commas
+  
+  char *args[MAXFUNCARG];
+  int narg = parse_args(contents,args);
 
   // special functions that operate on global vectors
 
@@ -3265,16 +3326,17 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     Compute *compute = NULL;
     Fix *fix = NULL;
     int index,nvec,nstride;
+    char *ptr1,*ptr2;
 
-    if (strstr(arg1,"c_") == arg1) {
-      ptr1 = strchr(arg1,'[');
+    if (strstr(args[0],"c_") == args[0]) {
+      ptr1 = strchr(args[0],'[');
       if (ptr1) {
         ptr2 = ptr1;
         index = int_between_brackets(ptr2,0);
         *ptr1 = '\0';
       } else index = 0;
 
-      int icompute = modify->find_compute(&arg1[2]);
+      int icompute = modify->find_compute(&args[0][2]);
       if (icompute < 0) 
         error->all(FLERR,"Invalid compute ID in variable formula");
       compute = modify->compute[icompute];
@@ -3305,15 +3367,15 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         nstride = compute->size_array_cols;
       } else error->all(FLERR,"Mismatched compute in variable formula");
 
-    } else if (strstr(arg1,"f_") == arg1) {
-      ptr1 = strchr(arg1,'[');
+    } else if (strstr(args[0],"f_") == args[0]) {
+      ptr1 = strchr(args[0],'[');
       if (ptr1) {
         ptr2 = ptr1;
         index = int_between_brackets(ptr2,0);
         *ptr1 = '\0';
       } else index = 0;
 
-      int ifix = modify->find_fix(&arg1[2]);
+      int ifix = modify->find_fix(&args[0][2]);
       if (ifix < 0) error->all(FLERR,"Invalid fix ID in variable formula");
       fix = modify->fix[ifix];
       if (index == 0 && fix->vector_flag) {
@@ -3407,7 +3469,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       Tree *newtree = new Tree();
       newtree->type = VALUE;
       newtree->value = value;
-      newtree->left = newtree->middle = newtree->right = NULL;
+      newtree->first = newtree->second = NULL;
+      newtree->nextra = 0;
       treestack[ntreestack++] = newtree;
     } else argstack[nargstack++] = value;
 
@@ -3419,14 +3482,15 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     if (narg != 1) 
       error->all(FLERR,"Invalid special function in variable formula");
 
-    int igroup = group->find(arg1);
+    int igroup = group->find(args[0]);
     if (igroup == -1)
       error->all(FLERR,"Group ID in variable formula does not exist");
 
     Tree *newtree = new Tree();
     newtree->type = GMASK;
     newtree->ivalue1 = group->bitmask[igroup];
-    newtree->left = newtree->middle = newtree->right = NULL;
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     treestack[ntreestack++] = newtree;
 
   } else if (strcmp(word,"rmask") == 0) {
@@ -3435,13 +3499,14 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     if (narg != 1) 
       error->all(FLERR,"Invalid special function in variable formula");
 
-    int iregion = region_function(arg1);
+    int iregion = region_function(args[0]);
     domain->regions[iregion]->prematch();
 
     Tree *newtree = new Tree();
     newtree->type = RMASK;
     newtree->ivalue1 = iregion;
-    newtree->left = newtree->middle = newtree->right = NULL;
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     treestack[ntreestack++] = newtree;
 
   } else if (strcmp(word,"grmask") == 0) {
@@ -3450,17 +3515,18 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     if (narg != 2) 
       error->all(FLERR,"Invalid special function in variable formula");
 
-    int igroup = group->find(arg1);
+    int igroup = group->find(args[0]);
     if (igroup == -1)
       error->all(FLERR,"Group ID in variable formula does not exist");
-    int iregion = region_function(arg2);
+    int iregion = region_function(args[1]);
     domain->regions[iregion]->prematch();
 
     Tree *newtree = new Tree();
     newtree->type = GRMASK;
     newtree->ivalue1 = group->bitmask[igroup];
     newtree->ivalue2 = iregion;
-    newtree->left = newtree->middle = newtree->right = NULL;
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     treestack[ntreestack++] = newtree;
 
   // special function for file-style or atomfile-style variables
@@ -3469,7 +3535,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     if (narg != 1) 
       error->all(FLERR,"Invalid special function in variable formula");
 
-    int ivar = find(arg1);
+    int ivar = find(args[0]);
     if (ivar == -1)
       error->all(FLERR,"Variable ID in variable formula does not exist");
 
@@ -3485,7 +3551,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         Tree *newtree = new Tree();
         newtree->type = VALUE;
         newtree->value = value;
-        newtree->left = newtree->middle = newtree->right = NULL;
+        newtree->first = newtree->second = NULL;
+        newtree->nextra = 0;
         treestack[ntreestack++] = newtree;
       } else argstack[nargstack++] = value;
 
@@ -3509,15 +3576,16 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       newtree->array = result;
       newtree->nstride = 1;
       newtree->selfalloc = 1;
-      newtree->left = newtree->middle = newtree->right = NULL;
+      newtree->first = newtree->second = NULL;
+      newtree->nextra = 0;
       treestack[ntreestack++] = newtree;
 
     } else error->all(FLERR,"Invalid variable style in special function next");
   }
 
-  delete [] arg1;
-  delete [] arg2;
-  delete [] arg3;
+  // delete stored args
+
+  for (int i = 0; i < narg; i++) delete [] args[i];
 
   return 1;
 }
@@ -3585,7 +3653,8 @@ void Variable::peratom2global(int flag, char *word,
     Tree *newtree = new Tree();
     newtree->type = VALUE;
     newtree->value = value;
-    newtree->left = newtree->middle = newtree->right = NULL;
+    newtree->first = newtree->second = NULL;
+    newtree->nextra = 0;
     treestack[ntreestack++] = newtree;
   } else argstack[nargstack++] = value;
 }
@@ -3634,7 +3703,8 @@ void Variable::atom_vector(char *word, Tree **tree,
   newtree->type = ATOMARRAY;
   newtree->nstride = 3;
   newtree->selfalloc = 0;
-  newtree->left = newtree->middle = newtree->right = NULL;
+  newtree->first = newtree->second = NULL;
+  newtree->nextra = 0;
   treestack[ntreestack++] = newtree;
 
   if (strcmp(word,"id") == 0) {
@@ -3713,6 +3783,36 @@ double Variable::constant(char *word)
   return 0.0;
 }
 
+/* ----------------------------------------------------------------------
+   parse string for comma-separated args
+   store copy of each arg in args array
+   max allowed # of args = MAXFUNCARG
+------------------------------------------------------------------------- */
+
+int Variable::parse_args(char *str, char **args)
+{
+  int n;
+  char *ptrnext;
+
+  int narg = 0;
+  char *ptr = str;
+
+  while (ptr && narg < MAXFUNCARG) {
+    ptrnext = find_next_comma(ptr);
+    if (ptrnext) *ptrnext = '\0';
+    n = strlen(ptr) + 1;
+    args[narg] = new char[n];
+    strcpy(args[narg],ptr);
+    narg++;
+    ptr = ptrnext;
+    if (ptr) ptr++;
+  }
+
+  if (ptr) error->all(FLERR,"Too many args in variable function");
+  return narg;
+}
+
+
 /* ----------------------------------------------------------------------
    find next comma in str
    skip commas inside one or more nested parenthesis
@@ -3737,9 +3837,10 @@ char *Variable::find_next_comma(char *str)
 void Variable::print_tree(Tree *tree, int level)
 {
   printf("TREE %d: %d %g\n",level,tree->type,tree->value);
-  if (tree->left) print_tree(tree->left,level+1);
-  if (tree->middle) print_tree(tree->middle,level+1);
-  if (tree->right) print_tree(tree->right,level+1);
+  if (tree->first) print_tree(tree->first,level+1);
+  if (tree->second) print_tree(tree->second,level+1);
+  if (tree->nextra)
+    for (int i = 0; i < tree->nextra; i++) print_tree(tree->extra[i],level+1);
   return;
 }
 
diff --git a/src/variable.h b/src/variable.h
index f04a04f9d4..84b6d04992 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -72,7 +72,9 @@ class Variable : protected Pointers {
     int nstride;           // stride between atoms if array is a 2d array
     int selfalloc;         // 1 if array is allocated here, else 0
     int ivalue1,ivalue2;   // extra values for needed for gmask,rmask,grmask
-    Tree *left,*middle,*right;    // ptrs further down tree
+    int nextra;            // # of additional args beyond first 2
+    Tree *first,*second;   // ptrs further down tree for first 2 args
+    Tree **extra;          // ptrs further down tree for nextra args
   };
 
   void remove(int);
@@ -94,6 +96,7 @@ class Variable : protected Pointers {
   void atom_vector(char *, Tree **, Tree **, int &);
   int is_constant(char *);
   double constant(char *);
+  int parse_args(char *, char **);
   char *find_next_comma(char *);
   void print_tree(Tree *, int);
 };

From 1591ed8379a6bdfc73a2d3f575a5d5f01eaf01ae Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 15:07:55 +0000
Subject: [PATCH 02/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12728
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/variable.html | 25 +++++++++++++++++++++----
 doc/variable.txt  | 25 +++++++++++++++++++++----
 2 files changed, 42 insertions(+), 8 deletions(-)

diff --git a/doc/variable.html b/doc/variable.html
index 0ff65ea228..666b971dd9 100644
--- a/doc/variable.html
+++ b/doc/variable.html
@@ -54,7 +54,7 @@
     math functions = sqrt(x), exp(x), ln(x), log(x), abs(x),
                      sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x),
                      random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x)
-                     ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
+                     ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
     group functions = count(group), mass(group), charge(group),
 		      xcm(group,dim), vcm(group,dim), fcm(group,dim),
 		      bound(group,dir), gyration(group), ke(group),
@@ -378,7 +378,7 @@ references to other variables.
 <TR><TD >Thermo keywords</TD><TD > vol, pe, ebond, etc</TD></TR>
 <TR><TD >Math operators</TD><TD > (), -x, x+y, x-y, x*y, x/y, x^y, x%y, </TD></TR>
 <TR><TD >Math operators</TD><TD > (), -x, x+y, x-y, x*y, x/y, x^y, x%y, x == y, x != y, x &lt y, x &lt= y, x &gt y, x &gt= y, x && y, x || y, !x</TD></TR>
-<TR><TD >Math functions</TD><TD > sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)</TD></TR>
+<TR><TD >Math functions</TD><TD > sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)</TD></TR>
 <TR><TD >Group functions</TD><TD > count(ID), mass(ID), charge(ID), xcm(ID,dim),                  vcm(ID,dim), fcm(ID,dim), bound(ID,dir), 		 gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), 		 inertia(ID,dimdim), omega(ID,dim)</TD></TR>
 <TR><TD >Region functions</TD><TD > count(ID,IDR), mass(ID,IDR), charge(ID,IDR), 		  xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), 		  bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), 		  angmom(ID,dim,IDR), torque(ID,dim,IDR),                   inertia(ID,dimdim,IDR), omega(ID,dim,IDR)</TD></TR>
 <TR><TD >Special functions</TD><TD > sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x)</TD></TR>
@@ -542,16 +542,33 @@ output timesteps:
 </PRE>
 <P>The stride(x,y,z) function uses the current timestep to generate a new
 timestep.  X,y >= 0 and z > 0 and x <= y are required.  The generated
-timesteps increase in increments of z, from x to y, I.e. it generates
+timesteps increase in increments of z, from x to y, i.e. it generates
 the sequece x,x+z,x+2z,...,y.  If y-x is not a multiple of z, then
 similar to the way a for loop operates, the last value will be one
 that does not exceed y.  For any current timestep, the next timestep
-in the sequence is returned.  Thus if stagger(1000,2000,100) is used
+in the sequence is returned.  Thus if stride(1000,2000,100) is used
 in a variable by the <A HREF = "dump_modify.html">dump_modify every</A> command, it
 will generate the sequence of output timesteps:
 </P>
 <PRE>1000,1100,1200, ... ,1900,2000 
 </PRE>
+<P>The stride2(x,y,z,a,b,c) function is similar to the stride() function
+except it generates two sets of strided timesteps, one at a coarser
+level and one at a finer level.  Thus it is useful for debugging,
+e.g. to produce output every timestep at the point in simulation when
+a problem occurs.  X,y >= 0 and z > 0 and x <= y are required, as are
+a,b >= 0 and c > 0 and a < b.  Also, a >= x and b <= y are required so
+that the second stride is inside the first.  The generated timesteps
+increase in increments of z, starting at x, until a is reached.  At
+that point the timestep increases in increments of c, from a to b,
+then after b, increments by z are resumed until y is reached.  For any
+current timestep, the next timestep in the sequence is returned.  Thus
+if stride(1000,2000,100,1350,1360,1) is used in a variable by the
+<A HREF = "dump_modify.html">dump_modify every</A> command, it will generate the
+sequence of output timesteps:
+</P>
+<PRE>1000,1100,1200,1300,1350,1351,1352, ... 1359,1360,1400,1500, ... ,2000 
+</PRE>
 <P>The vdisplace(x,y) function takes 2 arguments: x = value0 and y =
 velocity, and uses the elapsed time to change the value by a linear
 displacement due to the applied velocity over the course of a run,
diff --git a/doc/variable.txt b/doc/variable.txt
index c2805fed8a..868c0d70b3 100644
--- a/doc/variable.txt
+++ b/doc/variable.txt
@@ -49,7 +49,7 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
     math functions = sqrt(x), exp(x), ln(x), log(x), abs(x),
                      sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x),
                      random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x)
-                     ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
+                     ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
     group functions = count(group), mass(group), charge(group),
 		      xcm(group,dim), vcm(group,dim), fcm(group,dim),
 		      bound(group,dir), gyration(group), ke(group),
@@ -371,7 +371,7 @@ Constant: PI
 Thermo keywords: vol, pe, ebond, etc
 Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x%y, 
 Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x%y, x == y, x != y, x &lt y, x &lt= y, x &gt y, x &gt= y, x && y, x || y, !x
-Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
+Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
 Group functions: count(ID), mass(ID), charge(ID), xcm(ID,dim), \
                  vcm(ID,dim), fcm(ID,dim), bound(ID,dir), \
 		 gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), \
@@ -542,16 +542,33 @@ output timesteps:
 
 The stride(x,y,z) function uses the current timestep to generate a new
 timestep.  X,y >= 0 and z > 0 and x <= y are required.  The generated
-timesteps increase in increments of z, from x to y, I.e. it generates
+timesteps increase in increments of z, from x to y, i.e. it generates
 the sequece x,x+z,x+2z,...,y.  If y-x is not a multiple of z, then
 similar to the way a for loop operates, the last value will be one
 that does not exceed y.  For any current timestep, the next timestep
-in the sequence is returned.  Thus if stagger(1000,2000,100) is used
+in the sequence is returned.  Thus if stride(1000,2000,100) is used
 in a variable by the "dump_modify every"_dump_modify.html command, it
 will generate the sequence of output timesteps:
 
 1000,1100,1200, ... ,1900,2000 :pre
 
+The stride2(x,y,z,a,b,c) function is similar to the stride() function
+except it generates two sets of strided timesteps, one at a coarser
+level and one at a finer level.  Thus it is useful for debugging,
+e.g. to produce output every timestep at the point in simulation when
+a problem occurs.  X,y >= 0 and z > 0 and x <= y are required, as are
+a,b >= 0 and c > 0 and a < b.  Also, a >= x and b <= y are required so
+that the second stride is inside the first.  The generated timesteps
+increase in increments of z, starting at x, until a is reached.  At
+that point the timestep increases in increments of c, from a to b,
+then after b, increments by z are resumed until y is reached.  For any
+current timestep, the next timestep in the sequence is returned.  Thus
+if stride(1000,2000,100,1350,1360,1) is used in a variable by the
+"dump_modify every"_dump_modify.html command, it will generate the
+sequence of output timesteps:
+
+1000,1100,1200,1300,1350,1351,1352, ... 1359,1360,1400,1500, ... ,2000 :pre
+
 The vdisplace(x,y) function takes 2 arguments: x = value0 and y =
 velocity, and uses the elapsed time to change the value by a linear
 displacement due to the applied velocity over the course of a run,

From 74adae232eb463a934296a942eb54bbc509c198c Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 15:25:25 +0000
Subject: [PATCH 03/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12729
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/KOKKOS/atom_kokkos.cpp               |  2 +-
 src/KOKKOS/atom_vec_angle_kokkos.cpp     | 32 +++++++-------
 src/KOKKOS/atom_vec_bond_kokkos.cpp      | 24 +++++-----
 src/KOKKOS/atom_vec_charge_kokkos.cpp    |  8 ++--
 src/KOKKOS/atom_vec_full_kokkos.cpp      | 56 ++++++++++++------------
 src/KOKKOS/atom_vec_molecular_kokkos.cpp | 48 ++++++++++----------
 6 files changed, 85 insertions(+), 85 deletions(-)

diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp
index ec6dc607bd..910c0f9996 100644
--- a/src/KOKKOS/atom_kokkos.cpp
+++ b/src/KOKKOS/atom_kokkos.cpp
@@ -205,7 +205,7 @@ void AtomKokkos::sort()
 
 void AtomKokkos::grow(unsigned int mask){
 
-  if (mask && SPECIAL_MASK){
+  if (mask & SPECIAL_MASK){
     memory->destroy_kokkos(k_special, special);
     sync(Device, mask);
     modified(Device, mask);
diff --git a/src/KOKKOS/atom_vec_angle_kokkos.cpp b/src/KOKKOS/atom_vec_angle_kokkos.cpp
index 0687e7dfb9..e1f73064d7 100644
--- a/src/KOKKOS/atom_vec_angle_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_angle_kokkos.cpp
@@ -1777,17 +1777,17 @@ void AtomVecAngleKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPDeviceType>();
       atomKK->k_special.sync<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPDeviceType>();
       atomKK->k_bond_type.sync<LMPDeviceType>();
       atomKK->k_bond_atom.sync<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPDeviceType>();
       atomKK->k_angle_type.sync<LMPDeviceType>();
       atomKK->k_angle_atom1.sync<LMPDeviceType>();
@@ -1802,17 +1802,17 @@ void AtomVecAngleKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPHostType>();
       atomKK->k_special.sync<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPHostType>();
       atomKK->k_bond_type.sync<LMPHostType>();
       atomKK->k_bond_atom.sync<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPHostType>();
       atomKK->k_angle_type.sync<LMPHostType>();
       atomKK->k_angle_atom1.sync<LMPHostType>();
@@ -1834,17 +1834,17 @@ void AtomVecAngleKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPDeviceType>();
       atomKK->k_special.modify<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPDeviceType>();
       atomKK->k_bond_type.modify<LMPDeviceType>();
       atomKK->k_bond_atom.modify<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPDeviceType>();
       atomKK->k_angle_type.modify<LMPDeviceType>();
       atomKK->k_angle_atom1.modify<LMPDeviceType>();
@@ -1859,17 +1859,17 @@ void AtomVecAngleKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPHostType>();
       atomKK->k_special.modify<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPHostType>();
       atomKK->k_bond_type.modify<LMPHostType>();
       atomKK->k_bond_atom.modify<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPHostType>();
       atomKK->k_angle_type.modify<LMPHostType>();
       atomKK->k_angle_atom1.modify<LMPHostType>();
diff --git a/src/KOKKOS/atom_vec_bond_kokkos.cpp b/src/KOKKOS/atom_vec_bond_kokkos.cpp
index edac0bcfdc..a5ed6163a7 100644
--- a/src/KOKKOS/atom_vec_bond_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_bond_kokkos.cpp
@@ -1636,12 +1636,12 @@ void AtomVecBondKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPDeviceType>();
       atomKK->k_special.sync<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPDeviceType>();
       atomKK->k_bond_type.sync<LMPDeviceType>();
       atomKK->k_bond_atom.sync<LMPDeviceType>();
@@ -1654,12 +1654,12 @@ void AtomVecBondKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPHostType>();
       atomKK->k_special.sync<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPHostType>();
       atomKK->k_bond_type.sync<LMPHostType>();
       atomKK->k_bond_atom.sync<LMPHostType>();
@@ -1679,12 +1679,12 @@ void AtomVecBondKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPDeviceType>();
       atomKK->k_special.modify<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPDeviceType>();
       atomKK->k_bond_type.modify<LMPDeviceType>();
       atomKK->k_bond_atom.modify<LMPDeviceType>();
@@ -1697,12 +1697,12 @@ void AtomVecBondKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPHostType>();
       atomKK->k_special.modify<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPHostType>();
       atomKK->k_bond_type.modify<LMPHostType>();
       atomKK->k_bond_atom.modify<LMPHostType>();
diff --git a/src/KOKKOS/atom_vec_charge_kokkos.cpp b/src/KOKKOS/atom_vec_charge_kokkos.cpp
index 3b2f48bb66..0bcc594f75 100644
--- a/src/KOKKOS/atom_vec_charge_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_charge_kokkos.cpp
@@ -1475,7 +1475,7 @@ void AtomVecChargeKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPDeviceType>();
-    if (mask && Q_MASK) atomKK->k_q.sync<LMPDeviceType>();
+    if (mask & Q_MASK) atomKK->k_q.sync<LMPDeviceType>();
   } else {
     if (mask & X_MASK) atomKK->k_x.sync<LMPHostType>();
     if (mask & V_MASK) atomKK->k_v.sync<LMPHostType>();
@@ -1484,7 +1484,7 @@ void AtomVecChargeKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPHostType>();
-    if (mask && Q_MASK) atomKK->k_q.sync<LMPHostType>();
+    if (mask & Q_MASK) atomKK->k_q.sync<LMPHostType>();
   }
 }
 
@@ -1500,7 +1500,7 @@ void AtomVecChargeKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPDeviceType>();
-    if (mask && Q_MASK) atomKK->k_q.modify<LMPDeviceType>();
+    if (mask & Q_MASK) atomKK->k_q.modify<LMPDeviceType>();
   } else {
     if (mask & X_MASK) atomKK->k_x.modify<LMPHostType>();
     if (mask & V_MASK) atomKK->k_v.modify<LMPHostType>();
@@ -1509,6 +1509,6 @@ void AtomVecChargeKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPHostType>();
-    if (mask && Q_MASK) atomKK->k_q.modify<LMPHostType>();
+    if (mask & Q_MASK) atomKK->k_q.modify<LMPHostType>();
   }
 }
diff --git a/src/KOKKOS/atom_vec_full_kokkos.cpp b/src/KOKKOS/atom_vec_full_kokkos.cpp
index 2cb87f9bcb..6623fa2f25 100644
--- a/src/KOKKOS/atom_vec_full_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_full_kokkos.cpp
@@ -2144,25 +2144,25 @@ void AtomVecFullKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPDeviceType>();
-    if (mask && Q_MASK) atomKK->k_q.sync<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & Q_MASK) atomKK->k_q.sync<LMPDeviceType>();
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPDeviceType>();
       atomKK->k_special.sync<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPDeviceType>();
       atomKK->k_bond_type.sync<LMPDeviceType>();
       atomKK->k_bond_atom.sync<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPDeviceType>();
       atomKK->k_angle_type.sync<LMPDeviceType>();
       atomKK->k_angle_atom1.sync<LMPDeviceType>();
       atomKK->k_angle_atom2.sync<LMPDeviceType>();
       atomKK->k_angle_atom3.sync<LMPDeviceType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.sync<LMPDeviceType>();
       atomKK->k_dihedral_type.sync<LMPDeviceType>();
       atomKK->k_dihedral_atom1.sync<LMPDeviceType>();
@@ -2170,7 +2170,7 @@ void AtomVecFullKokkos::sync(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.sync<LMPDeviceType>();
       atomKK->k_dihedral_atom4.sync<LMPDeviceType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.sync<LMPDeviceType>();
       atomKK->k_improper_type.sync<LMPDeviceType>();
       atomKK->k_improper_atom1.sync<LMPDeviceType>();
@@ -2186,25 +2186,25 @@ void AtomVecFullKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPHostType>();
-    if (mask && Q_MASK) atomKK->k_q.sync<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & Q_MASK) atomKK->k_q.sync<LMPHostType>();
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPHostType>();
       atomKK->k_special.sync<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPHostType>();
       atomKK->k_bond_type.sync<LMPHostType>();
       atomKK->k_bond_atom.sync<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPHostType>();
       atomKK->k_angle_type.sync<LMPHostType>();
       atomKK->k_angle_atom1.sync<LMPHostType>();
       atomKK->k_angle_atom2.sync<LMPHostType>();
       atomKK->k_angle_atom3.sync<LMPHostType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.sync<LMPHostType>();
       atomKK->k_dihedral_type.sync<LMPHostType>();
       atomKK->k_dihedral_atom1.sync<LMPHostType>();
@@ -2212,7 +2212,7 @@ void AtomVecFullKokkos::sync(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.sync<LMPHostType>();
       atomKK->k_dihedral_atom4.sync<LMPHostType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.sync<LMPHostType>();
       atomKK->k_improper_type.sync<LMPHostType>();
       atomKK->k_improper_atom1.sync<LMPHostType>();
@@ -2235,25 +2235,25 @@ void AtomVecFullKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPDeviceType>();
-    if (mask && Q_MASK) atomKK->k_q.modify<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & Q_MASK) atomKK->k_q.modify<LMPDeviceType>();
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPDeviceType>();
       atomKK->k_special.modify<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPDeviceType>();
       atomKK->k_bond_type.modify<LMPDeviceType>();
       atomKK->k_bond_atom.modify<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPDeviceType>();
       atomKK->k_angle_type.modify<LMPDeviceType>();
       atomKK->k_angle_atom1.modify<LMPDeviceType>();
       atomKK->k_angle_atom2.modify<LMPDeviceType>();
       atomKK->k_angle_atom3.modify<LMPDeviceType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.modify<LMPDeviceType>();
       atomKK->k_dihedral_type.modify<LMPDeviceType>();
       atomKK->k_dihedral_atom1.modify<LMPDeviceType>();
@@ -2261,7 +2261,7 @@ void AtomVecFullKokkos::modified(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.modify<LMPDeviceType>();
       atomKK->k_dihedral_atom4.modify<LMPDeviceType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.modify<LMPDeviceType>();
       atomKK->k_improper_type.modify<LMPDeviceType>();
       atomKK->k_improper_atom1.modify<LMPDeviceType>();
@@ -2277,25 +2277,25 @@ void AtomVecFullKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPHostType>();
-    if (mask && Q_MASK) atomKK->k_q.modify<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & Q_MASK) atomKK->k_q.modify<LMPHostType>();
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPHostType>();
       atomKK->k_special.modify<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPHostType>();
       atomKK->k_bond_type.modify<LMPHostType>();
       atomKK->k_bond_atom.modify<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPHostType>();
       atomKK->k_angle_type.modify<LMPHostType>();
       atomKK->k_angle_atom1.modify<LMPHostType>();
       atomKK->k_angle_atom2.modify<LMPHostType>();
       atomKK->k_angle_atom3.modify<LMPHostType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.modify<LMPHostType>();
       atomKK->k_dihedral_type.modify<LMPHostType>();
       atomKK->k_dihedral_atom1.modify<LMPHostType>();
@@ -2303,7 +2303,7 @@ void AtomVecFullKokkos::modified(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.modify<LMPHostType>();
       atomKK->k_dihedral_atom4.modify<LMPHostType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.modify<LMPHostType>();
       atomKK->k_improper_type.modify<LMPHostType>();
       atomKK->k_improper_atom1.modify<LMPHostType>();
diff --git a/src/KOKKOS/atom_vec_molecular_kokkos.cpp b/src/KOKKOS/atom_vec_molecular_kokkos.cpp
index 4a45fb86d6..7c48b2dc85 100644
--- a/src/KOKKOS/atom_vec_molecular_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_molecular_kokkos.cpp
@@ -2063,24 +2063,24 @@ void AtomVecMolecularKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPDeviceType>();
       atomKK->k_special.sync<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPDeviceType>();
       atomKK->k_bond_type.sync<LMPDeviceType>();
       atomKK->k_bond_atom.sync<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPDeviceType>();
       atomKK->k_angle_type.sync<LMPDeviceType>();
       atomKK->k_angle_atom1.sync<LMPDeviceType>();
       atomKK->k_angle_atom2.sync<LMPDeviceType>();
       atomKK->k_angle_atom3.sync<LMPDeviceType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.sync<LMPDeviceType>();
       atomKK->k_dihedral_type.sync<LMPDeviceType>();
       atomKK->k_dihedral_atom1.sync<LMPDeviceType>();
@@ -2088,7 +2088,7 @@ void AtomVecMolecularKokkos::sync(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.sync<LMPDeviceType>();
       atomKK->k_dihedral_atom4.sync<LMPDeviceType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.sync<LMPDeviceType>();
       atomKK->k_improper_type.sync<LMPDeviceType>();
       atomKK->k_improper_atom1.sync<LMPDeviceType>();
@@ -2104,24 +2104,24 @@ void AtomVecMolecularKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.sync<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.sync<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.sync<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.sync<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.sync<LMPHostType>();
       atomKK->k_special.sync<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.sync<LMPHostType>();
       atomKK->k_bond_type.sync<LMPHostType>();
       atomKK->k_bond_atom.sync<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.sync<LMPHostType>();
       atomKK->k_angle_type.sync<LMPHostType>();
       atomKK->k_angle_atom1.sync<LMPHostType>();
       atomKK->k_angle_atom2.sync<LMPHostType>();
       atomKK->k_angle_atom3.sync<LMPHostType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.sync<LMPHostType>();
       atomKK->k_dihedral_type.sync<LMPHostType>();
       atomKK->k_dihedral_atom1.sync<LMPHostType>();
@@ -2129,7 +2129,7 @@ void AtomVecMolecularKokkos::sync(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.sync<LMPHostType>();
       atomKK->k_dihedral_atom4.sync<LMPHostType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.sync<LMPHostType>();
       atomKK->k_improper_type.sync<LMPHostType>();
       atomKK->k_improper_atom1.sync<LMPHostType>();
@@ -2152,24 +2152,24 @@ void AtomVecMolecularKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPDeviceType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPDeviceType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPDeviceType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPDeviceType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPDeviceType>();
       atomKK->k_special.modify<LMPDeviceType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPDeviceType>();
       atomKK->k_bond_type.modify<LMPDeviceType>();
       atomKK->k_bond_atom.modify<LMPDeviceType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPDeviceType>();
       atomKK->k_angle_type.modify<LMPDeviceType>();
       atomKK->k_angle_atom1.modify<LMPDeviceType>();
       atomKK->k_angle_atom2.modify<LMPDeviceType>();
       atomKK->k_angle_atom3.modify<LMPDeviceType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.modify<LMPDeviceType>();
       atomKK->k_dihedral_type.modify<LMPDeviceType>();
       atomKK->k_dihedral_atom1.modify<LMPDeviceType>();
@@ -2177,7 +2177,7 @@ void AtomVecMolecularKokkos::modified(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.modify<LMPDeviceType>();
       atomKK->k_dihedral_atom4.modify<LMPDeviceType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.modify<LMPDeviceType>();
       atomKK->k_improper_type.modify<LMPDeviceType>();
       atomKK->k_improper_atom1.modify<LMPDeviceType>();
@@ -2193,24 +2193,24 @@ void AtomVecMolecularKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & TYPE_MASK) atomKK->k_type.modify<LMPHostType>();
     if (mask & MASK_MASK) atomKK->k_mask.modify<LMPHostType>();
     if (mask & IMAGE_MASK) atomKK->k_image.modify<LMPHostType>();
-    if (mask && MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
-    if (mask && SPECIAL_MASK) {
+    if (mask & MOLECULE_MASK) atomKK->k_molecule.modify<LMPHostType>();
+    if (mask & SPECIAL_MASK) {
       atomKK->k_nspecial.modify<LMPHostType>();
       atomKK->k_special.modify<LMPHostType>();
     }
-    if (mask && BOND_MASK) {
+    if (mask & BOND_MASK) {
       atomKK->k_num_bond.modify<LMPHostType>();
       atomKK->k_bond_type.modify<LMPHostType>();
       atomKK->k_bond_atom.modify<LMPHostType>();
     }
-    if (mask && ANGLE_MASK) {
+    if (mask & ANGLE_MASK) {
       atomKK->k_num_angle.modify<LMPHostType>();
       atomKK->k_angle_type.modify<LMPHostType>();
       atomKK->k_angle_atom1.modify<LMPHostType>();
       atomKK->k_angle_atom2.modify<LMPHostType>();
       atomKK->k_angle_atom3.modify<LMPHostType>();
     }
-    if (mask && DIHEDRAL_MASK) {
+    if (mask & DIHEDRAL_MASK) {
       atomKK->k_num_dihedral.modify<LMPHostType>();
       atomKK->k_dihedral_type.modify<LMPHostType>();
       atomKK->k_dihedral_atom1.modify<LMPHostType>();
@@ -2218,7 +2218,7 @@ void AtomVecMolecularKokkos::modified(ExecutionSpace space, unsigned int mask)
       atomKK->k_dihedral_atom3.modify<LMPHostType>();
       atomKK->k_dihedral_atom4.modify<LMPHostType>();
     }
-    if (mask && IMPROPER_MASK) {
+    if (mask & IMPROPER_MASK) {
       atomKK->k_num_improper.modify<LMPHostType>();
       atomKK->k_improper_type.modify<LMPHostType>();
       atomKK->k_improper_atom1.modify<LMPHostType>();

From aaa76b5649ae90d7cfe5eae190d5dcc4756a515e Mon Sep 17 00:00:00 2001
From: stamoor <stamoor@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 15:56:24 +0000
Subject: [PATCH 04/17] Fixing bug

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12730 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
 src/kspace.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/kspace.cpp b/src/kspace.cpp
index 2e7f9a8736..5942a5d0c1 100644
--- a/src/kspace.cpp
+++ b/src/kspace.cpp
@@ -357,7 +357,7 @@ void KSpace::x2lamdaT(double *v, double *lamda)
 
 void KSpace::lamda2xT(double *lamda, double *v)
 {
-  double h[5];
+  double h[6];
   h[0] = domain->h[0];
   h[1] = domain->h[1];
   h[2] = domain->h[2];

From 526b81b955c0884028f4dfe4c1e7a92bb7785320 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 17:28:24 +0000
Subject: [PATCH 05/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12731
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index 172ce61278..0bd9cbf0da 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "20 Nov 2014"
+#define LAMMPS_VERSION "21 Nov 2014"

From 0617d87047978aca4593aa80c9068e7b9fe81d02 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 17:28:25 +0000
Subject: [PATCH 06/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12732
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Manual.html | 4 ++--
 doc/Manual.txt  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/Manual.html b/doc/Manual.html
index 147827f673..479053766f 100644
--- a/doc/Manual.html
+++ b/doc/Manual.html
@@ -1,7 +1,7 @@
 <HTML>
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="20 Nov 2014 version">
+<META NAME="docnumber" CONTENT="21 Nov 2014 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@@ -22,7 +22,7 @@
 
 <CENTER><H3>LAMMPS Documentation 
 </H3></CENTER>
-<CENTER><H4>20 Nov 2014 version 
+<CENTER><H4>21 Nov 2014 version 
 </H4></CENTER>
 <H4>Version info: 
 </H4>
diff --git a/doc/Manual.txt b/doc/Manual.txt
index e7726ffb08..c7bb138f20 100644
--- a/doc/Manual.txt
+++ b/doc/Manual.txt
@@ -1,6 +1,6 @@
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="20 Nov 2014 version">
+<META NAME="docnumber" CONTENT="21 Nov 2014 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@@ -18,7 +18,7 @@
 <H1></H1>
 
 LAMMPS Documentation :c,h3
-20 Nov 2014 version :c,h4
+21 Nov 2014 version :c,h4
 
 Version info: :h4
 

From a7c5a1661d2510fceb6abc0efdb3a0b8b5d93df5 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 17:42:32 +0000
Subject: [PATCH 07/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12734
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/GRANULAR/fix_pour.cpp                   | 1 +
 src/RIGID/fix_rigid_small.cpp               | 2 +-
 src/USER-FEP/pair_lj_cut_coul_long_soft.cpp | 4 ++--
 src/USER-FEP/pair_lj_cut_soft.cpp           | 4 ++--
 4 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 9947d9a84a..f2dd25b2f9 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -255,6 +255,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
   }
 
   nper = static_cast<int> (volfrac*volume/volume_one);
+  if (nper == 0) error->all(FLERR,"Fix pour insertion count per timestep is 0");
   int nfinal = update->ntimestep + 1 + (ninsert-1)/nper * nfreq;
 
   // print stats
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 9dbf6e8266..03ea47671b 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -421,7 +421,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   mass_body = NULL;
   nmax_mass = 0;
 
-  staticflag = 1;
+  staticflag = 0;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
index 8e599f2886..5b60535218 100644
--- a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
+++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
@@ -759,9 +759,9 @@ double PairLJCutCoulLongSoft::init_one(int i, int j)
     double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+    etail_ij = 8.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
-    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
   }
 
diff --git a/src/USER-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp
index 22af520a2a..8ec4a002fc 100644
--- a/src/USER-FEP/pair_lj_cut_soft.cpp
+++ b/src/USER-FEP/pair_lj_cut_soft.cpp
@@ -620,9 +620,9 @@ double PairLJCutSoft::init_one(int i, int j)
     double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
     double rc6 = rc3*rc3;
     double rc9 = rc3*rc6;
-    etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+    etail_ij = 8.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
       sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
-    ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+    ptail_ij = 16.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
       sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
   }
 

From 23a756c86e7e70e20034872345a6990b623bda5e Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 17:58:47 +0000
Subject: [PATCH 08/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12735
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/compute_force_molecule.cpp | 106 ++++++++++++++++++++++++
 src/compute_force_molecule.h   |  63 ++++++++++++++
 src/compute_vcm_molecule.cpp   | 147 +++++++++++++++++++++++++++++++++
 src/compute_vcm_molecule.h     |  64 ++++++++++++++
 4 files changed, 380 insertions(+)
 create mode 100644 src/compute_force_molecule.cpp
 create mode 100644 src/compute_force_molecule.h
 create mode 100644 src/compute_vcm_molecule.cpp
 create mode 100644 src/compute_vcm_molecule.h

diff --git a/src/compute_force_molecule.cpp b/src/compute_force_molecule.cpp
new file mode 100644
index 0000000000..d05b162e3f
--- /dev/null
+++ b/src/compute_force_molecule.cpp
@@ -0,0 +1,106 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_force_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeForceMolecule::ComputeForceMolecule(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal compute force/molecule command");
+
+  if (atom->molecular == 0)
+    error->all(FLERR,"Compute force/molecule requires molecular atom style");
+
+  array_flag = 1;
+  size_array_cols = 3;
+  extarray = 0;
+
+  // setup molecule-based data
+
+  nmolecules = molecules_in_group(idlo,idhi);
+  size_array_rows = nmolecules;
+
+  memory->create(force,nmolecules,3,"force/molecule:force");
+  memory->create(forceall,nmolecules,3,"force/molecule:forceall");
+  array = forceall;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeForceMolecule::~ComputeForceMolecule()
+{
+  memory->destroy(force);
+  memory->destroy(forceall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeForceMolecule::init()
+{
+  int ntmp = molecules_in_group(idlo,idhi);
+  if (ntmp != nmolecules)
+    error->all(FLERR,"Molecule count changed in compute force/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeForceMolecule::compute_array()
+{
+  tagint imol;
+  double massone;
+  double unwrap[3];
+
+  invoked_array = update->ntimestep;
+
+  for (int i = 0; i < nmolecules; i++)
+    force[i][0] = force[i][1] = force[i][2] = 0.0;
+
+  double **f = atom->f;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      force[imol][0] += f[i][0];
+      force[imol][1] += f[i][1];
+      force[imol][2] += f[i][2];
+    }
+
+  MPI_Allreduce(&force[0][0],&forceall[0][0],3*nmolecules,
+                MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeForceMolecule::memory_usage()
+{
+  double bytes = 0;
+  if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+  return bytes;
+}
diff --git a/src/compute_force_molecule.h b/src/compute_force_molecule.h
new file mode 100644
index 0000000000..ec90e3c9fe
--- /dev/null
+++ b/src/compute_force_molecule.h
@@ -0,0 +1,63 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(force/molecule,ComputeForceMolecule)
+
+#else
+
+#ifndef LMP_COMPUTE_FORCE_MOLECULE_H
+#define LMP_COMPUTE_FORCE_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeForceMolecule : public Compute {
+ public:
+  ComputeForceMolecule(class LAMMPS *, int, char **);
+  ~ComputeForceMolecule();
+  void init();
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int nmolecules;
+  tagint idlo,idhi;
+
+  double **force,**forceall;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute force/molecule requires molecular atom style
+
+Self-explanatory.
+
+E: Molecule count changed in compute force/molecule
+
+Number of molecules must remain constant over time.
+
+*/
diff --git a/src/compute_vcm_molecule.cpp b/src/compute_vcm_molecule.cpp
new file mode 100644
index 0000000000..ab86d2ee0a
--- /dev/null
+++ b/src/compute_vcm_molecule.cpp
@@ -0,0 +1,147 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_vcm_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVCMMolecule::ComputeVCMMolecule(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal compute vcm/molecule command");
+
+  if (atom->molecular == 0)
+    error->all(FLERR,"Compute vcm/molecule requires molecular atom style");
+
+  array_flag = 1;
+  size_array_cols = 3;
+  extarray = 0;
+
+  // setup molecule-based data
+
+  nmolecules = molecules_in_group(idlo,idhi);
+  size_array_rows = nmolecules;
+
+  memory->create(massproc,nmolecules,"vcm/molecule:massproc");
+  memory->create(masstotal,nmolecules,"vcm/molecule:masstotal");
+  memory->create(vcm,nmolecules,3,"vcm/molecule:vcm");
+  memory->create(vcmall,nmolecules,3,"vcm/molecule:vcmall");
+  array = vcmall;
+
+  // compute masstotal for each molecule
+
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  tagint imol;
+  double massone;
+
+  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      massproc[imol] += massone;
+    }
+
+  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVCMMolecule::~ComputeVCMMolecule()
+{
+  memory->destroy(massproc);
+  memory->destroy(masstotal);
+  memory->destroy(vcm);
+  memory->destroy(vcmall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVCMMolecule::init()
+{
+  int ntmp = molecules_in_group(idlo,idhi);
+  if (ntmp != nmolecules)
+    error->all(FLERR,"Molecule count changed in compute vcm/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVCMMolecule::compute_array()
+{
+  tagint imol;
+  double massone;
+  double unwrap[3];
+
+  invoked_array = update->ntimestep;
+
+  for (int i = 0; i < nmolecules; i++)
+    vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
+
+  double **v = atom->v;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int *type = atom->type;
+  imageint *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      vcm[imol][0] += v[i][0] * massone;
+      vcm[imol][1] += v[i][1] * massone;
+      vcm[imol][2] += v[i][2] * massone;
+    }
+
+  MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nmolecules,
+                MPI_DOUBLE,MPI_SUM,world);
+  for (int i = 0; i < nmolecules; i++) {
+    vcmall[i][0] /= masstotal[i];
+    vcmall[i][1] /= masstotal[i];
+    vcmall[i][2] /= masstotal[i];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeVCMMolecule::memory_usage()
+{
+  double bytes = (bigint) nmolecules * 2 * sizeof(double);
+  if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+  return bytes;
+}
diff --git a/src/compute_vcm_molecule.h b/src/compute_vcm_molecule.h
new file mode 100644
index 0000000000..b5bce320ca
--- /dev/null
+++ b/src/compute_vcm_molecule.h
@@ -0,0 +1,64 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(vcm/molecule,ComputeVCMMolecule)
+
+#else
+
+#ifndef LMP_COMPUTE_VCM_MOLECULE_H
+#define LMP_COMPUTE_VCM_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeVCMMolecule : public Compute {
+ public:
+  ComputeVCMMolecule(class LAMMPS *, int, char **);
+  ~ComputeVCMMolecule();
+  void init();
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int nmolecules;
+  tagint idlo,idhi;
+
+  double *massproc,*masstotal;
+  double **vcm,**vcmall;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute vcm/molecule requires molecular atom style
+
+Self-explanatory.
+
+E: Molecule count changed in compute vcm/molecule
+
+Number of molecules must remain constant over time.
+
+*/

From b6e72407abc29b8c5f17fd063c6d5dcf5440d04e Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 18:08:55 +0000
Subject: [PATCH 09/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12736
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Section_commands.html       | 16 ++++----
 doc/Section_commands.txt        |  2 +
 doc/compute_force_molecule.html | 63 ++++++++++++++++++++++++++++++++
 doc/compute_force_molecule.txt  | 58 +++++++++++++++++++++++++++++
 doc/compute_vcm_molecule.html   | 65 +++++++++++++++++++++++++++++++++
 doc/compute_vcm_molecule.txt    | 60 ++++++++++++++++++++++++++++++
 6 files changed, 256 insertions(+), 8 deletions(-)
 create mode 100644 doc/compute_force_molecule.html
 create mode 100644 doc/compute_force_molecule.txt
 create mode 100644 doc/compute_vcm_molecule.html
 create mode 100644 doc/compute_vcm_molecule.txt

diff --git a/doc/Section_commands.html b/doc/Section_commands.html
index f3551e16fb..ded5b62e2f 100644
--- a/doc/Section_commands.html
+++ b/doc/Section_commands.html
@@ -443,14 +443,14 @@ KOKKOS, o = USER-OMP, t = OPT.
 <TR ALIGN="center"><TD ><A HREF = "compute_angle_local.html">angle/local</A></TD><TD ><A HREF = "compute_atom_molecule.html">atom/molecule</A></TD><TD ><A HREF = "compute_body_local.html">body/local</A></TD><TD ><A HREF = "compute_bond_local.html">bond/local</A></TD><TD ><A HREF = "compute_centro_atom.html">centro/atom</A></TD><TD ><A HREF = "compute_cluster_atom.html">cluster/atom</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_cna_atom.html">cna/atom</A></TD><TD ><A HREF = "compute_com.html">com</A></TD><TD ><A HREF = "compute_com_molecule.html">com/molecule</A></TD><TD ><A HREF = "compute_contact_atom.html">contact/atom</A></TD><TD ><A HREF = "compute_coord_atom.html">coord/atom</A></TD><TD ><A HREF = "compute_damage_atom.html">damage/atom</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "compute_dihedral_local.html">dihedral/local</A></TD><TD ><A HREF = "compute_dilatation_atom.html">dilatation/atom</A></TD><TD ><A HREF = "compute_displace_atom.html">displace/atom</A></TD><TD ><A HREF = "compute_erotate_asphere.html">erotate/asphere</A></TD><TD ><A HREF = "compute_erotate_rigid.html">erotate/rigid</A></TD><TD ><A HREF = "compute_erotate_sphere.html">erotate/sphere</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_erotate_sphere_atom.html">erotate/sphere/atom</A></TD><TD ><A HREF = "compute_event_displace.html">event/displace</A></TD><TD ><A HREF = "compute_group_group.html">group/group</A></TD><TD ><A HREF = "compute_gyration.html">gyration</A></TD><TD ><A HREF = "compute_gyration_molecule.html">gyration/molecule</A></TD><TD ><A HREF = "compute_heat_flux.html">heat/flux</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_improper_local.html">improper/local</A></TD><TD ><A HREF = "compute_inertia_molecule.html">inertia/molecule</A></TD><TD ><A HREF = "compute_ke.html">ke</A></TD><TD ><A HREF = "compute_ke_atom.html">ke/atom</A></TD><TD ><A HREF = "compute_ke_rigid.html">ke/rigid</A></TD><TD ><A HREF = "compute_msd.html">msd</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_msd_molecule.html">msd/molecule</A></TD><TD ><A HREF = "compute_msd_nongauss.html">msd/nongauss</A></TD><TD ><A HREF = "compute_pair.html">pair</A></TD><TD ><A HREF = "compute_pair_local.html">pair/local</A></TD><TD ><A HREF = "compute_pe.html">pe (c)</A></TD><TD ><A HREF = "compute_pe_atom.html">pe/atom</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_plasticity_atom.html">plasticity/atom</A></TD><TD ><A HREF = "compute_pressure.html">pressure (c)</A></TD><TD ><A HREF = "compute_property_atom.html">property/atom</A></TD><TD ><A HREF = "compute_property_local.html">property/local</A></TD><TD ><A HREF = "compute_property_molecule.html">property/molecule</A></TD><TD ><A HREF = "compute_rdf.html">rdf</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_reduce.html">reduce</A></TD><TD ><A HREF = "compute_reduce.html">reduce/region</A></TD><TD ><A HREF = "compute_slice.html">slice</A></TD><TD ><A HREF = "compute_sna.html">sna/atom</A></TD><TD ><A HREF = "compute_sna.html">snad/atom</A></TD><TD ><A HREF = "compute_sna.html">snav/atom</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_stress_atom.html">stress/atom</A></TD><TD ><A HREF = "compute_temp.html">temp (c)</A></TD><TD ><A HREF = "compute_temp_asphere.html">temp/asphere</A></TD><TD ><A HREF = "compute_temp_com.html">temp/com</A></TD><TD ><A HREF = "compute_temp_deform.html">temp/deform</A></TD><TD ><A HREF = "compute_temp_partial.html">temp/partial (c)</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_temp_profile.html">temp/profile</A></TD><TD ><A HREF = "compute_temp_ramp.html">temp/ramp</A></TD><TD ><A HREF = "compute_temp_region.html">temp/region</A></TD><TD ><A HREF = "compute_temp_sphere.html">temp/sphere</A></TD><TD ><A HREF = "compute_ti.html">ti</A></TD><TD ><A HREF = "compute_vacf.html">vacf</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "compute_voronoi_atom.html">voronoi/atom</A> 
+<TR ALIGN="center"><TD ><A HREF = "compute_erotate_sphere_atom.html">erotate/sphere/atom</A></TD><TD ><A HREF = "compute_event_displace.html">event/displace</A></TD><TD ><A HREF = "compute_force_molecule.html">force/molecule</A></TD><TD ><A HREF = "compute_group_group.html">group/group</A></TD><TD ><A HREF = "compute_gyration.html">gyration</A></TD><TD ><A HREF = "compute_gyration_molecule.html">gyration/molecule</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_heat_flux.html">heat/flux</A></TD><TD ><A HREF = "compute_improper_local.html">improper/local</A></TD><TD ><A HREF = "compute_inertia_molecule.html">inertia/molecule</A></TD><TD ><A HREF = "compute_ke.html">ke</A></TD><TD ><A HREF = "compute_ke_atom.html">ke/atom</A></TD><TD ><A HREF = "compute_ke_rigid.html">ke/rigid</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_msd.html">msd</A></TD><TD ><A HREF = "compute_msd_molecule.html">msd/molecule</A></TD><TD ><A HREF = "compute_msd_nongauss.html">msd/nongauss</A></TD><TD ><A HREF = "compute_pair.html">pair</A></TD><TD ><A HREF = "compute_pair_local.html">pair/local</A></TD><TD ><A HREF = "compute_pe.html">pe (c)</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_pe_atom.html">pe/atom</A></TD><TD ><A HREF = "compute_plasticity_atom.html">plasticity/atom</A></TD><TD ><A HREF = "compute_pressure.html">pressure (c)</A></TD><TD ><A HREF = "compute_property_atom.html">property/atom</A></TD><TD ><A HREF = "compute_property_local.html">property/local</A></TD><TD ><A HREF = "compute_property_molecule.html">property/molecule</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_rdf.html">rdf</A></TD><TD ><A HREF = "compute_reduce.html">reduce</A></TD><TD ><A HREF = "compute_reduce.html">reduce/region</A></TD><TD ><A HREF = "compute_slice.html">slice</A></TD><TD ><A HREF = "compute_sna.html">sna/atom</A></TD><TD ><A HREF = "compute_sna.html">snad/atom</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_sna.html">snav/atom</A></TD><TD ><A HREF = "compute_stress_atom.html">stress/atom</A></TD><TD ><A HREF = "compute_temp.html">temp (c)</A></TD><TD ><A HREF = "compute_temp_asphere.html">temp/asphere</A></TD><TD ><A HREF = "compute_temp_com.html">temp/com</A></TD><TD ><A HREF = "compute_temp_deform.html">temp/deform</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_temp_partial.html">temp/partial (c)</A></TD><TD ><A HREF = "compute_temp_profile.html">temp/profile</A></TD><TD ><A HREF = "compute_temp_ramp.html">temp/ramp</A></TD><TD ><A HREF = "compute_temp_region.html">temp/region</A></TD><TD ><A HREF = "compute_temp_sphere.html">temp/sphere</A></TD><TD ><A HREF = "compute_ti.html">ti</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "compute_vacf.html">vacf</A></TD><TD ><A HREF = "compute_vcm_molecule.html">vcm/molecule</A></TD><TD ><A HREF = "compute_voronoi_atom.html">voronoi/atom</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are additional compute styles in USER packages, which can be
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index 24e420484e..33000c9b96 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -639,6 +639,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "erotate/sphere"_compute_erotate_sphere.html,
 "erotate/sphere/atom"_compute_erotate_sphere_atom.html,
 "event/displace"_compute_event_displace.html,
+"force/molecule"_compute_force_molecule.html,
 "group/group"_compute_group_group.html,
 "gyration"_compute_gyration.html,
 "gyration/molecule"_compute_gyration_molecule.html,
@@ -679,6 +680,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "temp/sphere"_compute_temp_sphere.html,
 "ti"_compute_ti.html,
 "vacf"_compute_vacf.html,
+"vcm/molecule"_compute_vcm_molecule.html,
 "voronoi/atom"_compute_voronoi_atom.html :tb(c=6,ea=c)
 
 These are additional compute styles in USER packages, which can be
diff --git a/doc/compute_force_molecule.html b/doc/compute_force_molecule.html
new file mode 100644
index 0000000000..0ef54b3feb
--- /dev/null
+++ b/doc/compute_force_molecule.html
@@ -0,0 +1,63 @@
+<HTML>
+<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
+</CENTER>
+
+
+
+
+
+
+<HR>
+
+<H3>compute force/molecule command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>compute ID group-ID force/molecule 
+</PRE>
+<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
+<LI>force/molecule = style name of this compute command 
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>compute 1 fluid force/molecule 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>Define a computation that calculates the total force on individual
+molecules, as a sum over the forces on atoms in the molecule.  The
+x,y,z components of the force on each molecule are computed.
+</P>
+<P>The force on a particular molecule is only computed if one or more of
+its atoms are in the specified group.  Normally all atoms in the
+molecule should be in the group, however this is not required.  LAMMPS
+will warn you if this is not the case.  Only atoms in the group
+contribute to the force on the molecule.
+</P>
+<P>The ordering of per-molecule quantities produced by this compute is
+consistent with the ordering produced by other compute commands that
+generate per-molecule datums.  Conceptually, them molecule IDs will be
+in ascending order for any molecule with one or more of its atoms in
+the specified group.
+</P>
+<P><B>Output info:</B>
+</P>
+<P>This compute calculates a global array where the number of rows =
+Nmolecules and the number of columns = 3 for the fx,fy,fz components
+of force on each molecule.  These values can be accessed by any
+command that uses global array values from a compute as input.  See
+<A HREF = "Section_howto.html#howto_15">Section_howto 15</A> for an overview of
+LAMMPS output options.
+</P>
+<P>The array values are "intensive".  The array values will be in
+distance <A HREF = "units.html">units</A>.
+</P>
+<P><B>Restrictions:</B> none
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "compute_vcm_molecule.html">compute vcm/molecule</A>
+</P>
+<P><B>Default:</B> none
+</P>
+</HTML>
diff --git a/doc/compute_force_molecule.txt b/doc/compute_force_molecule.txt
new file mode 100644
index 0000000000..5021788158
--- /dev/null
+++ b/doc/compute_force_molecule.txt
@@ -0,0 +1,58 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+compute force/molecule command :h3
+
+[Syntax:]
+
+compute ID group-ID force/molecule :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+force/molecule = style name of this compute command :ul
+
+[Examples:]
+
+compute 1 fluid force/molecule :pre
+
+[Description:]
+
+Define a computation that calculates the total force on individual
+molecules, as a sum over the forces on atoms in the molecule.  The
+x,y,z components of the force on each molecule are computed.
+
+The force on a particular molecule is only computed if one or more of
+its atoms are in the specified group.  Normally all atoms in the
+molecule should be in the group, however this is not required.  LAMMPS
+will warn you if this is not the case.  Only atoms in the group
+contribute to the force on the molecule.
+
+The ordering of per-molecule quantities produced by this compute is
+consistent with the ordering produced by other compute commands that
+generate per-molecule datums.  Conceptually, them molecule IDs will be
+in ascending order for any molecule with one or more of its atoms in
+the specified group.
+
+[Output info:]
+
+This compute calculates a global array where the number of rows =
+Nmolecules and the number of columns = 3 for the fx,fy,fz components
+of force on each molecule.  These values can be accessed by any
+command that uses global array values from a compute as input.  See
+"Section_howto 15"_Section_howto.html#howto_15 for an overview of
+LAMMPS output options.
+
+The array values are "intensive".  The array values will be in
+distance "units"_units.html.
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute vcm/molecule"_compute_vcm_molecule.html
+
+[Default:] none
diff --git a/doc/compute_vcm_molecule.html b/doc/compute_vcm_molecule.html
new file mode 100644
index 0000000000..0bc881b1d7
--- /dev/null
+++ b/doc/compute_vcm_molecule.html
@@ -0,0 +1,65 @@
+<HTML>
+<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
+</CENTER>
+
+
+
+
+
+
+<HR>
+
+<H3>compute vcm/molecule command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>compute ID group-ID vcm/molecule 
+</PRE>
+<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
+<LI>vcm/molecule = style name of this compute command 
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>compute 1 fluid vcm/molecule 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>Define a computation that calculates the center-of-mass velocity of
+individual molecules.  The x,y,z components of the velocity of each
+molecule are computed.  This is calcualted by summing mass*velocity
+for each atom in the molecule and dividing the sum by the total mass
+of the molecule.
+</P>
+<P>The velocity of a particular molecule is only computed if one or more
+of its atoms are in the specified group.  Normally all atoms in the
+molecule should be in the group, however this is not required.  LAMMPS
+will warn you if this is not the case.  Only atoms in the group
+contribute to the velocity calculation for the molecule.
+</P>
+<P>The ordering of per-molecule quantities produced by this compute is
+consistent with the ordering produced by other compute commands that
+generate per-molecule datums.  Conceptually, them molecule IDs will be
+in ascending order for any molecule with one or more of its atoms in
+the specified group.
+</P>
+<P><B>Output info:</B>
+</P>
+<P>This compute calculates a global array where the number of rows =
+Nmolecules and the number of columns = 3 for the vx,vy,vz components
+of the center-of-mass velocity of each molecule.  These values can be
+accessed by any command that uses global array values from a compute
+as input.  See <A HREF = "Section_howto.html#howto_15">Section_howto 15</A> for an
+overview of LAMMPS output options.
+</P>
+<P>The array values are "intensive".  The array values will be in
+distance <A HREF = "units.html">units</A>.
+</P>
+<P><B>Restrictions:</B> none
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "compute_force_molecule.html">compute force/molecule</A>
+</P>
+<P><B>Default:</B> none
+</P>
+</HTML>
diff --git a/doc/compute_vcm_molecule.txt b/doc/compute_vcm_molecule.txt
new file mode 100644
index 0000000000..24a8fc5031
--- /dev/null
+++ b/doc/compute_vcm_molecule.txt
@@ -0,0 +1,60 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+compute vcm/molecule command :h3
+
+[Syntax:]
+
+compute ID group-ID vcm/molecule :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+vcm/molecule = style name of this compute command :ul
+
+[Examples:]
+
+compute 1 fluid vcm/molecule :pre
+
+[Description:]
+
+Define a computation that calculates the center-of-mass velocity of
+individual molecules.  The x,y,z components of the velocity of each
+molecule are computed.  This is calcualted by summing mass*velocity
+for each atom in the molecule and dividing the sum by the total mass
+of the molecule.
+
+The velocity of a particular molecule is only computed if one or more
+of its atoms are in the specified group.  Normally all atoms in the
+molecule should be in the group, however this is not required.  LAMMPS
+will warn you if this is not the case.  Only atoms in the group
+contribute to the velocity calculation for the molecule.
+
+The ordering of per-molecule quantities produced by this compute is
+consistent with the ordering produced by other compute commands that
+generate per-molecule datums.  Conceptually, them molecule IDs will be
+in ascending order for any molecule with one or more of its atoms in
+the specified group.
+
+[Output info:]
+
+This compute calculates a global array where the number of rows =
+Nmolecules and the number of columns = 3 for the vx,vy,vz components
+of the center-of-mass velocity of each molecule.  These values can be
+accessed by any command that uses global array values from a compute
+as input.  See "Section_howto 15"_Section_howto.html#howto_15 for an
+overview of LAMMPS output options.
+
+The array values are "intensive".  The array values will be in
+distance "units"_units.html.
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute force/molecule"_compute_force_molecule.html
+
+[Default:] none

From edb95fc5cd25bb25f32648584f404b9039fafa97 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 18:20:56 +0000
Subject: [PATCH 10/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12737
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index 0bd9cbf0da..fe32dc77b6 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "21 Nov 2014"
+#define LAMMPS_VERSION "22 Nov 2014"

From e835f4b6f1c7adcd04928440cf7842c6dc4f816e Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 18:20:56 +0000
Subject: [PATCH 11/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12738
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Manual.html | 4 ++--
 doc/Manual.txt  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/Manual.html b/doc/Manual.html
index 479053766f..9a8ac0e264 100644
--- a/doc/Manual.html
+++ b/doc/Manual.html
@@ -1,7 +1,7 @@
 <HTML>
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="21 Nov 2014 version">
+<META NAME="docnumber" CONTENT="22 Nov 2014 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@@ -22,7 +22,7 @@
 
 <CENTER><H3>LAMMPS Documentation 
 </H3></CENTER>
-<CENTER><H4>21 Nov 2014 version 
+<CENTER><H4>22 Nov 2014 version 
 </H4></CENTER>
 <H4>Version info: 
 </H4>
diff --git a/doc/Manual.txt b/doc/Manual.txt
index c7bb138f20..432aff8ee5 100644
--- a/doc/Manual.txt
+++ b/doc/Manual.txt
@@ -1,6 +1,6 @@
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="21 Nov 2014 version">
+<META NAME="docnumber" CONTENT="22 Nov 2014 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@@ -18,7 +18,7 @@
 <H1></H1>
 
 LAMMPS Documentation :c,h3
-21 Nov 2014 version :c,h4
+22 Nov 2014 version :c,h4
 
 Version info: :h4
 

From 65030577a92f12caaa3704763404be47b1496e19 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 19:10:20 +0000
Subject: [PATCH 12/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12740
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/lattice.html | 97 +++++++++++++++++++++++++++++-------------------
 doc/lattice.txt  | 97 +++++++++++++++++++++++++++++-------------------
 2 files changed, 118 insertions(+), 76 deletions(-)

diff --git a/doc/lattice.html b/doc/lattice.html
index e3dd79db8a..2e75e2157d 100644
--- a/doc/lattice.html
+++ b/doc/lattice.html
@@ -82,6 +82,15 @@ are the edge vectors of the unit cell.  This is the nomenclature for
 unit cell they determine does not have to be a "primitive cell" of
 minimum volume.
 </P>
+<P>Note that the lattice command can be used multiple times in an input
+script.  Each time it is invoked, the lattice attributes are
+re-defined and are used for all subsequent commands (that use lattice
+attributes).  For example, a sequence of lattice,
+<A HREF = "region.html">region</A>, and <A HREF = "create_atoms.html">create_atoms</A> commands
+can be repeated multiple times to build a poly-crystalline model with
+different geometric regions populated with atoms in different lattice
+orientations.
+</P>
 <HR>
 
 <P>A lattice of style <I>none</I> does not define a unit cell and basis set,
@@ -168,61 +177,73 @@ mapping it into the simulation box.  The <I>dim</I> argument is one of the
 the crystallographic direction in the lattice that you want to orient
 along that axis, specified as integers.  E.g. "orient x 2 1 0" means
 the x-axis in the simulation box will be the [210] lattice
-direction.  The 3 lattice directions you specify must be mutually
+direction, and similarly for y and z.  The 3 lattice directions you
+specify do not have to be unit vectors, but they must be mutually
 orthogonal and obey the right-hand rule, i.e. (X cross Y) points in
-the Z direction.  Note that this description is really only valid for
-orthogonal lattices.  If you are using the more general lattice style
-<I>custom</I> with non-orthogonal a1,a2,a3 vectors, then think of the 3
-<I>orient</I> options as creating a 3x3 rotation matrix which is applied to
-a1,a2,a3 to rotate the original unit cell to a new orientation in the
-simulation box.
+the Z direction.
+</P>
+<P>IMPORTANT NOTE: The preceding paragraph describing lattice directions
+is only valid for orthogonal cubic unit cells (or square in 2d).  If
+you are using a <I>hcp</I> or <I>hex</I> lattice or the more general lattice
+style <I>custom</I> with non-orthogonal a1,a2,a3 vectors, then you should
+think of the 3 <I>orient</I> vectors as creating a 3x3 rotation matrix
+which is applied to a1,a2,a3 to rotate the original unit cell to a new
+orientation in the simulation box.
 </P>
 <HR>
 
 <P>Several LAMMPS commands have the option to use distance units that are
-inferred from "lattice spacing" in the x,y,z box directions.  E.g. the
-<A HREF = "region.html">region</A> command can create a block of size 10x20x20,
-where 10 means 10 lattice spacings in the x direction.
+inferred from "lattice spacings" in the x,y,z box directions.
+E.g. the <A HREF = "region.html">region</A> command can create a block of size
+10x20x20, where 10 means 10 lattice spacings in the x direction.  
 </P>
-<P>The <I>spacing</I> option sets the 3 lattice spacings directly.  All must
-be non-zero (use 1.0 for dz in a 2d simulation).  The specified values
-are multiplied by the multiplicative factor described above that is
-associated with the scale factor.  Thus a spacing of 1.0 means one
-unit cell independent of the scale factor.  This option can be useful
-if the spacings LAMMPS computes are inconvenient to use in subsequent
-commands, which can be the case for non-orthogonal or rotated
-lattices.
+<P>IMPORTANT NOTE: Though they are called lattice spacings, all the
+commands that have a "units lattice" option, simply use the 3 values
+as scale factors on the distance units defined by the
+<A HREF = "units.html">units</A> command.  Thus if you do not like the lattice
+spacings computed by LAMMPS (e.g. for a non-orthogonal or rotated unit
+cell), you can define the 3 values to be whatever you wish, via the
+<I>spacing</I> option.
 </P>
 <P>If the <I>spacing</I> option is not specified, the lattice spacings are
 computed by LAMMPS in the following way.  A unit cell of the lattice
-is mapped into the simulation box (scaled, shifted, rotated), so that
-it now has (perhaps) a modified size and orientation.  The lattice
-spacing in X is defined as the difference between the min/max extent
-of the x coordinates of the 8 corner points of the modified unit cell.
-Similarly, the Y and Z lattice spacings are defined as the difference
-in the min/max of the y and z coordinates.
+is mapped into the simulation box (scaled and rotated), so that it now
+has (perhaps) a modified size and orientation.  The lattice spacing in
+X is defined as the difference between the min/max extent of the x
+coordinates of the 8 corner points of the modified unit cell (4 in
+2d).  Similarly, the Y and Z lattice spacings are defined as the
+difference in the min/max of the y and z coordinates.
 </P>
-<P>Note that if the unit cell is orthogonal with axis-aligned edges (not
-rotated via the <I>orient</I> keyword), then the lattice spacings in each
+<P>Note that if the unit cell is orthogonal with axis-aligned edges (no
+rotation via the <I>orient</I> keyword), then the lattice spacings in each
 dimension are simply the scale factor (described above) multiplied by
 the length of a1,a2,a3.  Thus a <I>hex</I> style lattice with a scale
 factor of 3.0 Angstroms, would have a lattice spacing of 3.0 in x and
 3*sqrt(3.0) in y.
 </P>
 <P>IMPORTANT NOTE: For non-orthogonal unit cells and/or when a rotation
-is applied via the <I>orient</I> keyword, then the lattice spacings may be
-less intuitive.  In particular, in these cases, there is no guarantee
-that the lattice spacing is an integer multiple of the periodicity of
-the lattice in that direction.  Thus, if you create an orthogonal
-periodic simulation box whose size in a dimension is a multiple of the
-lattice spacing, and then fill it with atoms via the
-<A HREF = "create_atoms.html">create_atoms</A> command, you will NOT necessarily
-create a periodic system.  I.e. atoms may overlap incorrectly at the
-faces of the simulation box.
+is applied via the <I>orient</I> keyword, then the lattice spacings
+computed by LAMMPS are typically less intuitive.  In particular, in
+these cases, there is no guarantee that a particular lattice spacing
+is an integer multiple of the periodicity of the lattice in that
+direction.  Thus, if you create an orthogonal periodic simulation box
+whose size in a dimension is a multiple of the lattice spacing, and
+then fill it with atoms via the <A HREF = "create_atoms.html">create_atoms</A>
+command, you will NOT necessarily create a periodic system.
+I.e. atoms may overlap incorrectly at the faces of the simulation box.
 </P>
-<P>Regardless of these issues, the values of the lattice spacings LAMMPS
-calculates are printed out, so their effect in commands that use the
-spacings should be decipherable.
+<P>The <I>spacing</I> option sets the 3 lattice spacings directly.  All must
+be non-zero (use 1.0 for dz in a 2d simulation).  The specified values
+are multiplied by the multiplicative factor described above that is
+associated with the scale factor.  Thus a spacing of 1.0 means one
+unit cell edge length independent of the scale factor.  As mentioned
+above, this option can be useful if the spacings LAMMPS computes are
+inconvenient to use in subsequent commands, which can be the case for
+non-orthogonal or rotated lattices.
+</P>
+<P>Note that whenever the lattice command is used, the values of the
+lattice spacings LAMMPS calculates are printed out.  Thus their effect
+in commands that use the spacings should be decipherable.
 </P>
 <HR>
 
diff --git a/doc/lattice.txt b/doc/lattice.txt
index 38f38ae58b..776998f27b 100644
--- a/doc/lattice.txt
+++ b/doc/lattice.txt
@@ -74,6 +74,15 @@ are the edge vectors of the unit cell.  This is the nomenclature for
 unit cell they determine does not have to be a "primitive cell" of
 minimum volume.
 
+Note that the lattice command can be used multiple times in an input
+script.  Each time it is invoked, the lattice attributes are
+re-defined and are used for all subsequent commands (that use lattice
+attributes).  For example, a sequence of lattice,
+"region"_region.html, and "create_atoms"_create_atoms.html commands
+can be repeated multiple times to build a poly-crystalline model with
+different geometric regions populated with atoms in different lattice
+orientations.
+
 :line
 
 A lattice of style {none} does not define a unit cell and basis set,
@@ -160,61 +169,73 @@ mapping it into the simulation box.  The {dim} argument is one of the
 the crystallographic direction in the lattice that you want to orient
 along that axis, specified as integers.  E.g. "orient x 2 1 0" means
 the x-axis in the simulation box will be the \[210\] lattice
-direction.  The 3 lattice directions you specify must be mutually
+direction, and similarly for y and z.  The 3 lattice directions you
+specify do not have to be unit vectors, but they must be mutually
 orthogonal and obey the right-hand rule, i.e. (X cross Y) points in
-the Z direction.  Note that this description is really only valid for
-orthogonal lattices.  If you are using the more general lattice style
-{custom} with non-orthogonal a1,a2,a3 vectors, then think of the 3
-{orient} options as creating a 3x3 rotation matrix which is applied to
-a1,a2,a3 to rotate the original unit cell to a new orientation in the
-simulation box.
+the Z direction.
+
+IMPORTANT NOTE: The preceding paragraph describing lattice directions
+is only valid for orthogonal cubic unit cells (or square in 2d).  If
+you are using a {hcp} or {hex} lattice or the more general lattice
+style {custom} with non-orthogonal a1,a2,a3 vectors, then you should
+think of the 3 {orient} vectors as creating a 3x3 rotation matrix
+which is applied to a1,a2,a3 to rotate the original unit cell to a new
+orientation in the simulation box.
 
 :line
 
 Several LAMMPS commands have the option to use distance units that are
-inferred from "lattice spacing" in the x,y,z box directions.  E.g. the
-"region"_region.html command can create a block of size 10x20x20,
-where 10 means 10 lattice spacings in the x direction.
+inferred from "lattice spacings" in the x,y,z box directions.
+E.g. the "region"_region.html command can create a block of size
+10x20x20, where 10 means 10 lattice spacings in the x direction.  
 
-The {spacing} option sets the 3 lattice spacings directly.  All must
-be non-zero (use 1.0 for dz in a 2d simulation).  The specified values
-are multiplied by the multiplicative factor described above that is
-associated with the scale factor.  Thus a spacing of 1.0 means one
-unit cell independent of the scale factor.  This option can be useful
-if the spacings LAMMPS computes are inconvenient to use in subsequent
-commands, which can be the case for non-orthogonal or rotated
-lattices.
+IMPORTANT NOTE: Though they are called lattice spacings, all the
+commands that have a "units lattice" option, simply use the 3 values
+as scale factors on the distance units defined by the
+"units"_units.html command.  Thus if you do not like the lattice
+spacings computed by LAMMPS (e.g. for a non-orthogonal or rotated unit
+cell), you can define the 3 values to be whatever you wish, via the
+{spacing} option.
 
 If the {spacing} option is not specified, the lattice spacings are
 computed by LAMMPS in the following way.  A unit cell of the lattice
-is mapped into the simulation box (scaled, shifted, rotated), so that
-it now has (perhaps) a modified size and orientation.  The lattice
-spacing in X is defined as the difference between the min/max extent
-of the x coordinates of the 8 corner points of the modified unit cell.
-Similarly, the Y and Z lattice spacings are defined as the difference
-in the min/max of the y and z coordinates.
+is mapped into the simulation box (scaled and rotated), so that it now
+has (perhaps) a modified size and orientation.  The lattice spacing in
+X is defined as the difference between the min/max extent of the x
+coordinates of the 8 corner points of the modified unit cell (4 in
+2d).  Similarly, the Y and Z lattice spacings are defined as the
+difference in the min/max of the y and z coordinates.
 
-Note that if the unit cell is orthogonal with axis-aligned edges (not
-rotated via the {orient} keyword), then the lattice spacings in each
+Note that if the unit cell is orthogonal with axis-aligned edges (no
+rotation via the {orient} keyword), then the lattice spacings in each
 dimension are simply the scale factor (described above) multiplied by
 the length of a1,a2,a3.  Thus a {hex} style lattice with a scale
 factor of 3.0 Angstroms, would have a lattice spacing of 3.0 in x and
 3*sqrt(3.0) in y.
 
 IMPORTANT NOTE: For non-orthogonal unit cells and/or when a rotation
-is applied via the {orient} keyword, then the lattice spacings may be
-less intuitive.  In particular, in these cases, there is no guarantee
-that the lattice spacing is an integer multiple of the periodicity of
-the lattice in that direction.  Thus, if you create an orthogonal
-periodic simulation box whose size in a dimension is a multiple of the
-lattice spacing, and then fill it with atoms via the
-"create_atoms"_create_atoms.html command, you will NOT necessarily
-create a periodic system.  I.e. atoms may overlap incorrectly at the
-faces of the simulation box.
+is applied via the {orient} keyword, then the lattice spacings
+computed by LAMMPS are typically less intuitive.  In particular, in
+these cases, there is no guarantee that a particular lattice spacing
+is an integer multiple of the periodicity of the lattice in that
+direction.  Thus, if you create an orthogonal periodic simulation box
+whose size in a dimension is a multiple of the lattice spacing, and
+then fill it with atoms via the "create_atoms"_create_atoms.html
+command, you will NOT necessarily create a periodic system.
+I.e. atoms may overlap incorrectly at the faces of the simulation box.
 
-Regardless of these issues, the values of the lattice spacings LAMMPS
-calculates are printed out, so their effect in commands that use the
-spacings should be decipherable.
+The {spacing} option sets the 3 lattice spacings directly.  All must
+be non-zero (use 1.0 for dz in a 2d simulation).  The specified values
+are multiplied by the multiplicative factor described above that is
+associated with the scale factor.  Thus a spacing of 1.0 means one
+unit cell edge length independent of the scale factor.  As mentioned
+above, this option can be useful if the spacings LAMMPS computes are
+inconvenient to use in subsequent commands, which can be the case for
+non-orthogonal or rotated lattices.
+
+Note that whenever the lattice command is used, the values of the
+lattice spacings LAMMPS calculates are printed out.  Thus their effect
+in commands that use the spacings should be decipherable.
 
 :line
 

From 59c9f7af480c7334f9ba77c0950c37a83aec87af Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 19:39:09 +0000
Subject: [PATCH 13/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12741
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/QEQ/fix_qeq_dynamic.cpp  | 3 ++-
 src/QEQ/fix_qeq_point.cpp    | 1 +
 src/QEQ/fix_qeq_shielded.cpp | 1 +
 3 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp
index 848ca501c3..d7501f2236 100644
--- a/src/QEQ/fix_qeq_dynamic.cpp
+++ b/src/QEQ/fix_qeq_dynamic.cpp
@@ -213,7 +213,8 @@ double FixQEqDynamic::compute_eneg()
 
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
-
+        j &= NEIGHMASK;
+      
         delr[0] = x[i][0] - x[j][0];
         delr[1] = x[i][1] - x[j][1];
         delr[2] = x[i][2] - x[j][2];
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index 5f25332c44..f6f9cbe005 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -146,6 +146,7 @@ void FixQEqPoint::compute_H()
 
       for( jj = 0; jj < jnum; jj++ ) {
         j = jlist[jj];
+        j &= NEIGHMASK;
 
         dx = x[j][0] - x[i][0];
         dy = x[j][1] - x[i][1];
diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp
index a2c8611af7..ac91798f73 100644
--- a/src/QEQ/fix_qeq_shielded.cpp
+++ b/src/QEQ/fix_qeq_shielded.cpp
@@ -190,6 +190,7 @@ void FixQEqShielded::compute_H()
 
       for( jj = 0; jj < jnum; jj++ ) {
         j = jlist[jj];
+        j &= NEIGHMASK;
 
         dx = x[j][0] - x[i][0];
         dy = x[j][1] - x[i][1];

From 2d840c556ba00d5047f26f3bce0144bf7f057519 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 19:41:23 +0000
Subject: [PATCH 14/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12742
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/QEQ/fix_qeq.cpp          | 11 +++++------
 src/QEQ/fix_qeq.h            |  2 +-
 src/QEQ/fix_qeq_dynamic.cpp  |  6 ++----
 src/QEQ/fix_qeq_point.cpp    | 10 ++++------
 src/QEQ/fix_qeq_shielded.cpp |  5 +----
 5 files changed, 13 insertions(+), 21 deletions(-)

diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index 2e735ef10c..fc8d9be09c 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -262,17 +262,13 @@ void FixQEq::reallocate_matrix()
 void FixQEq::init_list(int id, NeighList *ptr)
 {
   list = ptr;
-  if (force->kspace) force->kspace->qsum_update_flag = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixQEq::setup_pre_force(int vflag)
 {
-  if (force->newton_pair == 0)
-    error->all(FLERR,"QEQ with 'newton pair off' not supported");
-
-  neighbor->build_one(list);
+  neighbor->build_one(list->index);
 
   deallocate_storage();
   allocate_storage();
@@ -408,9 +404,11 @@ void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
 {
   int i, j, itr_j;
   int nn, NN;
+  int *ilist;
 
   nn = atom->nlocal;
   NN = atom->nlocal + atom->nghost;
+  ilist = list->ilist;
 
   for( i = 0; i < nn; ++i ) {
     if (atom->mask[i] & groupbit)
@@ -693,11 +691,12 @@ void FixQEq::vector_add( double* dest, double c, double* v, int k )
   }
 }
 
+/* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 
 void FixQEq::read_file(char *file)
 {
-  int itype,ntypes;
+  int i,itype,ntypes;
   int params_per_line = 6;
   char **words = new char*[params_per_line+1];
 
diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h
index 9128d8e16f..68176fd668 100644
--- a/src/QEQ/fix_qeq.h
+++ b/src/QEQ/fix_qeq.h
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp
index d7501f2236..ce0000576f 100644
--- a/src/QEQ/fix_qeq_dynamic.cpp
+++ b/src/QEQ/fix_qeq_dynamic.cpp
@@ -30,7 +30,6 @@
 #include "force.h"
 #include "group.h"
 #include "pair.h"
-#include "kspace.h"
 #include "respa.h"
 #include "memory.h"
 #include "error.h"
@@ -168,7 +167,6 @@ void FixQEqDynamic::pre_force(int vflag)
     }
   }
 
-  if (force->kspace) force->kspace->setup();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -213,8 +211,8 @@ double FixQEqDynamic::compute_eneg()
 
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
-        j &= NEIGHMASK;
-      
+	j &= NEIGHMASK;
+
         delr[0] = x[i][0] - x[j][0];
         delr[1] = x[i][1] - x[j][1];
         delr[2] = x[i][2] - x[j][2];
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index f6f9cbe005..3ea187b8bd 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -28,7 +28,6 @@
 #include "neigh_request.h"
 #include "update.h"
 #include "force.h"
-#include "kspace.h"
 #include "group.h"
 #include "respa.h"
 #include "memory.h"
@@ -53,8 +52,8 @@ void FixQEqPoint::init()
   int irequest = neighbor->request(this);
   neighbor->requests[irequest]->pair = 0;
   neighbor->requests[irequest]->fix  = 1;
-  neighbor->requests[irequest]->half = 1;
-  neighbor->requests[irequest]->full = 0;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
 
   int ntypes = atom->ntypes;
   memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding");
@@ -68,6 +67,7 @@ void FixQEqPoint::init()
 
 void FixQEqPoint::pre_force(int vflag)
 {
+
   if (update->ntimestep % nevery) return;
 
   n = atom->nlocal;
@@ -83,8 +83,6 @@ void FixQEqPoint::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
-  if (force->kspace) force->kspace->setup();
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -146,7 +144,7 @@ void FixQEqPoint::compute_H()
 
       for( jj = 0; jj < jnum; jj++ ) {
         j = jlist[jj];
-        j &= NEIGHMASK;
+	j &= NEIGHMASK;
 
         dx = x[j][0] - x[i][0];
         dy = x[j][1] - x[i][1];
diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp
index ac91798f73..ab9ec47a6d 100644
--- a/src/QEQ/fix_qeq_shielded.cpp
+++ b/src/QEQ/fix_qeq_shielded.cpp
@@ -28,7 +28,6 @@
 #include "neigh_request.h"
 #include "update.h"
 #include "force.h"
-#include "kspace.h"
 #include "group.h"
 #include "respa.h"
 #include "memory.h"
@@ -126,8 +125,6 @@ void FixQEqShielded::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
-  if (force->kspace) force->kspace->setup();
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -190,7 +187,7 @@ void FixQEqShielded::compute_H()
 
       for( jj = 0; jj < jnum; jj++ ) {
         j = jlist[jj];
-        j &= NEIGHMASK;
+	j &= NEIGHMASK;
 
         dx = x[j][0] - x[i][0];
         dy = x[j][1] - x[i][1];

From 7eca5edcae587b93588dca6ce242d023377fd9d0 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 19:51:53 +0000
Subject: [PATCH 15/17] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@12743
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/QEQ/fix_qeq.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index fc8d9be09c..a02ea16e73 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -268,7 +268,7 @@ void FixQEq::init_list(int id, NeighList *ptr)
 
 void FixQEq::setup_pre_force(int vflag)
 {
-  neighbor->build_one(list->index);
+  neighbor->build_one(list);
 
   deallocate_storage();
   allocate_storage();

From 72a1dc1552f153c6f4087f428250961bb4621a64 Mon Sep 17 00:00:00 2001
From: stamoor <stamoor@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 21:05:07 +0000
Subject: [PATCH 16/17] Fixing bug in qeq

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12748 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
 src/QEQ/fix_qeq_dynamic.cpp  | 3 +++
 src/QEQ/fix_qeq_point.cpp    | 3 +++
 src/QEQ/fix_qeq_shielded.cpp | 3 +++
 3 files changed, 9 insertions(+)

diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp
index ce0000576f..6099560245 100644
--- a/src/QEQ/fix_qeq_dynamic.cpp
+++ b/src/QEQ/fix_qeq_dynamic.cpp
@@ -30,6 +30,7 @@
 #include "force.h"
 #include "group.h"
 #include "pair.h"
+#include "kspace.h"
 #include "respa.h"
 #include "memory.h"
 #include "error.h"
@@ -167,6 +168,8 @@ void FixQEqDynamic::pre_force(int vflag)
     }
   }
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index 3ea187b8bd..97007485bb 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -29,6 +29,7 @@
 #include "update.h"
 #include "force.h"
 #include "group.h"
+#include "kspace.h"
 #include "respa.h"
 #include "memory.h"
 #include "error.h"
@@ -83,6 +84,8 @@ void FixQEqPoint::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp
index ab9ec47a6d..d4b061575c 100644
--- a/src/QEQ/fix_qeq_shielded.cpp
+++ b/src/QEQ/fix_qeq_shielded.cpp
@@ -29,6 +29,7 @@
 #include "update.h"
 #include "force.h"
 #include "group.h"
+#include "kspace.h"
 #include "respa.h"
 #include "memory.h"
 #include "error.h"
@@ -125,6 +126,8 @@ void FixQEqShielded::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */

From f51b6d103b774fb2234d1dd9afd4e9261d99130a Mon Sep 17 00:00:00 2001
From: stamoor <stamoor@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 19 Nov 2014 22:43:21 +0000
Subject: [PATCH 17/17] Restoring qeq files

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12749 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
 src/QEQ/fix_qeq.cpp       | 9 +++++----
 src/QEQ/fix_qeq_point.cpp | 1 -
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index a02ea16e73..2e735ef10c 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -262,12 +262,16 @@ void FixQEq::reallocate_matrix()
 void FixQEq::init_list(int id, NeighList *ptr)
 {
   list = ptr;
+  if (force->kspace) force->kspace->qsum_update_flag = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixQEq::setup_pre_force(int vflag)
 {
+  if (force->newton_pair == 0)
+    error->all(FLERR,"QEQ with 'newton pair off' not supported");
+
   neighbor->build_one(list);
 
   deallocate_storage();
@@ -404,11 +408,9 @@ void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
 {
   int i, j, itr_j;
   int nn, NN;
-  int *ilist;
 
   nn = atom->nlocal;
   NN = atom->nlocal + atom->nghost;
-  ilist = list->ilist;
 
   for( i = 0; i < nn; ++i ) {
     if (atom->mask[i] & groupbit)
@@ -691,12 +693,11 @@ void FixQEq::vector_add( double* dest, double c, double* v, int k )
   }
 }
 
-/* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 
 void FixQEq::read_file(char *file)
 {
-  int i,itype,ntypes;
+  int itype,ntypes;
   int params_per_line = 6;
   char **words = new char*[params_per_line+1];
 
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index 97007485bb..28499b5f8d 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -68,7 +68,6 @@ void FixQEqPoint::init()
 
 void FixQEqPoint::pre_force(int vflag)
 {
-
   if (update->ntimestep % nevery) return;
 
   n = atom->nlocal;