git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14696 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute_slice.cpp

@@ -18,13 +18,15 @@
 #include "modify.h"
 #include "fix.h"
 #include "group.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 #include "force.h"
 
 using namespace LAMMPS_NS;
 
-enum{COMPUTE,FIX};
+enum{COMPUTE,FIX,VARIABLE};
 
 #define INVOKED_VECTOR 2
 #define INVOKED_ARRAY 4
@@ -55,9 +57,11 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
 
   for (int iarg = 6; iarg < narg; iarg++) {
     if (strncmp(arg[iarg],"c_",2) == 0 ||
-        strncmp(arg[iarg],"f_",2) == 0) {
+        strncmp(arg[iarg],"f_",2) == 0 ||
+        strncmp(arg[iarg],"v_",2) == 0) {
       if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
       else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+      else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
 
       int n = strlen(arg[iarg]);
       char *suffix = new char[n];
@@ -89,36 +93,53 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Compute ID for compute slice does not exist");
       if (modify->compute[icompute]->vector_flag) {
         if (argindex[i])
-          error->all(FLERR,"Compute slice compute does not calculate a global array");
+          error->all(FLERR,"Compute slice compute does not "
+                     "calculate a global array");
         if (nstop > modify->compute[icompute]->size_vector)
-          error->all(FLERR,"Compute slice compute vector is accessed out-of-range");
+          error->all(FLERR,"Compute slice compute vector is "
+                     "accessed out-of-range");
       } else if (modify->compute[icompute]->array_flag) {
         if (argindex[i] == 0)
-          error->all(FLERR,"Compute slice compute does not calculate a global vector");
+          error->all(FLERR,"Compute slice compute does not "
+                     "calculate a global vector");
         if (argindex[i] > modify->compute[icompute]->size_array_cols)
-          error->all(FLERR,"Compute slice compute array is accessed out-of-range");
+          error->all(FLERR,"Compute slice compute array is "
+                     "accessed out-of-range");
         if (nstop > modify->compute[icompute]->size_array_rows)
-          error->all(FLERR,"Compute slice compute array is accessed out-of-range");
+          error->all(FLERR,"Compute slice compute array is "
+                     "accessed out-of-range");
       } else error->all(FLERR,"Compute slice compute does not calculate "
                         "global vector or array");
+
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
         error->all(FLERR,"Fix ID for compute slice does not exist");
       if (modify->fix[ifix]->vector_flag) {
         if (argindex[i])
-          error->all(FLERR,"Compute slice fix does not calculate a global array");
+          error->all(FLERR,"Compute slice fix does not "
+                     "calculate a global array");
         if (nstop > modify->fix[ifix]->size_vector)
           error->all(FLERR,"Compute slice fix vector is accessed out-of-range");
       } else if (modify->fix[ifix]->array_flag) {
         if (argindex[i] == 0)
-          error->all(FLERR,"Compute slice fix does not calculate a global vector");
+          error->all(FLERR,"Compute slice fix does not "
+                     "calculate a global vector");
         if (argindex[i] > modify->fix[ifix]->size_array_cols)
           error->all(FLERR,"Compute slice fix array is accessed out-of-range");
         if (nstop > modify->fix[ifix]->size_array_rows)
           error->all(FLERR,"Compute slice fix array is accessed out-of-range");
       } else error->all(FLERR,"Compute slice fix does not calculate "
                         "global vector or array");
+
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+        error->all(FLERR,"Variable name for compute slice does not exist");
+      if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,"Compute slice variable is not vector-style variable");
+      if (argindex[i])
+        error->all(FLERR,"Compute slice vector variable cannot be indexed");
     }
   }
 
@@ -157,6 +178,8 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
             extlist[j++] = modify->fix[ifix]->extlist[i-1];
         }
       } else extvector = modify->fix[ifix]->extarray;
+    } else if (which[0] == VARIABLE) {
+      extvector = 0;
     }
 
   } else {
@@ -189,6 +212,8 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
         } else {
           if (modify->fix[ifix]->extarray) extarray = 1;
         }
+      } else if (which[i] == VARIABLE) {
+        // variable is always intensive, does not change extarray
       }
     }
   }
@@ -225,6 +250,11 @@ void ComputeSlice::init()
       if (ifix < 0)
         error->all(FLERR,"Fix ID for compute slice does not exist");
       value2index[m] = ifix;
+    } else if (which[m] == VARIABLE) {
+      int ivariable = input->variable->find(ids[m]);
+      if (ivariable < 0)
+        error->all(FLERR,"Variable name for compute slice does not exist");
+      value2index[m] = ivariable;
     }
   }
 }
@@ -292,7 +322,8 @@ void ComputeSlice::extract_one(int m, double *vec, int stride)
 
   } else if (which[m] == FIX) {
     if (update->ntimestep % modify->fix[value2index[m]]->global_freq)
-      error->all(FLERR,"Fix used in compute slice not computed at compatible time");
+      error->all(FLERR,"Fix used in compute slice not "
+                 "computed at compatible time");
     Fix *fix = modify->fix[value2index[m]];
 
     if (argindex[m] == 0) {
@@ -309,5 +340,18 @@ void ComputeSlice::extract_one(int m, double *vec, int stride)
         j += stride;
       }
     }
+
+    // invoke vector-style variable
+
+  } else if (which[m] == VARIABLE) {
+    double *varvec;
+    int nvec = input->variable->compute_vector(value2index[m],&varvec);
+    if (nvec < nstop) 
+      error->all(FLERR,"Compute slice variable is not long enough");
+    j = 0;
+    for (i = nstart; i < nstop; i += nskip) {
+      vec[j] = varvec[i-1];
+      j += stride;
+    }
   }
 }

src/domain.cpp

@@ -498,6 +498,17 @@ void Domain::pbc()
   int *mask = atom->mask;
   imageint *image = atom->image;
 
+  // verify owned atoms all have valid numerical coords
+  // may not if computed pairwise force between 2 atoms at same location
+
+  int flag = 0;
+  for (i = 0; i < nlocal; i++)
+    if (!ISFINITE(x[i][0]) || !ISFINITE(x[i][1]) || !ISFINITE(x[i][2]))
+      flag = 1;
+  if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable");
+
+  // setup for PBC checks
+
   if (triclinic == 0) {
     lo = boxlo;
     hi = boxhi;
@@ -508,6 +519,8 @@ void Domain::pbc()
     period = prd_lamda;
   }
 
+  // apply PBC to each owned atom
+
   for (i = 0; i < nlocal; i++) {
     if (xperiodic) {
       if (x[i][0] < lo[0]) {

src/fix_ave_correlate.cpp

@@ -216,9 +216,12 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
         error->all(FLERR,"Variable name for fix ave/correlate does not exist");
-      if (input->variable->equalstyle(ivariable) == 0)
+      if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
         error->all(FLERR,
                    "Fix ave/correlate variable is not equal-style variable");
+      if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,
+                   "Fix ave/correlate variable is not vector-style variable");
     }
   }
 
@@ -443,10 +446,19 @@ void FixAveCorrelate::end_of_step()
       else
         scalar = modify->fix[m]->compute_vector(argindex[i]-1);
 
-    // evaluate equal-style variable
+    // evaluate equal-style or vector-style variable
 
-    } else if (which[i] == VARIABLE)
+    } else if (which[i] == VARIABLE) {
-      scalar = input->variable->compute_equal(m);
+      if (argindex[i] == 0)
+        scalar = input->variable->compute_equal(m);
+      else {
+        double *varvec;
+        int nvec = input->variable->compute_vector(m,&varvec);
+        int index = argindex[i];
+        if (nvec < index) scalar = 0.0;
+        else scalar = varvec[index-1];
+      }
+    }
 
     values[lastindex][i] = scalar;
   }

src/fix_ave_histo.cpp

@@ -422,15 +422,32 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,
                    "Fix for fix ave/histo not computed at compatible time");
 
-    } else if (which[i] == VARIABLE && kind == GLOBAL) {
+    } else if (which[i] == VARIABLE && kind == GLOBAL && mode == SCALAR) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
         error->all(FLERR,"Variable name for fix ave/histo does not exist");
+      if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
+        error->all(FLERR,"Fix ave/histo variable is not equal-style variable");
+      if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,"Fix ave/histo variable is not vector-style variable");
+
+    } else if (which[i] == VARIABLE && kind == GLOBAL && mode == VECTOR) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+        error->all(FLERR,"Variable name for fix ave/histo does not exist");
+      if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,"Fix ave/histo variable is not vector-style variable");
+      if (argindex[i]) 
+        error->all(FLERR,"Fix ave/histo variable cannot be indexed");
 
     } else if (which[i] == VARIABLE && kind == PERATOM) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
         error->all(FLERR,"Variable name for fix ave/histo does not exist");
+      if (argindex[i] == 0 && input->variable->atomstyle(ivariable) == 0)
+        error->all(FLERR,"Fix ave/histo variable is not atom-style variable");
+      if (argindex[i]) 
+        error->all(FLERR,"Fix ave/histo variable cannot be indexed");
     }
   }
 
@@ -719,19 +736,32 @@ void FixAveHisto::end_of_step()
                      fix->size_local_cols);
       }
 
-    // evaluate equal-style or atom-style variable
+    // evaluate equal-style or vector-style or atom-style variable
 
-    } else if (which[i] == VARIABLE && kind == GLOBAL) {
+    } else if (which[i] == VARIABLE) {
-      bin_one(input->variable->compute_equal(m));
+      if (kind == GLOBAL && mode == SCALAR) {
+        if (j == 0) bin_one(input->variable->compute_equal(m));
+        else {
+          double *varvec;
+          int nvec = input->variable->compute_vector(m,&varvec);
+          if (nvec < j) bin_one(0.0);
+          else bin_one(varvec[j-1]);
+        }
 
-    } else if (which[i] == VARIABLE && kind == PERATOM) {
+      } else if (kind == GLOBAL && mode == VECTOR) {
-      if (atom->nlocal > maxatom) {
+        double *varvec;
-        memory->destroy(vector);
+        int nvec = input->variable->compute_vector(m,&varvec);
-        maxatom = atom->nmax;
+        bin_vector(nvec,varvec,1);
-        memory->create(vector,maxatom,"ave/histo:vector");
+
+      } else if (which[i] == VARIABLE && kind == PERATOM) {
+        if (atom->nlocal > maxatom) {
+          memory->destroy(vector);
+          maxatom = atom->nmax;
+          memory->create(vector,maxatom,"ave/histo:vector");
+        }
+        input->variable->compute_atom(m,igroup,vector,1,0);
+        bin_atoms(vector,1);
       }
-      input->variable->compute_atom(m,igroup,vector,1,0);
-      bin_atoms(vector,1);
     }
   }
 

src/fix_ave_time.cpp

@@ -243,14 +243,24 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,
                    "Fix for fix ave/time not computed at compatible time");
 
-    } else if (which[i] == VARIABLE) {
+    } else if (which[i] == VARIABLE && mode == SCALAR) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
         error->all(FLERR,"Variable name for fix ave/time does not exist");
-      if (input->variable->equalstyle(ivariable) == 0)
+      if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
         error->all(FLERR,"Fix ave/time variable is not equal-style variable");
-      if (mode == VECTOR)
+      if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
-        error->all(FLERR,"Fix ave/time cannot use variable with vector mode");
+        error->all(FLERR,"Fix ave/time variable is not vector-style variable");
+
+    } else if (which[i] == VARIABLE && mode == VECTOR) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+        error->all(FLERR,"Variable name for fix ave/time does not exist");
+      if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,"Fix ave/time variable is not vector-style variable");
+      if (argindex[i])
+        error->all(FLERR,"Fix ave/time mode vector variable cannot be indexed");
+      varlen[i] = 1;
     }
   }
 
@@ -282,7 +292,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   if (any_variable_length &&
       (nrepeat > 1 || ave == RUNNING || ave == WINDOW)) {
     for (int i = 0; i < nvalues; i++)
-      if (varlen[i]) {
+      if (varlen[i] && which[i] == COMPUTE) {
         int icompute = modify->find_compute(ids[i]);
         modify->compute[icompute]->lock_enable();
       }
@@ -358,8 +368,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
         if (argindex[0] == 0) extscalar = fix->extscalar;
         else if (fix->extvector >= 0) extscalar = fix->extvector;
         else extscalar = fix->extlist[argindex[0]-1];
-      } else if (which[0] == VARIABLE)
+      } else if (which[0] == VARIABLE) {
         extscalar = 0;
+      }
 
     } else {
       vector_flag = 1;
@@ -377,8 +388,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
           if (argindex[i] == 0) extlist[i] = fix->extscalar;
           else if (fix->extvector >= 0) extlist[i] = fix->extvector;
           else extlist[i] = fix->extlist[argindex[i]-1];
-        } else if (which[i] == VARIABLE)
+        } else if (which[i] == VARIABLE) {
           extlist[i] = 0;
+        }
       }
     }
 
@@ -405,6 +417,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
             for (int i = 0; i < nrows; i++) extlist[i] = fix->extlist[i];
           }
         } else extvector = fix->extarray;
+      } else if (which[0] == VARIABLE) {
+        extlist = new int[nrows];
+        for (int i = 0; i < nrows; i++) extlist[i] = 0;
       }
 
     } else {
@@ -422,6 +437,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
           Fix *fix = modify->fix[modify->find_fix(ids[i])];
           if (argindex[i] == 0) value = fix->extvector;
           else value = fix->extarray;
+        } else if (which[i] == VARIABLE) {
+          value = 0;
         }
         if (value == -1)
           error->all(FLERR,"Fix ave/time cannot set output array "
@@ -518,13 +535,11 @@ void FixAveTime::init()
       if (icompute < 0)
         error->all(FLERR,"Compute ID for fix ave/time does not exist");
       value2index[i] = icompute;
-
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
         error->all(FLERR,"Fix ID for fix ave/time does not exist");
       value2index[i] = ifix;
-
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
@@ -576,8 +591,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
   double scalar;
 
   // zero if first sample within single Nfreq epoch
-  // NOTE: doc this
+  // if any input is variable length, initialize current length
-  // are not checking for returned length, just initialize it
   // check for exceeding length is done below
 
   if (irepeat == 0) {
@@ -599,6 +613,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
     m = value2index[i];
 
     // invoke compute if not previously invoked
+    // insure no out-of-range access to variable-length compute vector
 
     if (which[i] == COMPUTE) {
       Compute *compute = modify->compute[m];
@@ -614,9 +629,6 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
           compute->compute_vector();
           compute->invoked_flag |= INVOKED_VECTOR;
         }
-
-        // insure no out-of-range access to variable-length compute vector
-
         if (varlen[i] && compute->size_vector < argindex[i]) scalar = 0.0;
         else scalar = compute->vector[argindex[i]-1];
       }
@@ -629,10 +641,19 @@ void FixAveTime::invoke_scalar(bigint ntimestep)
       else
         scalar = modify->fix[m]->compute_vector(argindex[i]-1);
 
-    // evaluate equal-style variable
+    // evaluate equal-style or vector-style variable
+    // insure no out-of-range access to vector-style variable
 
-    } else if (which[i] == VARIABLE)
+    } else if (which[i] == VARIABLE) {
-      scalar = input->variable->compute_equal(m);
+      if (argindex[i] == 0)
+        scalar = input->variable->compute_equal(m);
+      else {
+        double *varvec;
+        int nvec = input->variable->compute_vector(m,&varvec);
+        if (nvec < argindex[i]) scalar = 0.0;
+        else scalar = varvec[argindex[i]-1];
+      }
+    }
 
     // add value to vector or just set directly if offcol is set
 
@@ -750,7 +771,7 @@ void FixAveTime::invoke_vector(bigint ntimestep)
       bigint ntimestep = update->ntimestep;
       int lockforever_flag = 0;
       for (i = 0; i < nvalues; i++) {
-        if (!varlen[i]) continue;
+        if (!varlen[i] || which[i] != COMPUTE) continue;
         if (nrepeat > 1 && ave == ONE) {
           Compute *compute = modify->compute[value2index[i]];
           compute->lock(this,ntimestep,ntimestep+(nrepeat-1)*nevery);
@@ -812,6 +833,18 @@ void FixAveTime::invoke_vector(bigint ntimestep)
         for (i = 0; i < nrows; i++)
           column[i] = fix->compute_array(i,icol);
       }
+
+    // evaluate vector-style variable
+    // insure nvec = nrows, else error
+    // could be different on this timestep than when column_length(1) set nrows
+
+    } else if (which[j] == VARIABLE) {
+      double *varvec;
+      int nvec = input->variable->compute_vector(m,&varvec);
+      if (nvec != nrows) 
+        error->all(FLERR,"Fix ave/time vector-style variable changed length");
+      for (i = 0; i < nrows; i++)
+        column[i] = varvec[i];
     }
 
     // add columns of values to array or just set directly if offcol is set
@@ -935,7 +968,9 @@ int FixAveTime::column_length(int dynamic)
         int ifix = modify->find_fix(ids[i]);
         if (argindex[i] == 0) lengthone = modify->fix[ifix]->size_vector;
         else lengthone = modify->fix[ifix]->size_array_rows;
-      }
+      } else if (which[i] == VARIABLE) {
+        // variables are always varlen = 1, so dynamic
+      } 
       if (length == 0) length = lengthone;
       else if (lengthone != length)
         error->all(FLERR,"Fix ave/time columns are inconsistent lengths");
@@ -945,15 +980,20 @@ int FixAveTime::column_length(int dynamic)
   // determine new nrows for dynamic values
   // either all must be the same
   // or must match other static values
-  // don't need to check if not MODE = VECTOR, just inovke lock_length()
+  // don't need to check if not MODE = VECTOR, just invoke lock_length()
 
   if (dynamic) {
     length = 0;
     for (int i = 0; i < nvalues; i++) {
       if (varlen[i] == 0) continue;
       m = value2index[i];
-      Compute *compute = modify->compute[m];
+      if (which[i] == COMPUTE) {
-      lengthone = compute->lock_length();
+        Compute *compute = modify->compute[m];
+        lengthone = compute->lock_length();
+      } else if (which[i] == VARIABLE) {
+        double *varvec;
+        lengthone = input->variable->compute_vector(m,&varvec);
+      }
       if (mode == SCALAR) continue;
       if (all_variable_length) {
         if (length == 0) length = lengthone;

src/fix_vector.cpp

@@ -111,8 +111,10 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
         error->all(FLERR,"Variable name for fix vector does not exist");
-      if (input->variable->equalstyle(ivariable) == 0)
+      if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
         error->all(FLERR,"Fix vector variable is not equal-style variable");
+      if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
+        error->all(FLERR,"Fix vector variable is not vector-style variable");
     }
   }
 
@@ -290,10 +292,18 @@ void FixVector::end_of_step()
       else
         result[i] = modify->fix[m]->compute_vector(argindex[i]-1);
 
-    // evaluate equal-style variable
+    // evaluate equal-style or vector-style variable
 
     } else if (which[i] == VARIABLE)
-      result[i] = input->variable->compute_equal(m);
+      if (argindex[i] == 0)
+        result[i] = input->variable->compute_equal(m);
+      else {
+        double *varvec;
+        int nvec = input->variable->compute_vector(m,&varvec);
+        int index = argindex[i];
+        if (nvec < index) result[i] = 0.0;
+        else result[i] = varvec[index-1];
+      }
   }
 
   // trigger computes on next needed step

src/thermo.cpp

@@ -911,11 +911,14 @@ void Thermo::parse_fields(char *str)
         n = input->variable->find(id);
         if (n < 0)
           error->all(FLERR,"Could not find thermo custom variable name");
-        if (input->variable->equalstyle(n) == 0)
+        if (argindex1[nfield] == 0 && input->variable->equalstyle(n) == 0)
           error->all(FLERR,
                      "Thermo custom variable is not equal-style variable");
-        if (argindex1[nfield])
+        if (argindex1[nfield] && input->variable->vectorstyle(n) == 0)
-          error->all(FLERR,"Thermo custom variable cannot be indexed");
+          error->all(FLERR,
+                     "Thermo custom variable is not vector-style variable");
+        if (argindex2[nfield])
+          error->all(FLERR,"Thermo custom variable cannot have two indices");
 
         field2index[nfield] = add_variable(id);
         addfield(copy,&Thermo::compute_variable,FLOAT);
@@ -1474,7 +1477,17 @@ void Thermo::compute_fix()
 
 void Thermo::compute_variable()
 {
-  dvalue = input->variable->compute_equal(variables[field2index[ifield]]);
+  int iarg = argindex1[ifield];
+
+  if (iarg == 0)
+    dvalue = input->variable->compute_equal(variables[field2index[ifield]]);
+  else {
+    double *varvec;
+    int nvec = 
+      input->variable->compute_vector(variables[field2index[ifield]],&varvec);
+    if (nvec < iarg) dvalue = 0.0;
+    else dvalue = varvec[iarg-1];
+  }
 }
 
 /* ----------------------------------------------------------------------

src/variable.cpp

@@ -52,7 +52,7 @@ using namespace MathConst;
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
 enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV,
-     SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,PYTHON};
+     SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,VECTOR,PYTHON};
 enum{ARG,OP};
 
 // customize by adding a function
@@ -65,7 +65,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
      RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,LOGFREQ2,
      STRIDE,STRIDE2,VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
      IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,
-     VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY};
+     VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY};
 
 // customize by adding a special function
 
@@ -92,6 +92,7 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
   pad = NULL;
   reader = NULL;
   data = NULL;
+  vecs = NULL;
 
   eval_in_progress = NULL;
 
@@ -125,6 +126,7 @@ Variable::~Variable()
     if (style[i] == LOOP || style[i] == ULOOP) delete [] data[i][0];
     else for (int j = 0; j < num[i]; j++) delete [] data[i][j];
     delete [] data[i];
+    if (style[i] == VECTOR) memory->destroy(vecs[i].values);
   }
   memory->sfree(names);
   memory->destroy(style);
@@ -133,6 +135,7 @@ Variable::~Variable()
   memory->destroy(pad);
   memory->sfree(reader);
   memory->sfree(data);
+  memory->sfree(vecs);
 
   memory->destroy(eval_in_progress);
 
@@ -428,6 +431,30 @@ void Variable::set(int narg, char **arg)
       copy(1,&arg[2],data[nvar]);
     }
 
+  // VECTOR
+  // replace pre-existing var if also style VECTOR (allows it to be reset)
+  // num = 1, which = 1st value
+  // data = 1 value, string to eval
+
+  } else if (strcmp(arg[1],"vector") == 0) {
+    if (narg != 3) error->all(FLERR,"Illegal variable command");
+    int ivar = find(arg[0]);
+    if (ivar >= 0) {
+      if (style[ivar] != VECTOR)
+        error->all(FLERR,"Cannot redefine variable as a different style");
+      delete [] data[ivar][0];
+      copy(1,&arg[2],data[ivar]);
+      replaceflag = 1;
+    } else {
+      if (nvar == maxvar) grow();
+      style[nvar] = VECTOR;
+      num[nvar] = 1;
+      which[nvar] = 0;
+      pad[nvar] = 0;
+      data[nvar] = new char*[num[nvar]];
+      copy(1,&arg[2],data[nvar]);
+    }
+
   // PYTHON
   // replace pre-existing var if also style PYTHON (allows it to be reset)
   // num = 2, which = 1st value
@@ -530,12 +557,12 @@ int Variable::next(int narg, char **arg)
       error->all(FLERR,"All variables in next command must be same style");
   }
 
-  // invalid styles: STRING, EQUAL, WORLD, ATOM, GETENV, FORMAT, PYTHON
+  // invalid styles: STRING, EQUAL, WORLD, ATOM, VECTOR, GETENV, FORMAT, PYTHON
 
   int istyle = style[find(arg[0])];
   if (istyle == STRING || istyle == EQUAL || istyle == WORLD ||
-      istyle == GETENV || istyle == ATOM || istyle == FORMAT ||
+      istyle == GETENV || istyle == ATOM || istyle == VECTOR || 
-      istyle == PYTHON)
+      istyle == FORMAT || istyle == PYTHON)
     error->all(FLERR,"Invalid variable style with next command");
 
   // if istyle = UNIVERSE or ULOOP, insure all such variables are incremented
@@ -712,6 +739,17 @@ int Variable::atomstyle(int ivar)
   return 0;
 }
 
+/* ----------------------------------------------------------------------
+   return 1 if variable is VECTOR style, 0 if not
+   this is checked before call to compute_vector() to return a vector of doubles
+------------------------------------------------------------------------- */
+
+int Variable::vectorstyle(int ivar)
+{
+  if (style[ivar] == VECTOR) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    check if variable with name is PYTHON and matches funcname
    called by Python class before it invokes a Python function
@@ -737,7 +775,7 @@ char *Variable::pythonstyle(char *name, char *funcname)
    if FORMAT, evaluate its variable and put formatted result in str
    if GETENV, query environment and put result in str
    if PYTHON, evaluate Python function, it will put result in str
-   if ATOM or ATOMFILE, return NULL
+   if ATOM or ATOMFILE or VECTOR, return NULL
    return NULL if no variable with name, or which value is bad,
      caller must respond
 ------------------------------------------------------------------------- */
@@ -797,8 +835,9 @@ char *Variable::retrieve(char *name)
       error->all(FLERR,"Python variable does not match Python function");
     python->invoke_function(ifunc,data[ivar][1]);
     str = data[ivar][1];
-  } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return NULL;
+  } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE ||
-
+	     style[ivar] == VECTOR) return NULL;
+  
   eval_in_progress[ivar] = 0;
 
   return str;
@@ -859,6 +898,7 @@ void Variable::compute_atom(int ivar, int igroup,
   eval_in_progress[ivar] = 1;
 
   if (style[ivar] == ATOM) {
+    treetype = ATOM;
     evaluate(data[ivar][0],&tree);
     collapse_tree(tree);
   } else vstore = reader[ivar]->fixstore->vstore;
@@ -908,10 +948,58 @@ void Variable::compute_atom(int ivar, int igroup,
   }
 
   if (style[ivar] == ATOM) free_tree(tree);
-
   eval_in_progress[ivar] = 0;
 }
 
+/* ----------------------------------------------------------------------
+   compute result of vector-style variable evaluation
+   return length of vector and result pointer to vector values
+     if length == 0 or -1 (mismatch), generate an error
+   if variable already computed on this timestep, just return
+   else evaulate the formula and its length, store results in VecVar entry
+------------------------------------------------------------------------- */
+
+int Variable::compute_vector(int ivar, double **result)
+{
+  Tree *tree;
+  if (vecs[ivar].currentstep == update->ntimestep) {
+    *result = vecs[ivar].values;
+    return vecs[ivar].n;
+  }
+
+  if (eval_in_progress[ivar])
+    error->all(FLERR,"Variable has circular dependency");
+  eval_in_progress[ivar] = 1;
+
+  treetype = VECTOR;
+  evaluate(data[ivar][0],&tree);
+  collapse_tree(tree);
+  int nlen = size_tree_vector(tree);
+  if (nlen == 0) error->all(FLERR,"Vector-style variable has zero length");
+  if (nlen < 0) error->all(FLERR,
+                           "Inconsistent lengths in vector-style variable");
+
+  // (re)allocate space for results if necessary
+
+  if (nlen > vecs[ivar].nmax) {
+    memory->destroy(vecs[ivar].values);
+    vecs[ivar].nmax = nlen;
+    memory->create(vecs[ivar].values,vecs[ivar].nmax,"variable:values");
+  }
+
+  vecs[ivar].n = nlen;
+  vecs[ivar].currentstep = update->ntimestep;
+  double *vec = vecs[ivar].values;
+  for (int i = 0; i < nlen; i++)
+    vec[i] = eval_tree(tree,i);
+
+  free_tree(tree);
+  eval_in_progress[ivar] = 0;
+
+  *result = vec;
+  return nlen;
+}
+
 /* ----------------------------------------------------------------------
    save copy of EQUAL style ivar formula in copy
    allocate copy here, later equal_restore() call will free it
@@ -997,6 +1085,13 @@ void Variable::grow()
 
   data = (char ***) memory->srealloc(data,maxvar*sizeof(char **),"var:data");
 
+  vecs = (VecVar *) memory->srealloc(vecs,maxvar*sizeof(VecVar),"var:vecvar");
+  for (int i = old; i < maxvar; i++) {
+    vecs[i].nmax = 0;
+    vecs[i].currentstep = -1;
+    vecs[i].values = NULL;
+  }
+
   memory->grow(eval_in_progress,maxvar,"var:eval_in_progress");
   for (int i = 0; i < maxvar; i++) eval_in_progress[i] = 0;
 }
@@ -1017,7 +1112,8 @@ void Variable::copy(int narg, char **from, char **to)
 
 /* ----------------------------------------------------------------------
    recursive evaluation of a string str
-   str is an equal-style or atom-style formula containing one or more items:
+   str is an equal-style or atom-style or vector-style formula
+     containing one or more items:
      number = 0.0, -5.45, 2.8e-4, ...
      constant = PI, version, yes, no, on, off
      thermo keyword = ke, vol, atoms, ...
@@ -1034,7 +1130,7 @@ void Variable::copy(int narg, char **from, char **to)
      variable = v_name, v_name[i]
    equal-style variables passes in tree = NULL:
      evaluate the formula, return result as a double
-   atom-style variable passes in tree = non-NULL:
+   atom-style and vector-style variables pass in tree = non-NULL:
      parse the formula but do not evaluate it
      create a parse tree and return it
 ------------------------------------------------------------------------- */
@@ -1145,11 +1241,17 @@ double Variable::evaluate(char *str, Tree **tree)
       // compute
       // ----------------
 
-      if (strncmp(word,"c_",2) == 0) {
+      if (strncmp(word,"c_",2) == 0 || strncmp(word,"C_",2) == 0) {
         if (domain->box_exist == 0)
           error->all(FLERR,
                      "Variable evaluation before simulation box is defined");
 
+	// uppercase used to force access of 
+	// global vector vs global scalar, and global array vs global vector
+
+	int lowercase = 1;
+	if (word[0] == 'C') lowercase = 0;
+
         int icompute = modify->find_compute(word+2);
         if (icompute < 0)
           error->all(FLERR,"Invalid compute ID in variable formula");
@@ -1176,9 +1278,9 @@ double Variable::evaluate(char *str, Tree **tree)
           }
         }
 
-        // c_ID = scalar from global scalar
+        // c_ID = scalar from global scalar, must be lowercase
 
-        if (nbracket == 0 && compute->scalar_flag) {
+        if (nbracket == 0 && compute->scalar_flag && lowercase) {
 
           if (update->whichflag == 0) {
             if (compute->invoked_scalar != update->ntimestep)
@@ -1199,9 +1301,9 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
-        // c_ID[i] = scalar from global vector
+        // c_ID[i] = scalar from global vector, must be lowercase
 
-        } else if (nbracket == 1 && compute->vector_flag) {
+        } else if (nbracket == 1 && compute->vector_flag && lowercase) {
 
           if (index1 > compute->size_vector &&
               compute->size_vector_variable == 0)
@@ -1228,9 +1330,9 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
-        // c_ID[i][j] = scalar from global array
+        // c_ID[i][j] = scalar from global array, must be lowercase
 
-        } else if (nbracket == 2 && compute->array_flag) {
+        } else if (nbracket == 2 && compute->array_flag && lowercase) {
 
           if (index1 > compute->size_array_rows &&
               compute->size_array_rows_variable == 0)
@@ -1260,6 +1362,68 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
+        // c_ID = vector from global vector, lowercase or uppercase
+	
+        } else if (nbracket == 0 && compute->vector_flag) {
+
+          if (tree == NULL)
+            error->all(FLERR,
+                       "Compute global vector in equal-style variable formula");
+          if (treetype == ATOM)
+            error->all(FLERR,
+                       "Compute global vector in atom-style variable formula");
+	  if (compute->size_vector == 0) 
+            error->all(FLERR,"Variable formula compute vector is zero length");
+          if (update->whichflag == 0) {
+            if (compute->invoked_vector != update->ntimestep)
+              error->all(FLERR,"Compute used in variable between runs "
+                         "is not current");
+          } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+            compute->compute_vector();
+            compute->invoked_flag |= INVOKED_VECTOR;
+          }
+
+          Tree *newtree = new Tree();
+          newtree->type = VECTORARRAY;
+          newtree->array = compute->vector;
+	  newtree->nvector = compute->size_vector;
+          newtree->nstride = 1;
+          newtree->selfalloc = 0;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
+          treestack[ntreestack++] = newtree;
+
+        // c_ID[i] = vector from global array, lowercase or uppercase
+	
+        } else if (nbracket == 1 && compute->array_flag) {
+
+          if (tree == NULL)
+            error->all(FLERR,
+                       "Compute global vector in equal-style variable formula");
+          if (treetype == ATOM)
+            error->all(FLERR,
+                       "Compute global vector in atom-style variable formula");
+	  if (compute->size_array_rows == 0) 
+            error->all(FLERR,"Variable formula compute array is zero length");
+          if (update->whichflag == 0) {
+            if (compute->invoked_array != update->ntimestep)
+              error->all(FLERR,"Compute used in variable between runs "
+                         "is not current");
+          } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+            compute->compute_array();
+            compute->invoked_flag |= INVOKED_ARRAY;
+          }
+
+          Tree *newtree = new Tree();
+          newtree->type = VECTORARRAY;
+          newtree->array = &compute->array[0][index1-1];
+	  newtree->nvector = compute->size_array_rows;
+          newtree->nstride = compute->size_array_cols;
+          newtree->selfalloc = 0;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
+          treestack[ntreestack++] = newtree;
+
         // c_ID[i] = scalar from per-atom vector
 
         } else if (nbracket == 1 && compute->peratom_flag &&
@@ -1311,6 +1475,9 @@ double Variable::evaluate(char *str, Tree **tree)
           if (tree == NULL)
             error->all(FLERR,
                        "Per-atom compute in equal-style variable formula");
+          if (treetype == VECTOR)
+            error->all(FLERR,
+                       "Per-atom compute in vector-style variable formula");
           if (update->whichflag == 0) {
             if (compute->invoked_peratom != update->ntimestep)
               error->all(FLERR,"Compute used in variable between runs "
@@ -1337,6 +1504,9 @@ double Variable::evaluate(char *str, Tree **tree)
           if (tree == NULL)
             error->all(FLERR,
                        "Per-atom compute in equal-style variable formula");
+          if (treetype == VECTOR)
+            error->all(FLERR,
+                       "Per-atom compute in vector-style variable formula");
           if (index1 > compute->size_peratom_cols)
             error->all(FLERR,"Variable formula compute array "
                        "is accessed out-of-range");
@@ -1367,11 +1537,17 @@ double Variable::evaluate(char *str, Tree **tree)
       // fix
       // ----------------
 
-      } else if (strncmp(word,"f_",2) == 0) {
+      } else if (strncmp(word,"f_",2) == 0 || strncmp(word,"F_",2) == 0) {
         if (domain->box_exist == 0)
           error->all(FLERR,
                      "Variable evaluation before simulation box is defined");
 
+	// uppercase used to force access of 
+	// global vector vs global scalar, and global array vs global vector
+
+	int lowercase = 1;
+	if (word[0] == 'F') lowercase = 0;
+
         int ifix = modify->find_fix(word+2);
         if (ifix < 0) error->all(FLERR,"Invalid fix ID in variable formula");
         Fix *fix = modify->fix[ifix];
@@ -1397,9 +1573,9 @@ double Variable::evaluate(char *str, Tree **tree)
           }
         }
 
-        // f_ID = scalar from global scalar
+        // f_ID = scalar from global scalar, must be lowercase
 
-        if (nbracket == 0 && fix->scalar_flag) {
+        if (nbracket == 0 && fix->scalar_flag && lowercase) {
 
           if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
             error->all(FLERR,"Fix in variable not computed at compatible time");
@@ -1414,9 +1590,9 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
-        // f_ID[i] = scalar from global vector
+        // f_ID[i] = scalar from global vector, must be lowercase
 
-        } else if (nbracket == 1 && fix->vector_flag) {
+        } else if (nbracket == 1 && fix->vector_flag && lowercase) {
 
           if (index1 > fix->size_vector &&
               fix->size_vector_variable == 0)
@@ -1435,9 +1611,9 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
-        // f_ID[i][j] = scalar from global array
+        // f_ID[i][j] = scalar from global array, must be lowercase
 
-        } else if (nbracket == 2 && fix->array_flag) {
+        } else if (nbracket == 2 && fix->array_flag && lowercase) {
 
           if (index1 > fix->size_array_rows &&
               fix->size_array_rows_variable == 0)
@@ -1459,6 +1635,68 @@ double Variable::evaluate(char *str, Tree **tree)
             treestack[ntreestack++] = newtree;
           } else argstack[nargstack++] = value1;
 
+        // f_ID = vector from global vector, lowercase or uppercase
+
+        } else if (nbracket == 0 && fix->vector_flag) {
+
+          if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+            error->all(FLERR,"Fix in variable not computed at compatible time");
+          if (tree == NULL)
+            error->all(FLERR,"Fix global vector in "
+                       "equal-style variable formula");
+          if (treetype == ATOM)
+            error->all(FLERR,"Fix global vector in "
+                       "atom-style variable formula");
+	  if (fix->size_vector == 0) 
+            error->all(FLERR,"Variable formula fix vector is zero length");
+
+	  int nvec = fix->size_vector;
+	  double *vec;
+	  memory->create(vec,nvec,"variable:values");
+	  for (int m = 0; m < nvec; m++)
+	    vec[m] = fix->compute_vector(m);
+
+          Tree *newtree = new Tree();
+          newtree->type = VECTORARRAY;
+          newtree->array = vec;
+	  newtree->nvector = nvec;
+          newtree->nstride = 1;
+          newtree->selfalloc = 1;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
+          treestack[ntreestack++] = newtree;
+
+        // f_ID[i] = vector from global array, lowercase or uppercase
+
+        } else if (nbracket == 1 && fix->array_flag) {
+
+          if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+            error->all(FLERR,"Fix in variable not computed at compatible time");
+          if (tree == NULL)
+            error->all(FLERR,"Fix global vector in "
+                       "equal-style variable formula");
+          if (treetype == ATOM)
+            error->all(FLERR,"Fix global vector in "
+                       "atom-style variable formula");
+	  if (fix->size_array_rows == 0) 
+            error->all(FLERR,"Variable formula fix array is zero length");
+
+	  int nvec = fix->size_array_rows;
+	  double *vec;
+	  memory->create(vec,nvec,"variable:values");
+	  for (int m = 0; m < nvec; m++)
+	    vec[m] = fix->compute_array(m,index1-1);
+
+          Tree *newtree = new Tree();
+          newtree->type = VECTORARRAY;
+          newtree->array = vec;
+	  newtree->nvector = nvec;
+          newtree->nstride = 1;
+          newtree->selfalloc = 1;
+          newtree->first = newtree->second = NULL;
+          newtree->nextra = 0;
+          treestack[ntreestack++] = newtree;
+
         // f_ID[i] = scalar from per-atom vector
 
         } else if (nbracket == 1 && fix->peratom_flag &&
@@ -1568,9 +1806,10 @@ double Variable::evaluate(char *str, Tree **tree)
           i = ptr-str+1;
         }
 
-        // v_name = scalar from non atom/atomfile variable
+        // v_name = scalar from non atom/atomfile and non vector-style variable
 
-        if (nbracket == 0 && style[ivar] != ATOM && style[ivar] != ATOMFILE) {
+        if (nbracket == 0 && style[ivar] != ATOM && style[ivar] != ATOMFILE &&
+	    style[ivar] != VECTOR) {
 
           char *var = retrieve(word+2);
           if (var == NULL)
@@ -1592,6 +1831,10 @@ double Variable::evaluate(char *str, Tree **tree)
           if (tree == NULL)
             error->all(FLERR,
                        "Atom-style variable in equal-style variable formula");
+	  if (treetype == VECTOR)
+            error->all(FLERR,
+                       "Atom-style variable in vector-style variable formula");
+
           Tree *newtree;
           evaluate(data[ivar][0],&newtree);
           treestack[ntreestack++] = newtree;
@@ -1603,6 +1846,10 @@ double Variable::evaluate(char *str, Tree **tree)
           if (tree == NULL)
             error->all(FLERR,"Atomfile-style variable in "
                        "equal-style variable formula");
+	  if (treetype == VECTOR)
+            error->all(FLERR,
+                       "Atomfile-style variable in vector-style variable formula");
+
           Tree *newtree = new Tree();
           newtree->type = ATOMARRAY;
           newtree->array = reader[ivar]->fixstore->vstore;
@@ -1612,6 +1859,31 @@ double Variable::evaluate(char *str, Tree **tree)
           newtree->nextra = 0;
           treestack[ntreestack++] = newtree;
 
+        // v_name = vector from vector-style variable
+        // evaluate the vector-style variable, put result in newtree
+
+	} else if (nbracket == 0 && style[ivar] == VECTOR) {
+
+          if (tree == NULL)
+            error->all(FLERR,
+                       "Vector-style variable in equal-style variable formula");
+          if (treetype == ATOM)
+            error->all(FLERR,
+                       "Vector-style variable in atom-style variable formula");
+
+	  double *vec;
+	  int nvec = compute_vector(ivar,&vec);
+
+	  Tree *newtree = new Tree();
+	  newtree->type = VECTORARRAY;
+	  newtree->array = vec;
+	  newtree->nvector = nvec;
+	  newtree->nstride = 1;
+	  newtree->selfalloc = 0;
+	  newtree->first = newtree->second = NULL;
+	  newtree->nextra = 0;
+	  treestack[ntreestack++] = newtree;
+
         // v_name[N] = scalar from atom-style variable
         // compute the per-atom variable in result
         // use peratom2global to extract single value from result
@@ -1632,6 +1904,26 @@ double Variable::evaluate(char *str, Tree **tree)
           peratom2global(1,NULL,reader[ivar]->fixstore->vstore,1,index,
                          tree,treestack,ntreestack,argstack,nargstack);
 
+        // v_name[N] = scalar from vector-style variable
+        // compute the vector-style variable, extract single value
+
+        } else if (nbracket && style[ivar] == VECTOR) {
+
+	  double *vec;
+	  int nvec = compute_vector(ivar,&vec);
+	  if (index <= 0 || index > nvec)
+	    error->all(FLERR,"Invalid index into vector-style variable");
+	  int m = index;   // convert from tagint to int
+
+          if (tree) {
+            Tree *newtree = new Tree();
+            newtree->type = VALUE;
+            newtree->value = vec[m-1];
+            newtree->first = newtree->second = NULL;
+            newtree->nextra = 0;
+            treestack[ntreestack++] = newtree;
+          } else argstack[nargstack++] = vec[m-1];
+
         } else error->all(FLERR,"Mismatched variable in variable formula");
 
       // ----------------
@@ -1892,7 +2184,7 @@ double Variable::evaluate(char *str, Tree **tree)
    one-time collapse of an atom-style variable parse tree
    tree was created by one-time parsing of formula string via evaluate()
    only keep tree nodes that depend on
-     ATOMARRAY, TYPEARRAY, INTARRAY, BIGINTARRAY
+     ATOMARRAY, TYPEARRAY, INTARRAY, BIGINTARRAY, VECTOR
    remainder is converted to single VALUE
    this enables optimal eval_tree loop over atoms
    customize by adding a function:
@@ -1912,6 +2204,7 @@ double Variable::collapse_tree(Tree *tree)
   if (tree->type == TYPEARRAY) return 0.0;
   if (tree->type == INTARRAY) return 0.0;
   if (tree->type == BIGINTARRAY) return 0.0;
+  if (tree->type == VECTORARRAY) return 0.0;
 
   if (tree->type == ADD) {
     arg1 = collapse_tree(tree->first);
@@ -2403,7 +2696,8 @@ double Variable::collapse_tree(Tree *tree)
 }
 
 /* ----------------------------------------------------------------------
-   evaluate an atom-style variable parse tree for atom I
+   evaluate an atom-style or vector-style variable parse tree
+   index I = atom I or vector index I
    tree was created by one-time parsing of formula string via evaulate()
    customize by adding a function:
      sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
@@ -2422,6 +2716,7 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == TYPEARRAY) return tree->array[atom->type[i]];
   if (tree->type == INTARRAY) return (double) tree->iarray[i*tree->nstride];
   if (tree->type == BIGINTARRAY) return (double) tree->barray[i*tree->nstride];
+  if (tree->type == VECTORARRAY) return tree->array[i*tree->nstride];
 
   if (tree->type == ADD)
     return eval_tree(tree->first,i) + eval_tree(tree->second,i);
@@ -2728,6 +3023,40 @@ double Variable::eval_tree(Tree *tree, int i)
   return 0.0;
 }
 
+/* ----------------------------------------------------------------------
+   scan entire tree, find size of vectors for vector-style variable
+   return N for consistent vector size
+   return 0 for no vector size, caller flags as error
+   return -1 for inconsistent vector size, caller flags as error
+------------------------------------------------------------------------- */
+
+int Variable::size_tree_vector(Tree *tree)
+{
+  int nsize = 0;
+  if (tree->type == VECTORARRAY) nsize = tree->nvector;
+  if (tree->first) nsize = compare_tree_vector(nsize,size_tree_vector(tree->first));
+  if (tree->second) nsize = compare_tree_vector(nsize,size_tree_vector(tree->second));
+  if (tree->nextra) {
+    for (int i = 0; i < tree->nextra; i++) 
+      nsize = compare_tree_vector(nsize,size_tree_vector(tree->extra[i]));
+  }
+  return nsize;
+}
+
+/* ----------------------------------------------------------------------
+   compare size of two vectors for vector-style variable
+   return positive size if same or one has no size 0
+   return -1 error if one is already error or not same positive size
+------------------------------------------------------------------------- */
+
+int Variable::compare_tree_vector(int i, int j)
+{
+  if (i < 0 || j < 0) return -1;
+  if (i == 0 || j == 0) return MAX(i,j);
+  if (i != j) return -1;
+  return i;
+}
+
 /* ---------------------------------------------------------------------- */
 
 void Variable::free_tree(Tree *tree)
@@ -2739,9 +3068,7 @@ void Variable::free_tree(Tree *tree)
     delete [] tree->extra;
   }
 
-  if (tree->type == ATOMARRAY && tree->selfalloc)
+  if (tree->selfalloc) memory->destroy(tree->array);
-    memory->destroy(tree->array);
-
   delete tree;
 }
 
@@ -3532,9 +3859,12 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
     Compute *compute = NULL;
     Fix *fix = NULL;
+    int ivar = -1;
     int index,nvec,nstride;
     char *ptr1,*ptr2;
 
+    // argument is compute
+
     if (strstr(args[0],"c_") == args[0]) {
       ptr1 = strchr(args[0],'[');
       if (ptr1) {
@@ -3574,6 +3904,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         nstride = compute->size_array_cols;
       } else error->all(FLERR,"Mismatched compute in variable formula");
 
+    // argument is fix
+
     } else if (strstr(args[0],"f_") == args[0]) {
       ptr1 = strchr(args[0],'[');
       if (ptr1) {
@@ -3600,6 +3932,28 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         nstride = fix->size_array_cols;
       } else error->all(FLERR,"Mismatched fix in variable formula");
 
+    // argument is vector-style variable
+
+    } else if (strstr(args[0],"v_") == args[0]) {
+      ptr1 = strchr(args[0],'[');
+      if (ptr1) {
+        ptr2 = ptr1;
+        index = (int) int_between_brackets(ptr2,0);
+        *ptr1 = '\0';
+      } else index = 0;
+
+      if (index) 
+	error->all(FLERR,"Invalid special function in variable formula");
+      ivar = find(&args[0][2]);
+      if (ivar < 0) error->all(FLERR,"Invalid special function in variable formula");
+      if (style[ivar] != VECTOR) 
+	error->all(FLERR,"Mis-matched special function variable in variable formula");
+      if (eval_in_progress[ivar]) error->all(FLERR,"Variable has circular dependency");
+
+      double *vec;
+      nvec = compute_vector(ivar,&vec);
+      nstride = 1;
+
     } else error->all(FLERR,"Invalid special function in variable formula");
 
     value = 0.0;
@@ -3661,6 +4015,28 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       }
     }
 
+    if (ivar >= 0) {
+      double one;
+      double *vec = vecs[ivar].values;
+      for (int i = 0; i < nvec; i++) {
+	one = vec[i];
+        if (method == SUM) value += one;
+        else if (method == XMIN) value = MIN(value,one);
+        else if (method == XMAX) value = MAX(value,one);
+        else if (method == AVE) value += one;
+        else if (method == TRAP) value += one;
+        else if (method == SLOPE) {
+          if (nvec > 1) xvalue = (double) i / (nvec-1);
+          else xvalue = 0.0;
+          sx += xvalue;
+          sy += one;
+          sxx += xvalue*xvalue;
+          sxy += xvalue*one;
+        }
+      }
+      if (method == TRAP) value -= 0.5*vec[0] + 0.5*vec[nvec-1];
+    }
+
     if (method == AVE) value /= nvec;
 
     if (method == SLOPE) {

src/variable.h

@@ -35,12 +35,14 @@ class Variable : protected Pointers {
 
   int equalstyle(int);
   int atomstyle(int);
+  int vectorstyle(int);
   char *pythonstyle(char *, char *);
 
   char *retrieve(char *);
   double compute_equal(int);
   double compute_equal(char *);
   void compute_atom(int, int, double *, int, int);
+  int compute_vector(int, double **);
   tagint int_between_brackets(char *&, int);
   double evaluate_boolean(char *);
 
@@ -63,7 +65,15 @@ class Variable : protected Pointers {
   class VarReader **reader;   // variable that reads from file
   char ***data;            // str value of each variable's values
 
-  int *eval_in_progress;   // flag if evaluation of variable is in progress
+  struct VecVar {
+    int n,nmax;
+    bigint currentstep;
+    double *values;
+  };
+  VecVar *vecs;
+
+  int *eval_in_progress;       // flag if evaluation of variable is in progress
+  int treetype;                // ATOM or VECTOR flag for formula evaluation
 
   class RanMars *randomequal;   // random number generator for equal-style vars
   class RanMars *randomatom;    // random number generator for atom-style vars
@@ -73,12 +83,13 @@ class Variable : protected Pointers {
 
   class Python *python;    // ptr to embedded Python interpreter
 
-  struct Tree {            // parse tree for atom-style variables
+  struct Tree {            // parse tree for atom-style or vector-style variables
     double value;          // single scalar
     double *array;         // per-atom or per-type list of doubles
     int *iarray;           // per-atom list of ints
     bigint *barray;        // per-atom list of bigints
     int type;              // operation, see enum{} in variable.cpp
+    int nvector;           // length of array for vector-style variable
     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
@@ -94,6 +105,8 @@ class Variable : protected Pointers {
   double evaluate(char *, Tree **);
   double collapse_tree(Tree *);
   double eval_tree(Tree *, int);
+  int size_tree_vector(Tree *);
+  int compare_tree_vector(int, int);
   void free_tree(Tree *);
   int find_matching_paren(char *, int, char *&);
   int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);