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

src/ASPHERE/compute_temp_asphere.cpp

@@ -31,8 +31,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
 /* ---------------------------------------------------------------------- */
 
 ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :

src/compute.h

@@ -25,7 +25,7 @@ class Compute : protected Pointers {
 
   double scalar;            // computed global scalar
   double *vector;           // computed global vector
-  double *array;            // computed global array
+  double **array;           // computed global array
   double *vector_atom;      // computed per-atom vector
   double **array_atom;      // computed per-atom array
   double *vector_local;     // computed local vector
@@ -48,7 +48,7 @@ class Compute : protected Pointers {
   int extscalar;            // 0/1 if global scalar is intensive/extensive
   int extvector;            // 0/1/-1 if global vector is all int/ext/extlist
   int *extlist;             // list of 0/1 int/ext for each vec component
-  int extarray;             // 0/1 if global array is intensive/extensive
+  int extarray;             // 0/1 if global array is all intensive/extensive
 
   int tempflag;       // 1 if Compute can be used as temperature
                       // must have both compute_scalar, compute_vector

src/compute_heat_flux.cpp

@@ -35,7 +35,7 @@
 
 using namespace LAMMPS_NS;
 
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+#define INVOKED_PERATOM 8
 
 /* ---------------------------------------------------------------------- */
 
@@ -67,6 +67,7 @@ ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeHeatFlux::~ComputeHeatFlux()
 { 
+  delete [] id_atomPE;
   delete [] vector;
 }
 

src/compute_pressure.cpp

@@ -31,8 +31,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
 /* ---------------------------------------------------------------------- */
 
 ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :

src/compute_reduce.cpp

@@ -21,6 +21,7 @@
 #include "fix.h"
 #include "force.h"
 #include "comm.h"
+#include "group.h"
 #include "input.h"
 #include "variable.h"
 #include "memory.h"
@@ -28,9 +29,14 @@
 
 using namespace LAMMPS_NS;
 
-enum{SUM,MINN,MAXX};
+enum{SUM,MINN,MAXX,AVE};
 enum{X,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{GLOBAL,PERATOM,LOCAL};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+#define INVOKED_LOCAL 16
 
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
@@ -55,6 +61,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
   if (strcmp(arg[iarg],"sum") == 0) mode = SUM;
   else if (strcmp(arg[iarg],"min") == 0) mode = MINN;
   else if (strcmp(arg[iarg],"max") == 0) mode = MAXX;
+  else if (strcmp(arg[iarg],"ave") == 0) mode = AVE;
   else error->all("Illegal compute reduce command");
   iarg++;
 
@@ -62,6 +69,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
 
   which = new int[narg-4];
   argindex = new int[narg-4];
+  flavor = new int[narg-4];
   ids = new char*[narg-4];
   value2index = new int[narg-4];
   nvalues = 0;
@@ -136,9 +144,21 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
       int icompute = modify->find_compute(ids[i]);
       if (icompute < 0)
 	error->all("Compute ID for compute reduce does not exist");
-      if (modify->compute[icompute]->peratom_flag == 0)
+      if (modify->compute[icompute]->vector_flag || 
+	  modify->compute[icompute]->array_flag) {
+	flavor[i] = GLOBAL;
+	if (argindex[i] == 0 && 
+	    modify->compute[icompute]->vector_flag == 0)
 	  error->all("Compute reduce compute does not "
-		   "calculate per-atom values");
+		     "calculate a global vector");
+	if (argindex[i] && modify->compute[icompute]->array_flag == 0)
+	  error->all("Compute reduce compute does not "
+		     "calculate a global array");
+	if (argindex[i] && 
+	    argindex[i] > modify->compute[icompute]->size_array_cols)
+	  error->all("Compute reduce compute array is accessed out-of-range");
+      } else if (modify->compute[icompute]->peratom_flag) {
+	flavor[i] = PERATOM;
 	if (argindex[i] == 0 && 
 	    modify->compute[icompute]->size_peratom_cols != 0)
 	  error->all("Compute reduce compute does not "
@@ -146,16 +166,66 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
 	if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
 	  error->all("Compute reduce compute does not "
 		     "calculate a per-atom array");
+	if (argindex[i] && 
+	    argindex[i] > modify->compute[icompute]->size_peratom_cols)
+	  error->all("Compute reduce compute array is accessed out-of-range");
+      } else if (modify->compute[icompute]->local_flag) {
+	flavor[i] = LOCAL;
+	if (argindex[i] == 0 && 
+	    modify->compute[icompute]->size_local_cols != 0)
+	  error->all("Compute reduce compute does not "
+		     "calculate a local vector");
+	if (argindex[i] && modify->compute[icompute]->size_local_cols == 0)
+	  error->all("Compute reduce compute does not "
+		     "calculate a local array");
+	if (argindex[i] && 
+	    argindex[i] > modify->compute[icompute]->size_local_cols)
+	  error->all("Compute reduce compute array is accessed out-of-range");
+      }
+
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
 	error->all("Fix ID for compute reduce does not exist");
-      if (modify->fix[ifix]->peratom_flag == 0)
+      if (modify->fix[ifix]->vector_flag || 
-	error->all("Compute reduce fix does not calculate per-atom values");
+	  modify->fix[ifix]->array_flag) {
-      if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
+	flavor[i] = GLOBAL;
-	error->all("Compute reduce fix does not calculate a per-atom vector");
+	if (argindex[i] == 0 && 
+	    modify->fix[ifix]->vector_flag == 0)
+	  error->all("Compute reduce fix does not "
+		     "calculate a global vector");
+	if (argindex[i] && modify->fix[ifix]->array_flag == 0)
+	  error->all("Compute reduce fix does not "
+		     "calculate a global array");
+	if (argindex[i] && 
+	    argindex[i] > modify->fix[ifix]->size_array_cols)
+	  error->all("Compute reduce fix array is accessed out-of-range");
+      } else if (modify->fix[ifix]->peratom_flag) {
+	flavor[i] = PERATOM;
+	if (argindex[i] == 0 && 
+	    modify->fix[ifix]->size_peratom_cols != 0)
+	  error->all("Compute reduce fix does not "
+		     "calculate a per-atom vector");
 	if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
-	error->all("Compute reduce fix does not calculate a per-atom array");
+	  error->all("Compute reduce fix does not "
+		     "calculate a per-atom array");
+	if (argindex[i] && 
+	    argindex[i] > modify->fix[ifix]->size_peratom_cols)
+	  error->all("Compute reduce fix array is accessed out-of-range");
+      } else if (modify->fix[ifix]->local_flag) {
+	flavor[i] = LOCAL;
+	if (argindex[i] == 0 && 
+	    modify->fix[ifix]->size_local_cols != 0)
+	  error->all("Compute reduce fix does not "
+		     "calculate a local vector");
+	if (argindex[i] && modify->fix[ifix]->size_local_cols == 0)
+	  error->all("Compute reduce fix does not "
+		     "calculate a local array");
+	if (argindex[i] && 
+	    argindex[i] > modify->fix[ifix]->size_local_cols)
+	  error->all("Compute reduce fix array is accessed out-of-range");
+      }
+
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
@@ -192,6 +262,7 @@ ComputeReduce::~ComputeReduce()
 {
   delete [] which;
   delete [] argindex;
+  delete [] flavor;
   for (int m = 0; m < nvalues; m++) delete [] ids[m];
   delete [] ids;
   delete [] value2index;
@@ -238,12 +309,17 @@ double ComputeReduce::compute_scalar()
   invoked_scalar = update->ntimestep;
 
   double one = compute_one(0);
+
   if (mode == SUM)
     MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   else if (mode == MINN)
     MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
   else if (mode == MAXX)
     MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
+  else if (mode == AVE) {
+    MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+    scalar /= count(0);
+  }
   return scalar;
 }
 
@@ -254,12 +330,17 @@ void ComputeReduce::compute_vector()
   invoked_vector = update->ntimestep;
 
   for (int m = 0; m < nvalues; m++) onevec[m] = compute_one(m);
+
   if (mode == SUM)
     MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
   else if (mode == MINN)
     MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MIN,world);
   else if (mode == MAXX)
     MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MAX,world);
+  else if (mode == AVE) {
+    MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MAX,world);
+    for (int m = 0; m < nvalues; m++) vector[m] /= count(m);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -269,15 +350,15 @@ double ComputeReduce::compute_one(int m)
   int i;
 
   // invoke the appropriate attribute,compute,fix,variable
-  // compute scalar quantity by summing over atom scalars
+  // compute scalar quantity by summing over atom properties
-  // only include atoms in group
+  // only include atoms in group for atom properties and per-atom quantities
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
 
   int n = value2index[m];
   int j = argindex[m];
 
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
   double one;
   if (mode == SUM) one = 0.0;
   else if (mode == MINN) one = BIG;
@@ -300,6 +381,30 @@ double ComputeReduce::compute_one(int m)
 
   } else if (which[m] == COMPUTE) {
     Compute *compute = modify->compute[n];
+    
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) {
+	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+	  compute->compute_vector();
+	  compute->invoked_flag |= INVOKED_VECTOR;
+	}
+	double *compute_vector = compute->vector;
+	int n = compute->size_vector;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_vector[i]);
+      } else {
+	if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+	  compute->compute_array();
+	  compute->invoked_flag |= INVOKED_ARRAY;
+	}
+	double **compute_array = compute->array;
+	int n = compute->size_array_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_array[i][jm1]);
+      }
+
+    } else if (flavor[m] == PERATOM) {
       if (!(compute->invoked_flag & INVOKED_PERATOM)) {
 	compute->compute_peratom();
 	compute->invoked_flag |= INVOKED_PERATOM;
@@ -307,32 +412,85 @@ double ComputeReduce::compute_one(int m)
 
       if (j == 0) {
 	double *compute_vector = compute->vector_atom;
-      for (i = 0; i < nlocal; i++)
+	int n = nlocal;
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit) combine(one,compute_vector[i]);
       } else {
 	double **compute_array = compute->array_atom;
+	int n = nlocal;
 	int jm1 = j - 1;
-      for (i = 0; i < nlocal; i++)
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit) combine(one,compute_array[i][jm1]);
       }
 
-  // access fix fields, check if frequency is a match
+    } else if (flavor[m] == LOCAL) {
+      if (!(compute->invoked_flag & INVOKED_LOCAL)) {
+	compute->compute_local();
+	compute->invoked_flag |= INVOKED_LOCAL;
+      }
+
+      if (j == 0) {
+	double *compute_vector = compute->vector_local;
+	int n = compute->size_local_rows;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_vector[i]);
+      } else {
+	double **compute_array = compute->array_local;
+	int n = compute->size_local_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_array[i][jm1]);
+      }
+    }
+
+  // access fix fields, check if fix frequency is a match
 
   } else if (which[m] == FIX) {
     if (update->ntimestep % modify->fix[n]->peratom_freq)
       error->all("Fix used in compute reduce not computed at compatible time");
+    Fix *fix = modify->fix[n];
 
+    if (flavor[m] == GLOBAL) {
       if (j == 0) {
-      double *fix_vector = modify->fix[n]->vector_atom;
+	int n = fix->size_vector;
+	for (i = 0; i < n; i++)
+	  combine(one,fix->compute_vector(i));
+      } else {
+	int n = fix->size_array_rows;
+	int jm1 = j - 1;
 	for (i = 0; i < nlocal; i++)
+	  combine(one,fix->compute_array(i,jm1));
+      }
+
+    } else if (flavor[m] == PERATOM) {
+      if (j == 0) {
+	double *fix_vector = fix->vector_atom;
+	int n = nlocal;
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit) combine(one,fix_vector[i]);
       } else {
-      double **fix_array = modify->fix[n]->array_atom;
+	double **fix_array = fix->array_atom;
+	int n = nlocal;
 	int jm1 = j - 1;
 	for (i = 0; i < nlocal; i++)
 	  if (mask[i] & groupbit) combine(one,fix_array[i][jm1]);
       }
 
+    } else if (flavor[m] == LOCAL) {
+      if (j == 0) {
+	double *fix_vector = fix->vector_local;
+	int n = fix->size_local_rows;
+	for (i = 0; i < n; i++)
+	  combine(one,fix_vector[i]);
+      } else {
+	double **fix_array = fix->array_local;
+	int n = fix->size_local_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,fix_array[i][jm1]);
+      }
+    }
+    
   // evaluate atom-style variable
 
   } else if (which[m] == VARIABLE) {
@@ -351,13 +509,57 @@ double ComputeReduce::compute_one(int m)
   return one;
 }
 
+/* ---------------------------------------------------------------------- */
+
+double ComputeReduce::count(int m)
+{
+  int n = value2index[m];
+  int j = argindex[m];
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (which[m] == X || which[m] == V || which[m] == F)
+    return group->count(igroup);
+  else if (which[m] == COMPUTE) {
+    Compute *compute = modify->compute[n];
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) return(1.0*compute->size_vector);
+      else return(1.0*compute->size_array_rows);
+    } else if (flavor[m] == PERATOM) {
+      return group->count(igroup);
+    } else if (flavor[m] == LOCAL) {
+      double ncount = compute->size_local_rows;
+      double ncountall;
+      MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+      return ncountall;
+    }
+  } else if (which[m] == FIX) {
+    Fix *fix = modify->fix[n];
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) return(1.0*fix->size_vector);
+      else return(1.0*fix->size_array_rows);
+    } else if (flavor[m] == PERATOM) {
+      return group->count(igroup);
+    } else if (flavor[m] == LOCAL) {
+      double ncount = fix->size_local_rows;
+      double ncountall;
+      MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+      return ncountall;
+    }
+  } else if (which[m] == VARIABLE)
+    return group->count(igroup);
+
+  return 0.0;
+}
+
 /* ----------------------------------------------------------------------
    combine two values according to reduction mode
 ------------------------------------------------------------------------- */
 
 void ComputeReduce::combine(double &one, double two)
 {
-  if (mode == SUM) one += two;
+  if (mode == SUM || mode == AVE) one += two;
   else if (mode == MINN) one = MIN(one,two);
   else if (mode == MAXX) one = MAX(one,two);
 }

src/compute_reduce.h

@@ -29,7 +29,7 @@ class ComputeReduce : public Compute {
 
  protected:
   int mode,nvalues,iregion;
-  int *which,*argindex,*value2index;
+  int *which,*argindex,*flavor,*value2index;
   char **ids;
   double *onevec;
 
@@ -37,6 +37,7 @@ class ComputeReduce : public Compute {
   double *varatom;
 
   virtual double compute_one(int);
+  virtual double count(int);
   void combine(double &, double);
 };
 

src/compute_reduce_region.cpp

@@ -18,6 +18,7 @@
 #include "update.h"
 #include "modify.h"
 #include "domain.h"
+#include "group.h"
 #include "region.h"
 #include "fix.h"
 #include "input.h"
@@ -29,7 +30,12 @@ using namespace LAMMPS_NS;
 
 enum{SUM,MINN,MAXX};
 enum{X,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{GLOBAL,PERATOM,LOCAL};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+#define INVOKED_LOCAL 16
 
 #define BIG 1.0e20
 
@@ -81,6 +87,30 @@ double ComputeReduceRegion::compute_one(int m)
 
   } else if (which[m] == COMPUTE) {
     Compute *compute = modify->compute[n];
+    
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) {
+	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+	  compute->compute_vector();
+	  compute->invoked_flag |= INVOKED_VECTOR;
+	}
+	double *compute_vector = compute->vector;
+	int n = compute->size_vector;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_vector[i]);
+      } else {
+	if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+	  compute->compute_array();
+	  compute->invoked_flag |= INVOKED_ARRAY;
+	}
+	double **compute_array = compute->array;
+	int n = compute->size_array_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_array[i][jm1]);
+      }
+
+    } else if (flavor[m] == PERATOM) {
       if (!(compute->invoked_flag & INVOKED_PERATOM)) {
 	compute->compute_peratom();
 	compute->invoked_flag |= INVOKED_PERATOM;
@@ -88,36 +118,89 @@ double ComputeReduceRegion::compute_one(int m)
 
       if (j == 0) {
 	double *compute_vector = compute->vector_atom;
-      for (i = 0; i < nlocal; i++)
+	int n = nlocal;
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
 	    combine(one,compute_vector[i]);
       } else {
 	double **compute_array = compute->array_atom;
+	int n = nlocal;
 	int jm1 = j - 1;
-      for (i = 0; i < nlocal; i++)
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
 	    combine(one,compute_array[i][jm1]);
       }
 
-  // access fix fields, check if frequency is a match
+    } else if (flavor[m] == LOCAL) {
+      if (!(compute->invoked_flag & INVOKED_LOCAL)) {
+	compute->compute_local();
+	compute->invoked_flag |= INVOKED_LOCAL;
+      }
+
+      if (j == 0) {
+	double *compute_vector = compute->vector_local;
+	int n = compute->size_local_rows;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_vector[i]);
+      } else {
+	double **compute_array = compute->array_local;
+	int n = compute->size_local_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,compute_array[i][jm1]);
+      }
+    }
+
+  // check if fix frequency is a match
 
   } else if (which[m] == FIX) {
     if (update->ntimestep % modify->fix[n]->peratom_freq)
       error->all("Fix used in compute reduce not computed at compatible time");
+    Fix *fix = modify->fix[n];
 
+    if (flavor[m] == GLOBAL) {
       if (j == 0) {
-      double *fix_vector = modify->fix[n]->vector_atom;
+	int n = fix->size_vector;
+	for (i = 0; i < n; i++)
+	  combine(one,fix->compute_vector(i));
+      } else {
+	int n = fix->size_array_rows;
+	int jm1 = j - 1;
 	for (i = 0; i < nlocal; i++)
+	  combine(one,fix->compute_array(i,jm1));
+      }
+
+    } else if (flavor[m] == PERATOM) {
+      if (j == 0) {
+	double *fix_vector = fix->vector_atom;
+	int n = nlocal;
+	for (i = 0; i < n; i++)
 	  if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
 	    combine(one,fix_vector[i]);
       } else {
-      double **fix_array = modify->fix[n]->array_atom;
+	double **fix_array = fix->array_atom;
+	int n = nlocal;
 	int jm1 = j - 1;
 	for (i = 0; i < nlocal; i++)
 	  if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
 	    combine(one,fix_array[i][jm1]);
       }
 
+    } else if (flavor[m] == LOCAL) {
+      if (j == 0) {
+	double *fix_vector = fix->vector_local;
+	int n = fix->size_local_rows;
+	for (i = 0; i < n; i++)
+	  combine(one,fix_vector[i]);
+      } else {
+	double **fix_array = fix->array_local;
+	int n = fix->size_local_rows;
+	int jm1 = j - 1;
+	for (i = 0; i < n; i++)
+	  combine(one,fix_array[i][jm1]);
+      }
+    }
+    
   // evaluate atom-style variable
 
   } else if (which[m] == VARIABLE) {
@@ -136,3 +219,47 @@ double ComputeReduceRegion::compute_one(int m)
 
   return one;
 }
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeReduceRegion::count(int m)
+{
+  int n = value2index[m];
+  int j = argindex[m];
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (which[m] == X || which[m] == V || which[m] == F)
+    return group->count(igroup,iregion);
+  else if (which[m] == COMPUTE) {
+    Compute *compute = modify->compute[n];
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) return(1.0*compute->size_vector);
+      else return(1.0*compute->size_array_rows);
+    } else if (flavor[m] == PERATOM) {
+      return group->count(igroup,iregion);
+    } else if (flavor[m] == LOCAL) {
+      double ncount = compute->size_local_rows;
+      double ncountall;
+      MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+      return ncountall;
+    }
+  } else if (which[m] == FIX) {
+    Fix *fix = modify->fix[n];
+    if (flavor[m] == GLOBAL) {
+      if (j == 0) return(1.0*fix->size_vector);
+      else return(1.0*fix->size_array_rows);
+    } else if (flavor[m] == PERATOM) {
+      return group->count(igroup,iregion);
+    } else if (flavor[m] == LOCAL) {
+      double ncount = fix->size_local_rows;
+      double ncountall;
+      MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+      return ncountall;
+    }
+  } else if (which[m] == VARIABLE)
+    return group->count(igroup,iregion);
+
+  return 0.0;
+}

src/compute_reduce_region.h

@@ -25,6 +25,7 @@ class ComputeReduceRegion : public ComputeReduce {
 
  private:
   double compute_one(int);
+  double count(int);
 };
 
 }

src/compute_temp_sphere.cpp

@@ -26,8 +26,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
 #define INERTIA 0.4          // moment of inertia for sphere
 
 /* ---------------------------------------------------------------------- */

src/dump_custom.cpp

@@ -40,7 +40,8 @@ enum{ID,MOL,TYPE,MASS,
        COMPUTE,FIX,VARIABLE};
 enum{LT,LE,GT,GE,EQ,NEQ};
 enum{INT,DOUBLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
 
 /* ---------------------------------------------------------------------- */
 
@@ -960,7 +961,7 @@ void DumpCustom::parse_fields(int narg, char **arg)
       vtype[i] = DOUBLE;
 
     // compute value = c_ID
-    // if no trailing [], then arg is set to 0, else arg is between []
+    // if no trailing [], then arg is set to 0, else arg is int between []
 
     } else if (strncmp(arg[iarg],"c_",2) == 0) {
       pack_choice[i] = &DumpCustom::pack_compute;
@@ -981,14 +982,14 @@ void DumpCustom::parse_fields(int narg, char **arg)
       n = modify->find_compute(suffix);
       if (n < 0) error->all("Could not find dump custom compute ID");
       if (modify->compute[n]->peratom_flag == 0)
-	error->all("Dump custom compute ID does not compute per-atom info");
+	error->all("Dump custom compute does not compute per-atom info");
       if (argindex[i] == 0 && modify->compute[n]->size_peratom_cols > 0)
-	error->all("Dump custom compute ID does not compute per-atom vector");
+	error->all("Dump custom compute does not calculate per-atom vector");
       if (argindex[i] > 0 && modify->compute[n]->size_peratom_cols == 0)
-	error->all("Dump custom compute ID does not compute per-atom array");
+	error->all("Dump custom compute does not calculate per-atom array");
       if (argindex[i] > 0 && 
 	  argindex[i] > modify->compute[n]->size_peratom_cols)
-	error->all("Dump custom compute ID vector is not large enough");
+	error->all("Dump custom compute vector is accessed out-of-range");
 
       field2index[i] = add_compute(suffix);
       delete [] suffix;
@@ -1015,14 +1016,14 @@ void DumpCustom::parse_fields(int narg, char **arg)
       n = modify->find_fix(suffix);
       if (n < 0) error->all("Could not find dump custom fix ID");
       if (modify->fix[n]->peratom_flag == 0)
-	error->all("Dump custom fix ID does not compute per-atom info");
+	error->all("Dump custom fix does not compute per-atom info");
       if (argindex[i] == 0 && modify->fix[n]->size_peratom_cols > 0)
-	error->all("Dump custom fix ID does not compute per-atom vector");
+	error->all("Dump custom fix does not compute per-atom vector");
       if (argindex[i] > 0 && modify->fix[n]->size_peratom_cols == 0)
-	error->all("Dump custom fix ID does not compute per-atom array");
+	error->all("Dump custom fix does not compute per-atom array");
       if (argindex[i] > 0 && 
 	  argindex[i] > modify->fix[n]->size_peratom_cols)
-	error->all("Dump custom fix ID vector is not large enough");
+	error->all("Dump custom fix vector is accessed out-of-range");
 
       field2index[i] = add_fix(suffix);
       delete [] suffix;
@@ -1398,7 +1399,6 @@ void DumpCustom::pack_compute(int n)
   double *vector = compute[field2index[n]]->vector_atom;
   double **array = compute[field2index[n]]->array_atom;
   int index = argindex[n];
-
   int nlocal = atom->nlocal;
 
   if (index == 0) {
@@ -1424,7 +1424,6 @@ void DumpCustom::pack_fix(int n)
   double *vector = fix[field2index[n]]->vector_atom;
   double **array = fix[field2index[n]]->array_atom;
   int index = argindex[n];
-
   int nlocal = atom->nlocal;
 
   if (index == 0) {
@@ -1448,7 +1447,6 @@ void DumpCustom::pack_fix(int n)
 void DumpCustom::pack_variable(int n)
 {
   double *vector = vbuf[field2index[n]];
-
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++)

src/fix.h

@@ -136,7 +136,7 @@ class Fix : protected Pointers {
 
   virtual double compute_scalar() {return 0.0;}
   virtual double compute_vector(int) {return 0.0;}
-  virtual double compute_array(int) {return 0.0;}
+  virtual double compute_array(int,int) {return 0.0;}
 
   virtual int dof(int) {return 0;}
   virtual void deform(int) {}

src/fix_ave_atom.cpp

@@ -27,7 +27,8 @@
 using namespace LAMMPS_NS;
 
 enum{X,XU,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
 
 /* ---------------------------------------------------------------------- */
 
@@ -125,6 +126,7 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
   }
 
   // setup and error check
+  // for fix inputs, check that fix frequency is acceptable
 
   if (nevery <= 0) error->all("Illegal fix ave/atom command");
   if (peratom_freq < nevery || peratom_freq % nevery ||
@@ -147,19 +149,23 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
 		   "calculate a per-atom array");
       if (argindex[i] && 
 	  argindex[i] > modify->compute[icompute]->size_peratom_cols)
-	error->all("Fix ave/atom compute vector is accessed out-of-range");
+	error->all("Fix ave/atom compute array is accessed out-of-range");
+
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
 	error->all("Fix ID for fix ave/atom does not exist");
       if (modify->fix[ifix]->peratom_flag == 0)
 	error->all("Fix ave/atom fix does not calculate per-atom values");
-      if (argindex[i] && modify->fix[ifix]->size_peratom_cols != 0)
+      if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
 	error->all("Fix ave/atom fix does not calculate a per-atom vector");
       if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
 	error->all("Fix ave/atom fix does not calculate a per-atom array");
       if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
-	error->all("Fix ave/atom fix vector is accessed out-of-range");
+	error->all("Fix ave/atom fix array is accessed out-of-range");
+      if (nevery % modify->fix[ifix]->peratom_freq)
+	error->all("Fix for fix ave/atom not computed at compatible time");
+
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
@@ -240,7 +246,6 @@ int FixAveAtom::setmask()
 void FixAveAtom::init()
 {
   // set indices and check validity of all computes,fixes,variables
-  // check that fix frequency is acceptable
 
   for (int m = 0; m < nvalues; m++) {
     if (which[m] == COMPUTE) {
@@ -255,9 +260,6 @@ void FixAveAtom::init()
 	error->all("Fix ID for fix ave/atom does not exist");
       value2index[m] = ifix;
 
-      if (nevery % modify->fix[ifix]->peratom_freq)
-	error->all("Fix for fix ave/atom not computed at compatible time");
-
     } else if (which[m] == VARIABLE) {
       int ivariable = input->variable->find(ids[m]);
       if (ivariable < 0) 

src/fix_ave_spatial.cpp

@@ -36,7 +36,8 @@ enum{X,V,F,DENSITY_NUMBER,DENSITY_MASS,COMPUTE,FIX,VARIABLE};
 enum{SAMPLE,ALL};
 enum{BOX,LATTICE,REDUCED};
 enum{ONE,RUNNING,WINDOW};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
 
 #define BIG 1000000000
 
@@ -154,6 +155,9 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   scaleflag = LATTICE;
   fp = NULL;
   ave = ONE;
+  char *title1 = NULL;
+  char *title2 = NULL;
+  char *title3 = NULL;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"norm") == 0) {
@@ -193,6 +197,27 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
       }
       iarg += 2;
       if (ave == WINDOW) iarg++;
+    } else if (strcmp(arg[iarg],"title1") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+      delete [] title1;
+      int n = strlen(arg[iarg+1]) + 1;
+      title1 = new char[n];
+      strcpy(title1,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"title2") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+      delete [] title2;
+      int n = strlen(arg[iarg+1]) + 1;
+      title2 = new char[n];
+      strcpy(title2,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"title3") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+      delete [] title3;
+      int n = strlen(arg[iarg+1]) + 1;
+      title3 = new char[n];
+      strcpy(title3,arg[iarg+1]);
+      iarg += 2;
     } else error->all("Illegal fix ave/spatial command");
   }
 
@@ -243,27 +268,38 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
-  // print header into file
+  // print file comment lines
 
   if (fp && me == 0) {
-    fprintf(fp,"# Spatial-averaged data for fix %s and group %s\n",id,arg[1]);
+    if (title1) fprintf(fp,"%s\n",title1);
-    fprintf(fp,"# TimeStep Number-of-layers\n");
+    else fprintf(fp,"# Spatial-averaged data for fix %s and group %s\n",
-    fprintf(fp,"# Layer Coordinate Natoms");
+		 id,arg[1]);
-    for (int i = 0; i < nvalues; i++)
+    if (title2) fprintf(fp,"%s\n",title2);
+    else fprintf(fp,"# Timestep Number-of-layers\n");
+    if (title3) fprintf(fp,"%s\n",title3);
+    else {
+      fprintf(fp,"# Layer Coord Ncount");
+      for (int i = 0; i < nvalues; i++) {
 	if (which[i] == COMPUTE) fprintf(fp," c_%s",ids[i]);
 	else if (which[i] == FIX) fprintf(fp," f_%s",ids[i]);
 	else if (which[i] == VARIABLE) fprintf(fp," v_%s",ids[i]);
 	else fprintf(fp," %s",arg[9+i]);
+      }
       fprintf(fp,"\n");
     }
+  }
 
-  // this fix produces a global vector
+  delete [] title1;
-  // set size_vector to BIG since compute_vector() checks bounds on-the-fly
+  delete [] title2;
+  delete [] title3;
 
-  vector_flag = 1;
+  // this fix produces a global array
-  size_vector = BIG;
+
+  array_flag = 1;
+  size_local_rows = BIG;
+  size_local_cols = nvalues+2;
   global_freq = nfreq;
-  extvector = 0;
+  extarray = 0;
 
   // setup scaling
 
@@ -803,18 +839,18 @@ void FixAveSpatial::end_of_step()
 }
 
 /* ----------------------------------------------------------------------
-   return Nth vector value
+   return I,J array value
-   since values_sum is 2d array, map N into ilayer and ivalue
+   if I exceeds current layers, return 0.0 instead of generating an error
-   if ilayer exceeds current layers, return 0.0 instead of generate an error
+   1st column = bin coord, 2nd column = count, remaining columns = Nvalues
 ------------------------------------------------------------------------- */
 
-double FixAveSpatial::compute_vector(int n)
+double FixAveSpatial::compute_array(int i, int j)
 {
-  int ivalue = n % nvalues;
-  int ilayer = n / nvalues;
-  if (ilayer >= nlayers) return 0.0;
   if (values_total == NULL) return 0.0;
-  return values_total[ilayer][ivalue]/norm;
+  if (i >= nlayers) return 0.0;
+  if (j == 0) return coord[i];
+  if (j == 1) return count_total[i]/norm;
+  return values_total[i][j]/norm;
 }
 
 /* ----------------------------------------------------------------------

src/fix_ave_spatial.h

@@ -27,7 +27,7 @@ class FixAveSpatial : public Fix {
   void init();
   void setup(int);
   void end_of_step();
-  double compute_vector(int);
+  double compute_array(int,int);
   double memory_usage();
 
  private:
@@ -35,6 +35,7 @@ class FixAveSpatial : public Fix {
   int nrepeat,nfreq,nvalid,irepeat;
   int dim,originflag,normflag;
   double origin,delta;
+  char *tstring,*sstring;
   int *which,*argindex,*value2index;
   char **ids;
   FILE *fp;

src/fix_ave_time.cpp

@@ -31,7 +31,10 @@ using namespace LAMMPS_NS;
 
 enum{COMPUTE,FIX,VARIABLE};
 enum{ONE,RUNNING,WINDOW};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
 
 /* ---------------------------------------------------------------------- */
 
@@ -129,6 +132,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   }
 
   // setup and error check
+  // for fix inputs, check that fix frequency is acceptable
 
   if (nevery <= 0) error->all("Illegal fix ave/time command");
   if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
@@ -145,6 +149,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
 	error->all("Fix ave/time compute does not calculate a vector");
       if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
 	error->all("Fix ave/time compute vector is accessed out-of-range");
+
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
@@ -155,6 +160,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
 	error->all("Fix ave/time fix does not calculate a vector");
       if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
 	error->all("Fix ave/time fix vector is accessed out-of-range");
+      if (nevery % modify->fix[ifix]->global_freq)
+	error->all("Fix for fix ave/time not computed at compatible time");
+
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0)
@@ -285,8 +293,7 @@ int FixAveTime::setmask()
 
 void FixAveTime::init()
 {
-  // set indices and check validity of all computes,fixes,variables
+  // set current indices for all computes,fixes,variables
-  // check that fix frequency is acceptable
 
   for (int i = 0; i < nvalues; i++) {
     if (which[i] == COMPUTE) {
@@ -301,9 +308,6 @@ void FixAveTime::init()
 	error->all("Fix ID for fix ave/time does not exist");
       value2index[i] = ifix;
 
-      if (nevery % modify->fix[ifix]->global_freq)
-	error->all("Fix for fix ave/time not computed at compatible time");
-
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
       if (ivariable < 0) 
@@ -440,7 +444,6 @@ void FixAveTime::end_of_step()
 
 /* ----------------------------------------------------------------------
    return scalar value
-   could be ONE, RUNNING, or WINDOW value
 ------------------------------------------------------------------------- */
 
 double FixAveTime::compute_scalar()
@@ -450,12 +453,21 @@ double FixAveTime::compute_scalar()
 }
 
 /* ----------------------------------------------------------------------
-   return Nth vector value
+   return Ith vector value
-   could be ONE, RUNNING, or WINDOW value
 ------------------------------------------------------------------------- */
 
-double FixAveTime::compute_vector(int n)
+double FixAveTime::compute_vector(int i)
 {
-  if (norm) return vector_total[n]/norm;
+  if (norm) return vector_total[i]/norm;
+  return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   return I,J array value
+------------------------------------------------------------------------- */
+
+double FixAveTime::compute_array(int i, int j)
+{
+  if (norm) return array_total[i][j]/norm;
   return 0.0;
 }

src/fix_ave_time.h

@@ -29,6 +29,7 @@ class FixAveTime : public Fix {
   void end_of_step();
   double compute_scalar();
   double compute_vector(int);
+  double compute_array(int,int);
 
  private:
   int me,nvalues;
@@ -43,6 +44,7 @@ class FixAveTime : public Fix {
   int norm,iwindow,window_limit;
   double *vector_total;
   double **vector_list;
+  double **array_total;
 };
 
 }

src/style.h

@@ -90,6 +90,7 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_pe.h"
 #include "compute_pe_atom.h"
 #include "compute_pressure.h"
+#include "compute_rdf.h"
 #include "compute_reduce.h"
 #include "compute_reduce_region.h"
 #include "compute_erotate_sphere.h"
@@ -121,6 +122,7 @@ ComputeStyle(pe/atom,ComputePEAtom)
 ComputeStyle(pressure,ComputePressure)
 ComputeStyle(reduce,ComputeReduce)
 ComputeStyle(reduce/region,ComputeReduceRegion)
+ComputeStyle(rdf,ComputeRDF)
 ComputeStyle(erotate/sphere,ComputeERotateSphere)
 ComputeStyle(stress/atom,ComputeStressAtom)
 ComputeStyle(temp,ComputeTemp)

src/thermo.cpp

@@ -53,7 +53,11 @@ using namespace LAMMPS_NS;
 enum{IGNORE,WARN,ERROR};           // same as write_restart.cpp
 enum{ONELINE,MULTILINE};
 enum{INT,FLOAT};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{SCALAR,VECTOR,ARRAY};
+
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
 
 #define MAXLINE 1024
 #define DELTA 8
@@ -271,16 +275,21 @@ void Thermo::compute(int flag)
   // which = 0 is global scalar, which = 1 is global vector
 
   for (i = 0; i < ncompute; i++)
-    if (compute_which[i] == 0) {
+    if (compute_which[i] == SCALAR) {
       if (!(computes[i]->invoked_flag & INVOKED_SCALAR)) {
 	computes[i]->compute_scalar();
 	computes[i]->invoked_flag |= INVOKED_SCALAR;
       }
-    } else {
+    } else if (compute_which[i] == VECTOR) {
       if (!(computes[i]->invoked_flag & INVOKED_VECTOR)) {
 	computes[i]->compute_vector();
 	computes[i]->invoked_flag |= INVOKED_VECTOR;
       }
+    } else if (compute_which[i] == ARRAY) {
+      if (!(computes[i]->invoked_flag & INVOKED_ARRAY)) {
+	computes[i]->compute_array();
+	computes[i]->invoked_flag |= INVOKED_ARRAY;
+      }
     }
 
   // if lineflag = MULTILINE, prepend step/cpu header line
@@ -502,7 +511,8 @@ void Thermo::allocate()
   for (int i = 0; i < n; i++) format_user[i] = NULL;
 
   field2index = new int[n];
-  argindex = new int[n];
+  argindex1 = new int[n];
+  argindex2 = new int[n];
 
   // factor of 3 is max number of computes a single field can add
 
@@ -539,7 +549,8 @@ void Thermo::deallocate()
   delete [] format_user;
 
   delete [] field2index;
-  delete [] argindex;
+  delete [] argindex1;
+  delete [] argindex2;
 
   for (int i = 0; i < ncompute; i++) delete [] id_compute[i];
   delete [] id_compute;
@@ -580,56 +591,56 @@ void Thermo::parse_fields(char *str)
 
     } else if (strcmp(word,"temp") == 0) {
       addfield("Temp",&Thermo::compute_temp,FLOAT);
-      index_temp = add_compute(id_temp,0);
+      index_temp = add_compute(id_temp,SCALAR);
     } else if (strcmp(word,"press") == 0) {
       addfield("Press",&Thermo::compute_press,FLOAT);
-      index_press_scalar = add_compute(id_press,0);
+      index_press_scalar = add_compute(id_press,SCALAR);
     } else if (strcmp(word,"pe") == 0) {
       addfield("PotEng",&Thermo::compute_pe,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"ke") == 0) {
       addfield("KinEng",&Thermo::compute_ke,FLOAT);
-      index_temp = add_compute(id_temp,0);
+      index_temp = add_compute(id_temp,SCALAR);
     } else if (strcmp(word,"etotal") == 0) {
       addfield("TotEng",&Thermo::compute_etotal,FLOAT);
-      index_temp = add_compute(id_temp,0);
+      index_temp = add_compute(id_temp,SCALAR);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"enthalpy") == 0) {
       addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
-      index_temp = add_compute(id_temp,0);
+      index_temp = add_compute(id_temp,SCALAR);
-      index_press_scalar = add_compute(id_press,0);
+      index_press_scalar = add_compute(id_press,SCALAR);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
 
     } else if (strcmp(word,"evdwl") == 0) {
       addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"ecoul") == 0) {
       addfield("E_coul",&Thermo::compute_ecoul,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"epair") == 0) {
       addfield("E_pair",&Thermo::compute_epair,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"ebond") == 0) {
       addfield("E_bond",&Thermo::compute_ebond,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"eangle") == 0) {
       addfield("E_angle",&Thermo::compute_eangle,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"edihed") == 0) {
       addfield("E_dihed",&Thermo::compute_edihed,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"eimp") == 0) {
       addfield("E_impro",&Thermo::compute_eimp,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"emol") == 0) {
       addfield("E_mol",&Thermo::compute_emol,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"elong") == 0) {
       addfield("E_long",&Thermo::compute_elong,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
     } else if (strcmp(word,"etail") == 0) {
       addfield("E_tail",&Thermo::compute_etail,FLOAT);
-      index_pe = add_compute(id_pe,0);
+      index_pe = add_compute(id_pe,SCALAR);
 
     } else if (strcmp(word,"vol") == 0) {
       addfield("Volume",&Thermo::compute_vol,FLOAT);
@@ -662,25 +673,25 @@ void Thermo::parse_fields(char *str)
 
     } else if (strcmp(word,"pxx") == 0) {
       addfield("Pxx",&Thermo::compute_pxx,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
     } else if (strcmp(word,"pyy") == 0) {
       addfield("Pyy",&Thermo::compute_pyy,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
     } else if (strcmp(word,"pzz") == 0) {
       addfield("Pzz",&Thermo::compute_pzz,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
     } else if (strcmp(word,"pxy") == 0) {
       addfield("Pxy",&Thermo::compute_pxy,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
     } else if (strcmp(word,"pxz") == 0) {
       addfield("Pxz",&Thermo::compute_pxz,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
     } else if (strcmp(word,"pyz") == 0) {
       addfield("Pyz",&Thermo::compute_pyz,FLOAT);
-      index_press_vector = add_compute(id_press,1);
+      index_press_vector = add_compute(id_press,VECTOR);
 
     // compute value = c_ID, fix value = f_ID, variable value = v_ID
-    // if no trailing [], then arg is set to 0, else arg is between []
+    // count trailing [] and store int arguments
     // copy = at most 8 chars of ID to pass to addfield
 
     } else if ((strncmp(word,"c_",2) == 0) || (strncmp(word,"f_",2) == 0) ||
@@ -693,38 +704,66 @@ void Thermo::parse_fields(char *str)
       strncpy(copy,id,8);
       copy[8] = '\0';
 
+      // parse zero or one or two trailing brackets from ID
+      // argindex1,argindex2 = int inside each bracket pair, 0 if no bracket
+
       char *ptr = strchr(id,'[');
-      if (ptr) {
+      if (ptr == NULL) argindex1[nfield] = 0;
-	if (id[strlen(id)-1] != ']')
+      else {
-	  error->all("Invalid keyword in thermo_style custom command");
-	argindex[nfield] = atoi(ptr+1);
 	*ptr = '\0';
-      } else argindex[nfield] = 0;
+	argindex1[nfield] = input->variable->int_between_brackets(ptr);
+	ptr++;
+	if (*ptr == '[') {
+	  argindex2[nfield] = input->variable->int_between_brackets(ptr);
+	  ptr++;
+	} else argindex2[nfield] = 0;
+      }
 
       if (word[0] == 'c') {
 	n = modify->find_compute(id);
 	if (n < 0) error->all("Could not find thermo custom compute ID");
-	if (argindex[nfield] == 0 && modify->compute[n]->scalar_flag == 0)
+	if (argindex1[nfield] == 0 && modify->compute[n]->scalar_flag == 0)
-	  error->all("Thermo compute ID does not compute scalar info");
+	  error->all("Thermo compute does not compute scalar");
-	if (argindex[nfield] > 0 && modify->compute[n]->vector_flag == 0)
+	if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
-	  error->all("Thermo compute ID does not compute vector info");
+	  if (modify->compute[n]->vector_flag == 0)
-	if (argindex[nfield] > 0 && 
+	    error->all("Thermo compute does not compute vector");
-	    argindex[nfield] > modify->compute[n]->size_vector)
+	  if (argindex1[nfield] > modify->compute[n]->size_vector)
-	  error->all("Thermo compute ID vector is not large enough");
+	    error->all("Thermo compute vector is accessed out-of-range");
+	}
+	if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
+	  if (modify->compute[n]->array_flag == 0)
+	    error->all("Thermo compute does not compute array");
+	  if (argindex1[nfield] > modify->compute[n]->size_array_rows ||
+	      argindex2[nfield] > modify->compute[n]->size_array_cols)
+	    error->all("Thermo compute array is accessed out-of-range");
+	}
 
-	field2index[nfield] = add_compute(id,MIN(argindex[nfield],1));
+	if (argindex1[nfield] == 0)
+	  field2index[nfield] = add_compute(id,SCALAR);
+	else if (argindex2[nfield] == 0)
+	  field2index[nfield] = add_compute(id,VECTOR);
+	else 
+	  field2index[nfield] = add_compute(id,ARRAY);
 	addfield(copy,&Thermo::compute_compute,FLOAT);
 
       } else if (word[0] == 'f') {
 	n = modify->find_fix(id);
 	if (n < 0) error->all("Could not find thermo custom fix ID");
-	if (argindex[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
+	if (argindex1[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
-	  error->all("Thermo fix ID does not compute scalar info");
+	  error->all("Thermo fix does not compute scalar");
-	if (argindex[nfield] > 0 && modify->fix[n]->vector_flag == 0)
+	if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
-	  error->all("Thermo fix ID does not compute vector info");
+	  if (modify->fix[n]->vector_flag == 0)
-	if (argindex[nfield] > 0 && 
+	    error->all("Thermo fix does not compute vector");
-	    argindex[nfield] > modify->fix[n]->size_vector)
+	  if (argindex1[nfield] > modify->fix[n]->size_vector)
-	  error->all("Thermo fix ID vector is not large enough");
+	    error->all("Thermo fix vector is accessed out-of-range");
+	}
+	if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
+	  if (modify->fix[n]->array_flag == 0)
+	    error->all("Thermo fix does not compute array");
+	  if (argindex1[nfield] > modify->fix[n]->size_array_rows ||
+	      argindex2[nfield] > modify->fix[n]->size_array_cols)
+	    error->all("Thermo fix array is accessed out-of-range");
+	}
 
 	field2index[nfield] = add_fix(id);
 	addfield(copy,&Thermo::compute_fix,FLOAT);
@@ -734,6 +773,8 @@ void Thermo::parse_fields(char *str)
 	if (n < 0) error->all("Could not find thermo custom variable name");
 	if (input->variable->equalstyle(n) == 0)
 	  error->all("Thermo custom variable is not equal-style variable");
+	if (argindex1[nfield]) 
+	  error->all("Thermo custom variable cannot be indexed");
 
 	field2index[nfield] = add_variable(id);
 	addfield(copy,&Thermo::compute_variable,FLOAT);
@@ -1150,17 +1191,22 @@ int Thermo::evaluate_keyword(char *word, double *answer)
 
 void Thermo::compute_compute()
 {
-  Compute *compute = computes[field2index[ifield]];
+  int m = field2index[ifield];
-  if (argindex[ifield] == 0) {
+  Compute *compute = computes[m];
+
+  if (compute_which[m] == SCALAR) {
     dvalue = compute->scalar;
     if (normflag && compute->extscalar) dvalue /= natoms;
-  } else {
+  } else if (compute_which[m] == VECTOR) {
-    dvalue = compute->vector[argindex[ifield]-1];
+    dvalue = compute->vector[argindex1[ifield]-1];
     if (normflag) {
       if (compute->extvector == 0) return;
       else if (compute->extvector == 1) dvalue /= natoms;
-      else if (compute->extlist[argindex[ifield]-1]) dvalue /= natoms;
+      else if (compute->extlist[argindex1[ifield]-1]) dvalue /= natoms;
     }
+  } else {
+    dvalue = compute->array[argindex1[ifield]-1][argindex2[ifield]-1];
+    if (normflag && compute->extarray) dvalue /= natoms;
   }
 }
 
@@ -1168,17 +1214,22 @@ void Thermo::compute_compute()
 
 void Thermo::compute_fix()
 {
-  Fix *fix = fixes[field2index[ifield]];
+  int m = field2index[ifield];
-  if (argindex[ifield] == 0) {
+  Fix *fix = fixes[m];
+
+  if (argindex1[ifield] == 0) {
     dvalue = fix->compute_scalar();
     if (normflag && fix->extscalar) dvalue /= natoms;
-  } else {
+  } else if (argindex2[ifield] == 0) {
-    dvalue = fix->compute_vector(argindex[ifield]-1);
+    dvalue = fix->compute_vector(argindex1[ifield]-1);
     if (normflag) {
       if (fix->extvector == 0) return;
       else if (fix->extvector == 1) dvalue /= natoms;
-      else if (fix->extlist[argindex[ifield]-1]) dvalue /= natoms;
+      else if (fix->extlist[argindex1[ifield]-1]) dvalue /= natoms;
     }
+  } else {
+    dvalue = fix->compute_array(argindex1[ifield]-1,argindex2[ifield]-1);
+    if (normflag && fix->extarray) dvalue /= natoms;
   }
 }
 

src/thermo.h

@@ -62,8 +62,8 @@ class Thermo : protected Pointers {
   double dvalue,natoms;  // dvalue = double value to print
   int ifield;            // which field in thermo output is being computed
   int *field2index;      // which compute,fix,variable calcs this field
-  int *argindex;         // index into compute,fix scalar,vector
+  int *argindex1;        // indices into compute,fix scalar,vector
-
+  int *argindex2;
                          // data for keyword-specific Compute objects
                          // index = where they are in computes list
                          // id = ID of Compute objects

src/variable.cpp

@@ -37,13 +37,17 @@ using namespace LAMMPS_NS;
 
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
 enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP,ATOM};
 enum{ARG,OP};
 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
        SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,
        CEIL,FLOOR,ROUND,VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+
 /* ---------------------------------------------------------------------- */
 
 Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
@@ -524,11 +528,11 @@ void Variable::copy(int narg, char **from, char **to)
                       sqrt(x),exp(x),ln(x),log(x),
 		      sin(x),cos(x),tan(x),asin(x),acos(x),atan(x)
      group function = count(group), mass(group), xcm(group,x), ...
-     atom value = x[123], y[3], vx[34], ...
+     atom value = x[i], y[i], vx[i], ...
-     atom vector = x[], y[], vx[], ...
+     atom vector = x, y, vx, ...
-     compute = c_ID, c_ID[2], c_ID[], c_ID[][2], c_ID[N], c_ID[N][2]
+     compute = c_ID, c_ID[i], c_ID[i][j]
-     fix = f_ID, f_ID[2], f_ID[], f_ID[][2], f_ID[N], f_ID[N][2]
+     fix = f_ID, f_ID[i], f_ID[i][j]
-     variable = v_name, v_name[], v_name[N]
+     variable = v_name, v_name[i]
    if tree ptr passed in, return ptr to parse tree created for formula
    if no tree ptr passed in, return value of string as double
 ------------------------------------------------------------------------- */
@@ -538,6 +542,7 @@ double Variable::evaluate(char *str, Tree **tree)
   int op,opprevious;
   double value1,value2;
   char onechar;
+  char *ptr;
 
   double argstack[MAXLEVEL];
   Tree *treestack[MAXLEVEL];
@@ -613,7 +618,8 @@ double Variable::evaluate(char *str, Tree **tree)
       delete [] number;
 
     // ----------------
-    // letter: c_ID, f_ID, v_name, exp(), xcm(,), x[234], x[], vol
+    // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
+    //         v_name, v_name[], exp(), xcm(,), x, x[], vol
     // ----------------
 
     } else if (islower(onechar)) {
@@ -637,6 +643,9 @@ double Variable::evaluate(char *str, Tree **tree)
       // ----------------
 
       if (strncmp(word,"c_",2) == 0) {
+	if (domain->box_exist == 0)
+	  error->all("Variable evaluation before simulation box is defined");
+ 
 	n = strlen(word) - 2 + 1;
 	char *id = new char[n];
 	strcpy(id,&word[2]);
@@ -646,28 +655,27 @@ double Variable::evaluate(char *str, Tree **tree)
 	Compute *compute = modify->compute[icompute];
 	delete [] id;
 
-	if (domain->box_exist == 0)
-	  error->all("Variable evaluation before simulation box is defined");
- 
 	// parse zero or one or two trailing brackets
 	// point i beyond last bracket
 	// nbracket = # of bracket pairs
-	// index1,index2 = int inside each bracket pair, 0 if empty
+	// index1,index2 = int inside each bracket pair
 
 	int nbracket,index1,index2;
 	if (str[i] != '[') nbracket = 0;
 	else {
 	  nbracket = 1;
-	  i = int_between_brackets(str,i,index1,1);
+	  ptr = &str[i];
-	  i++;
+	  index1 = int_between_brackets(ptr);
+	  i = ptr-str+1;
 	  if (str[i] == '[') {
 	    nbracket = 2;
-	    i = int_between_brackets(str,i,index2,0);
+	    ptr = &str[i];
-	    i++;
+	    index2 = int_between_brackets(ptr);
+	    i = ptr-str+1;
 	  }
 	}
 
-        // c_ID = global scalar
+        // c_ID = scalar from global scalar
 
 	if (nbracket == 0 && compute->scalar_flag) {
 
@@ -689,12 +697,13 @@ double Variable::evaluate(char *str, Tree **tree)
 	    treestack[ntreestack++] = newtree;
 	  } else argstack[nargstack++] = value1;
 
-        // c_ID[2] = global vector
+        // c_ID[i] = scalar from global vector
 
-	} else if (nbracket == 1 && index1 > 0 && compute->vector_flag) {
+	} else if (nbracket == 1 && compute->vector_flag) {
 
 	  if (index1 > compute->size_vector)
-	      error->all("Compute vector in variable formula is too small");
+	    error->all("Variable formula compute vector "
+		       "is accessed out-of-range");
 	  if (update->whichflag == 0) {
 	    if (compute->invoked_vector != update->ntimestep)
 	      error->all("Compute used in variable between runs "
@@ -713,10 +722,76 @@ double Variable::evaluate(char *str, Tree **tree)
 	    treestack[ntreestack++] = newtree;
 	  } else argstack[nargstack++] = value1;
 
-        // c_ID[] = per-atom vector
+        // c_ID[i][j] = scalar from global array
 
-	} else if (nbracket == 1 && index1 == 0 && 
+	} else if (nbracket == 2 && compute->array_flag) {
-		   compute->peratom_flag && compute->size_peratom_cols == 0) {
+
+	  if (index1 > compute->size_array_rows)
+	    error->all("Variable formula compute vector "
+		       "is accessed out-of-range");
+	  if (index2 > compute->size_array_cols)
+	    error->all("Variable formula compute vector "
+		       "is accessed out-of-range");
+	  if (update->whichflag == 0) {
+	    if (compute->invoked_array != update->ntimestep)
+	      error->all("Compute used in variable between runs "
+			 "is not current");
+	  } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+	    compute->compute_array();
+	    compute->invoked_flag |= INVOKED_ARRAY;
+	  }
+
+	  value1 = compute->array[index1-1][index2-1];
+	  if (tree) {
+	    Tree *newtree = new Tree();
+	    newtree->type = VALUE;
+	    newtree->value = value1;
+	    newtree->left = newtree->right = NULL;
+	    treestack[ntreestack++] = newtree;
+	  } else argstack[nargstack++] = value1;
+
+        // c_ID[i] = scalar from per-atom vector
+
+	} else if (nbracket == 1 && compute->peratom_flag && 
+		   compute->size_peratom_cols == 0) {
+
+	  if (update->whichflag == 0) {
+	    if (compute->invoked_peratom != update->ntimestep)
+	      error->all("Compute used in variable between runs "
+			 "is not current");
+	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+	    compute->compute_peratom();
+	    compute->invoked_flag |= INVOKED_PERATOM;
+	  }
+
+	  peratom2global(1,NULL,compute->vector_atom,1,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // c_ID[i][j] = scalar from per-atom array
+
+	} else if (nbracket == 2 && compute->peratom_flag &&
+		   compute->size_peratom_cols > 0) {
+
+	  if (index2 > compute->size_peratom_cols)
+	    error->all("Variable formula compute vector "
+		       "is accessed out-of-range");
+	  if (update->whichflag == 0) {
+	    if (compute->invoked_peratom != update->ntimestep)
+	      error->all("Compute used in variable between runs "
+			 "is not current");
+	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+	    compute->compute_peratom();
+	    compute->invoked_flag |= INVOKED_PERATOM;
+	  }
+
+	  peratom2global(1,NULL,&compute->array_atom[0][index2-1],
+			 compute->size_peratom_cols,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // c_ID = vector from per-atom vector
+
+	} else if (nbracket == 0 && compute->peratom_flag && 
+		   compute->size_peratom_cols == 0) {
 
 	  if (tree == NULL)
 	    error->all("Per-atom compute in equal-style variable formula");
@@ -736,32 +811,16 @@ double Variable::evaluate(char *str, Tree **tree)
 	  newtree->left = newtree->right = NULL;
 	  treestack[ntreestack++] = newtree;
 
-        // c_ID[N] = global value from per-atom vector
+        // c_ID[i] = vector from per-atom array
 
-	} else if (nbracket == 1 && index1 > 0 && 
+	} else if (nbracket == 1 && compute->peratom_flag &&
-		   compute->peratom_flag && compute->size_peratom_cols == 0) {
+		   compute->size_peratom_cols > 0) {
-
-	  if (update->whichflag == 0) {
-	    if (compute->invoked_peratom != update->ntimestep)
-	      error->all("Compute used in variable between runs "
-			 "is not current");
-	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
-	    compute->compute_peratom();
-	    compute->invoked_flag |= INVOKED_PERATOM;
-	  }
-
-	  peratom2global(1,NULL,compute->vector_atom,1,index1,
-			 tree,treestack,ntreestack,argstack,nargstack);
-
-        // c_ID[][2] = per-atom array
-
-	} else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
-		   compute->peratom_flag) {
 
 	  if (tree == NULL)
 	    error->all("Per-atom compute in equal-style variable formula");
-	  if (index2 > compute->size_peratom_cols)
+	  if (index1 > compute->size_peratom_cols)
-	    error->all("Compute vector in variable formula is too small");
+	    error->all("Variable formula compute vector "
+		       "is accessed out-of-range");
 	  if (update->whichflag == 0) {
 	    if (compute->invoked_peratom != update->ntimestep)
 	      error->all("Compute used in variable between runs "
@@ -773,31 +832,11 @@ double Variable::evaluate(char *str, Tree **tree)
 
 	  Tree *newtree = new Tree();
 	  newtree->type = ATOMARRAY;
-	  newtree->array = &compute->array_atom[0][index2-1];
+	  newtree->array = &compute->array_atom[0][index1-1];
 	  newtree->nstride = compute->size_peratom_cols;
 	  newtree->left = newtree->right = NULL;
 	  treestack[ntreestack++] = newtree;
 
-        // c_ID[N][2] = global value from per-atom array
-
-	} else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
-		   compute->peratom_flag) {
-
-	  if (index2 > compute->size_peratom_cols)
-	    error->all("Compute vector in variable formula is too small");
-	  if (update->whichflag == 0) {
-	    if (compute->invoked_peratom != update->ntimestep)
-	      error->all("Compute used in variable between runs "
-			 "is not current");
-	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
-	    compute->compute_peratom();
-	    compute->invoked_flag |= INVOKED_PERATOM;
-	  }
-
-	  peratom2global(1,NULL,&compute->array_atom[0][index2-1],
-			 compute->size_peratom_cols,index1,
-			 tree,treestack,ntreestack,argstack,nargstack);
-
 	} else error->all("Mismatched compute in variable formula");
 
       // ----------------
@@ -805,6 +844,9 @@ double Variable::evaluate(char *str, Tree **tree)
       // ----------------
 
       } else if (strncmp(word,"f_",2) == 0) {
+	if (domain->box_exist == 0)
+	  error->all("Variable evaluation before simulation box is defined");
+ 
 	n = strlen(word) - 2 + 1;
 	char *id = new char[n];
 	strcpy(id,&word[2]);
@@ -814,33 +856,33 @@ double Variable::evaluate(char *str, Tree **tree)
 	Fix *fix = modify->fix[ifix];
 	delete [] id;
 
-	if (domain->box_exist == 0)
-	  error->all("Variable evaluation before simulation box is defined");
- 
 	// parse zero or one or two trailing brackets
 	// point i beyond last bracket
 	// nbracket = # of bracket pairs
-	// index1,index2 = int inside each bracket pair, 0 if empty
+	// index1,index2 = int inside each bracket pair
 
 	int nbracket,index1,index2;
 	if (str[i] != '[') nbracket = 0;
 	else {
 	  nbracket = 1;
-	  i = int_between_brackets(str,i,index1,1);
+	  ptr = &str[i];
-	  i++;
+	  index1 = int_between_brackets(ptr);
+	  i = ptr-str+1;
 	  if (str[i] == '[') {
 	    nbracket = 2;
-	    i = int_between_brackets(str,i,index2,0);
+	    ptr = &str[i];
-	    i++;
+	    index2 = int_between_brackets(ptr);
+	    i = ptr-str+1;
 	  }
 	}
 
-        // f_ID = global scalar
+        // f_ID = scalar from global scalar
 
 	if (nbracket == 0 && fix->scalar_flag) {
 
 	  if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
 	    error->all("Fix in variable not computed at compatible time");
+
 	  value1 = fix->compute_scalar();
 	  if (tree) {
 	    Tree *newtree = new Tree();
@@ -850,14 +892,15 @@ double Variable::evaluate(char *str, Tree **tree)
 	    treestack[ntreestack++] = newtree;
 	  } else argstack[nargstack++] = value1;
 
-        // f_ID[2] = global vector
+        // f_ID[i] = scalar from global vector
 
-	} else if (nbracket == 1 && index1 > 0 && fix->vector_flag) {
+	} else if (nbracket == 1 && fix->vector_flag) {
 
 	  if (index1 > fix->size_vector)
-	      error->all("Fix vector in variable formula is too small");
+	    error->all("Variable formula fix vector is accessed out-of-range");
 	  if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
 	    error->all("Fix in variable not computed at compatible time");
+
 	  value1 = fix->compute_vector(index1-1);
 	  if (tree) {
 	    Tree *newtree = new Tree();
@@ -867,16 +910,64 @@ double Variable::evaluate(char *str, Tree **tree)
 	    treestack[ntreestack++] = newtree;
 	  } else argstack[nargstack++] = value1;
 
-        // f_ID[] = per-atom vector
+        // f_ID[i][j] = scalar from global array
 
-	} else if (nbracket == 1 && index1 == 0 && 
+	} else if (nbracket == 2 && fix->array_flag) {
-		   fix->peratom_flag && fix->size_peratom_cols == 0) {
+
+	  if (index1 > fix->size_array_rows)
+	    error->all("Variable formula fix vector is accessed out-of-range");
+	  if (index2 > fix->size_array_cols)
+	    error->all("Variable formula fix vector is accessed out-of-range");
+	  if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+	    error->all("Fix in variable not computed at compatible time");
+
+	  value1 = fix->compute_array(index1-1,index2-1);
+	  if (tree) {
+	    Tree *newtree = new Tree();
+	    newtree->type = VALUE;
+	    newtree->value = value1;
+	    newtree->left = newtree->right = NULL;
+	    treestack[ntreestack++] = newtree;
+	  } else argstack[nargstack++] = value1;
+
+        // f_ID[i] = scalar from per-atom vector
+
+	} else if (nbracket == 1 && fix->peratom_flag && 
+		   fix->size_peratom_cols == 0) {
+
+	  if (update->whichflag > 0 && 
+	      update->ntimestep % fix->peratom_freq)
+	    error->all("Fix in variable not computed at compatible time");
+
+	  peratom2global(1,NULL,fix->vector_atom,1,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // f_ID[i][j] = scalar from per-atom array
+
+	} else if (nbracket == 2 && fix->peratom_flag &&
+		   fix->size_peratom_cols > 0) {
+
+	  if (index2 > fix->size_peratom_cols)
+	    error->all("Variable formula fix vector is accessed out-of-range");
+	  if (update->whichflag > 0 && 
+	      update->ntimestep % fix->peratom_freq)
+	    error->all("Fix in variable not computed at compatible time");
+
+	  peratom2global(1,NULL,&fix->array_atom[0][index2-1],
+			 fix->size_peratom_cols,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // f_ID = vector from per-atom vector
+
+	} else if (nbracket == 0 && fix->peratom_flag && 
+		   fix->size_peratom_cols == 0) {
 
 	  if (tree == NULL)
 	    error->all("Per-atom fix in equal-style variable formula");
 	  if (update->whichflag > 0 && 
 	      update->ntimestep % fix->peratom_freq)
 	    error->all("Fix in variable not computed at compatible time");
+
 	  Tree *newtree = new Tree();
 	  newtree->type = ATOMARRAY;
 	  newtree->array = fix->vector_atom;
@@ -884,49 +975,26 @@ double Variable::evaluate(char *str, Tree **tree)
 	  newtree->left = newtree->right = NULL;
 	  treestack[ntreestack++] = newtree;
 
-        // f_ID[N] = global value from per-atom vector
+        // f_ID[i] = vector from per-atom array
 
-	} else if (nbracket == 1 && index1 > 0 && 
+	} else if (nbracket == 1 && fix->peratom_flag &&
-		   fix->peratom_flag && fix->size_peratom_cols == 0) {
+		   fix->size_peratom_cols > 0) {
-
-	  if (update->whichflag > 0 && 
-	      update->ntimestep % fix->peratom_freq)
-	    error->all("Fix in variable not computed at compatible time");
-	  peratom2global(1,NULL,fix->vector_atom,1,index1,
-			 tree,treestack,ntreestack,argstack,nargstack);
-
-        // f_ID[][2] = per-atom array
-
-	} else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
-		   fix->peratom_flag) {
 
 	  if (tree == NULL)
 	    error->all("Per-atom fix in equal-style variable formula");
-	  if (index2 > fix->size_peratom_cols)
+	  if (index1 > fix->size_peratom_cols)
-	    error->all("Fix vector in variable formula is too small");
+	    error->all("Variable formula fix vector is accessed out-of-range");
 	  if (update->whichflag > 0 && 
 	      update->ntimestep % fix->peratom_freq)
 	    error->all("Fix in variable not computed at compatible time");
+
 	  Tree *newtree = new Tree();
 	  newtree->type = ATOMARRAY;
-	  newtree->array = &fix->array_atom[0][index2-1];
+	  newtree->array = &fix->array_atom[0][index1-1];
 	  newtree->nstride = fix->size_peratom_cols;
 	  newtree->left = newtree->right = NULL;
 	  treestack[ntreestack++] = newtree;
 
-        // f_ID[N][2] = global value from per-atom array
-
-	} else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
-		   fix->peratom_flag) {
-
-	  if (index2 > fix->size_peratom_cols)
-	    error->all("Fix vector in variable formula is too small");
-	  if (update->whichflag > 0 && 
-	      update->ntimestep % fix->peratom_freq)
-	    error->all("Fix in variable not computed at compatible time");
-	  peratom2global(1,NULL,&fix->array_atom[0][index2-1],
-			 fix->size_peratom_cols,index1,
-			 tree,treestack,ntreestack,argstack,nargstack);
 
 	} else error->all("Mismatched fix in variable formula");
 
@@ -945,14 +1013,15 @@ double Variable::evaluate(char *str, Tree **tree)
 	// parse zero or one trailing brackets
 	// point i beyond last bracket
 	// nbracket = # of bracket pairs
-	// index = int inside bracket, 0 if empty
+	// index = int inside bracket
 
 	int nbracket,index;
 	if (str[i] != '[') nbracket = 0;
 	else {
 	  nbracket = 1;
-	  i = int_between_brackets(str,i,index,1);
+	  ptr = &str[i];
-	  i++;
+	  index = int_between_brackets(ptr);
+	  i = ptr-str+1;
 	}
 
         // v_name = non atom-style variable = global value
@@ -970,9 +1039,9 @@ double Variable::evaluate(char *str, Tree **tree)
 	    treestack[ntreestack++] = newtree;
 	  } else argstack[nargstack++] = atof(var);
 
-        // v_name[] = atom-style variable
+        // v_name = atom-style variable
 
-	} else if (nbracket && index == 0 && style[ivar] == ATOM) {
+	} else if (nbracket == 0 && style[ivar] == ATOM) {
 
 	  if (tree == NULL)
 	    error->all("Atom-style variable in equal-style variable formula");
@@ -984,7 +1053,7 @@ double Variable::evaluate(char *str, Tree **tree)
 	// compute the per-atom variable in result
 	// use peratom2global to extract single value from result
 
-	} else if (nbracket && index > 0 && style[ivar] == ATOM) {
+	} else if (nbracket && style[ivar] == ATOM) {
 
 	  double *result = (double *)
 	    memory->smalloc(atom->nlocal*sizeof(double),"variable:result");
@@ -1003,7 +1072,9 @@ double Variable::evaluate(char *str, Tree **tree)
 
       } else {
 
+	// ----------------
 	// math or group function
+	// ----------------
 
 	if (str[i] == '(') {
 	  char *contents;
@@ -1019,24 +1090,28 @@ double Variable::evaluate(char *str, Tree **tree)
 	  delete [] contents;
 
 	// ----------------
-	// atom value or vector
+	// atom value
 	// ----------------
 
 	} else if (str[i] == '[') {
-	  int id;
-	  i = int_between_brackets(str,i,id,1);
-	  i++;
-
 	  if (domain->box_exist == 0)
 	    error->all("Variable evaluation before simulation box is defined");
 
-	  // ID between brackets exists: atom value
+	  ptr = &str[i];
-	  // empty brackets: atom vector
+	  int id = int_between_brackets(ptr);
+	  i = ptr-str+1;
  
-	  if (id > 0)
 	  peratom2global(0,word,NULL,0,id,
 			 tree,treestack,ntreestack,argstack,nargstack);
-	  else
+
+	// ----------------
+	// atom vector
+	// ----------------
+
+	} else if (is_atom_vector(word)) {
+	  if (domain->box_exist == 0)
+	    error->all("Variable evaluation before simulation box is defined");
+
 	  atom_vector(word,tree,treestack,ntreestack);
 
 	// ----------------
@@ -1270,37 +1345,33 @@ int Variable::find_matching_paren(char *str, int i,char *&contents)
 }
 
 /* ----------------------------------------------------------------------
-   find int between brackets and set index to its value
+   find int between brackets and return it
-   if emptyflag, then brackets can be empty (index = 0)
+   ptr initially points to left bracket
-   else if not emptyflag and brackets empty, is an error
+   return it pointing to right bracket
-   else contents of brackets must be positive int
+   error if no right bracket or brackets are empty
-   i = left bracket
+   error if any between-bracket chars are non-digits or value == 0
-   return loc of right bracket
 ------------------------------------------------------------------------- */
 
-int Variable::int_between_brackets(char *str, int i, int &index, int emptyflag)
+int Variable::int_between_brackets(char *&ptr)
 {
-  // istop = matching ']'
+  char *start = ++ptr;
 
-  int istart = i;
+  while (*ptr && *ptr != ']') {
-  while (!str[i] || str[i] != ']') i++;
+    if (!isdigit(*ptr)) 
-  if (!str[i]) error->all("Invalid syntax in variable formula");
+      error->all("Non digit character between brackets in input command");
-  int istop = i;
+    ptr++;
-
-  int n = istop - istart - 1;
-  char *arg = new char[n+1];
-  strncpy(arg,&str[istart+1],n);
-  arg[n] = '\0';
-
-  if (n == 0 && emptyflag) index = 0;
-  else if (n == 0) error->all("Empty brackets in variable formula");
-  else {
-    index = atoi(arg);
-    if (index <= 0) error->all("Invalid index in variable formula");
   }
-  delete [] arg;
 
-  return istop;
+  if (*ptr != ']') error->all("Mismatched brackets in input command");
+  if (ptr == start) error->all("Empty brackets in input command");
+
+  *ptr = '\0';
+  int index = atoi(start);
+  *ptr = ']';
+
+  if (index == 0) 
+    error->all("Index between input command brackets must be positive");
+  return index;
 }
 
 /* ----------------------------------------------------------------------
@@ -1645,6 +1716,28 @@ void Variable::peratom2global(int flag, char *word,
   } else argstack[nargstack++] = value;
 }
 
+/* ----------------------------------------------------------------------
+   check if word matches an atom vector
+   return 1 if yes, else 0
+   customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+------------------------------------------------------------------------- */
+
+int Variable::is_atom_vector(char *word)
+{
+  if (strcmp(word,"mass") == 0) return 1;
+  if (strcmp(word,"type") == 0) return 1;
+  if (strcmp(word,"x") == 0) return 1;
+  if (strcmp(word,"y") == 0) return 1;
+  if (strcmp(word,"z") == 0) return 1;
+  if (strcmp(word,"vx") == 0) return 1;
+  if (strcmp(word,"vy") == 0) return 1;
+  if (strcmp(word,"vz") == 0) return 1;
+  if (strcmp(word,"fx") == 0) return 1;
+  if (strcmp(word,"fy") == 0) return 1;
+  if (strcmp(word,"fz") == 0) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    process an atom vector in formula
    push result onto tree
@@ -1686,6 +1779,4 @@ void Variable::atom_vector(char *word, Tree **tree,
   else if (strcmp(word,"fx") == 0) newtree->array = &atom->f[0][0];
   else if (strcmp(word,"fy") == 0) newtree->array = &atom->f[0][1];
   else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
-  
-  else error->all("Invalid atom vector in variable formula");
 }

src/variable.h

@@ -31,6 +31,7 @@ class Variable : protected Pointers {
   char *retrieve(char *);
   double compute_equal(int);
   void compute_atom(int, int, double *, int, int);
+  int int_between_brackets(char *&);
 
  private:
   int me;
@@ -59,12 +60,12 @@ class Variable : protected Pointers {
   double eval_tree(Tree *, int);
   void free_tree(Tree *);
   int find_matching_paren(char *, int, char *&);
-  int int_between_brackets(char *, int, int &, int);
   int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);
   int group_function(char *, char *, Tree **, Tree **, int &, double *, int &);
   int region_function(char *);
   void peratom2global(int, char *, double *, int, int,
 		      Tree **, Tree **, int &, double *, int &);
+  int is_atom_vector(char *);
   void atom_vector(char *, Tree **, Tree **, int &);
 };