diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index 32ff36a97a..dcceb1ac5e 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
 
 enum{SUM,MINN,MAXX,AVE};
 enum{X,V,F,COMPUTE,FIX,VARIABLE};
-enum{GLOBAL,PERATOM,LOCAL};
+enum{PERATOM,LOCAL};
 
 #define INVOKED_VECTOR 2
 #define INVOKED_ARRAY 4
@@ -181,20 +181,7 @@ 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]->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 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) {
+      if (modify->compute[icompute]->peratom_flag) {
 	flavor[i] = PERATOM;
 	if (argindex[i] == 0 && 
 	    modify->compute[icompute]->size_peratom_cols != 0)
@@ -218,26 +205,13 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
 	if (argindex[i] && 
 	    argindex[i] > modify->compute[icompute]->size_local_cols)
 	  error->all("Compute reduce compute array is accessed out-of-range");
-      }
+      } else error->all("Compute reduce compute calculates global values");
 
     } 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]->vector_flag || 
-	  modify->fix[ifix]->array_flag) {
-	flavor[i] = GLOBAL;
-	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) {
+      if (modify->fix[ifix]->peratom_flag) {
 	flavor[i] = PERATOM;
 	if (argindex[i] == 0 && 
 	    modify->fix[ifix]->size_peratom_cols != 0)
@@ -261,7 +235,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
 	if (argindex[i] && 
 	    argindex[i] > modify->fix[ifix]->size_local_cols)
 	  error->all("Compute reduce fix array is accessed out-of-range");
-      }
+      } else error->all("Compute reduce fix calculates global values");
 
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);
@@ -354,25 +328,13 @@ double ComputeReduce::compute_scalar()
   double one = compute_one(0,-1);
 
   if (mode == SUM) {
-    if (flavor[0] == GLOBAL)
-      scalar = one;
-    else
-      MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   } else if (mode == MINN) {
-    if (flavor[0] == GLOBAL)
-      scalar = one;
-    else
-      MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
+    MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
   } else if (mode == MAXX) {
-    if (flavor[0] == GLOBAL)
-      scalar = one;
-    else
-      MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
+    MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
   } else if (mode == AVE) {
-    if (flavor[0] == GLOBAL)
-      scalar = one;
-    else
-      MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
     scalar /= count(0);
   }
 
@@ -392,34 +354,22 @@ void ComputeReduce::compute_vector()
     }
 
   if (mode == SUM) {
-    for (int m = 0; m < nvalues; m++) {
-      if (flavor[m] == GLOBAL)
-	vector[m] = onevec[m];
-      else
-	MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
-    }
+    for (int m = 0; m < nvalues; m++)
+      MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
 
   } else if (mode == MINN) {
     if (!replace) {
-      for (int m = 0; m < nvalues; m++) {
-	if (flavor[m] == GLOBAL)
-	  vector[m] = onevec[m];
-	else
-	  MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MIN,world);
-      }
+      for (int m = 0; m < nvalues; m++)
+	MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MIN,world);
+
     } else {
       for (int m = 0; m < nvalues; m++)
 	if (replace[m] < 0) {
-	  if (flavor[m] == GLOBAL) {
-	    vector[m] = onevec[m];
-	    owner[m] = me;
-	  } else {
-	    pairme.value = onevec[m];
-	    pairme.proc = me;
-	    MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
-	    vector[m] = pairall.value;
-	    owner[m] = pairall.proc;
-	  }
+	  pairme.value = onevec[m];
+	  pairme.proc = me;
+	  MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
+	  vector[m] = pairall.value;
+	  owner[m] = pairall.proc;
 	}
       for (int m = 0; m < nvalues; m++)
 	if (replace[m] >= 0) {
@@ -431,25 +381,17 @@ void ComputeReduce::compute_vector()
 
   } else if (mode == MAXX) {
     if (!replace) {
-      for (int m = 0; m < nvalues; m++) {
-	if (flavor[m] == GLOBAL)
-	  vector[m] = onevec[m];
-	else
-	  MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MAX,world);
-      }
+      for (int m = 0; m < nvalues; m++)
+	MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MAX,world);
+
     } else {
       for (int m = 0; m < nvalues; m++)
 	if (replace[m] < 0) {
-	  if (flavor[m] == GLOBAL) {
-	    vector[m] = onevec[m];
-	    owner[m] = me;
-	  } else {
-	    pairme.value = onevec[m];
-	    pairme.proc = me;
-	    MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
-	    vector[m] = pairall.value;
-	    owner[m] = pairall.proc;
-	  }
+	  pairme.value = onevec[m];
+	  pairme.proc = me;
+	  MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
+	  vector[m] = pairall.value;
+	  owner[m] = pairall.proc;
 	}
       for (int m = 0; m < nvalues; m++)
 	if (replace[m] >= 0) {
@@ -461,10 +403,7 @@ void ComputeReduce::compute_vector()
 
   } else if (mode == AVE) {
     for (int m = 0; m < nvalues; m++) {
-      if (flavor[m] == GLOBAL)
-	vector[m] = onevec[m];
-      else
-	MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
+      MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
       vector[m] /= count(m);
     }
   }
@@ -524,33 +463,7 @@ double ComputeReduce::compute_one(int m, int flag)
   } else if (which[m] == COMPUTE) {
     Compute *compute = modify->compute[vidx];
     
-    if (flavor[m] == GLOBAL) {
-      if (aidx == 0) {
-	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
-	  compute->compute_vector();
-	  compute->invoked_flag |= INVOKED_VECTOR;
-	}
-	double *comp_vec = compute->vector;
-	int n = compute->size_vector;
-	if (flag < 0) 
-	  for (i = 0; i < n; i++)
-	    combine(one,comp_vec[i],i);
-	else one = comp_vec[flag];
-      } else {
-	if (!(compute->invoked_flag & INVOKED_ARRAY)) {
-	  compute->compute_array();
-	  compute->invoked_flag |= INVOKED_ARRAY;
-	}
-	double **carray = compute->array;
-	int n = compute->size_array_rows;
-	int aidx_1 = aidx - 1;
-	if (flag < 0)
-	  for (i = 0; i < n; i++)
-	    combine(one,carray[i][aidx_1],i);
-	else one = carray[flag][aidx_1];
-      }
-
-    } else if (flavor[m] == PERATOM) {
+    if (flavor[m] == PERATOM) {
       if (!(compute->invoked_flag & INVOKED_PERATOM)) {
 	compute->compute_peratom();
 	compute->invoked_flag |= INVOKED_PERATOM;
@@ -604,22 +517,7 @@ double ComputeReduce::compute_one(int m, int flag)
       error->all("Fix used in compute reduce not computed at compatible time");
     Fix *fix = modify->fix[vidx];
 
-    if (flavor[m] == GLOBAL) {
-      if (aidx == 0) {
-	int n = fix->size_vector;
-	if (flag < 0)
-	  for (i = 0; i < n; i++)
-	    combine(one,fix->compute_vector(i),i);
-	else one = fix->compute_vector(flag);
-      } else {
-	int aidxm1 = aidx - 1;
-	if (flag < 0)
-	  for (i = 0; i < nlocal; i++)
-	    combine(one,fix->compute_array(i,aidxm1),i);
-	else one = fix->compute_array(flag,aidxm1);
-      }
-
-    } else if (flavor[m] == PERATOM) {
+    if (flavor[m] == PERATOM) {
       if (aidx == 0) {
 	double *fix_vector = fix->vector_atom;
 	int n = nlocal;
@@ -686,10 +584,7 @@ double ComputeReduce::count(int m)
     return group->count(igroup);
   else if (which[m] == COMPUTE) {
     Compute *compute = modify->compute[vidx];
-    if (flavor[m] == GLOBAL) {
-      if (aidx == 0) return(1.0*compute->size_vector);
-      else return(1.0*compute->size_array_rows);
-    } else if (flavor[m] == PERATOM) {
+    if (flavor[m] == PERATOM) {
       return group->count(igroup);
     } else if (flavor[m] == LOCAL) {
       double ncount = compute->size_local_rows;
@@ -699,10 +594,7 @@ double ComputeReduce::count(int m)
     }
   } else if (which[m] == FIX) {
     Fix *fix = modify->fix[vidx];
-    if (flavor[m] == GLOBAL) {
-      if (aidx == 0) return(1.0*fix->size_vector);
-      else return(1.0*fix->size_array_rows);
-    } else if (flavor[m] == PERATOM) {
+    if (flavor[m] == PERATOM) {
       return group->count(igroup);
     } else if (flavor[m] == LOCAL) {
       double ncount = fix->size_local_rows;
diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp
index 4d9e39581d..a543b95224 100644
--- a/src/fix_ave_atom.cpp
+++ b/src/fix_ave_atom.cpp
@@ -118,9 +118,9 @@ 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 ||
-      (nrepeat-1)*nevery >= peratom_freq)
+  if (nevery <= 0 || nrepeat <= 0 || peratom_freq <= 0)
+    error->all("Illegal fix ave/atom command");
+  if (peratom_freq % nevery || (nrepeat-1)*nevery >= peratom_freq)
     error->all("Illegal fix ave/atom command");
 
   for (int i = 0; i < nvalues; i++) {
diff --git a/src/fix_ave_correlate.cpp b/src/fix_ave_correlate.cpp
new file mode 100644
index 0000000000..fe3cbb46c7
--- /dev/null
+++ b/src/fix_ave_correlate.cpp
@@ -0,0 +1,574 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: Benoit Leblanc, (Materials Design)
+                         Reese Jones (Sandia)
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "string.h"
+#include "fix_ave_correlate.h"
+#include "update.h"
+#include "modify.h"
+#include "compute.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{COMPUTE,FIX,VARIABLE};
+enum{ONE,RUNNING};
+enum{AUTO,UPPER,LOWER,AUTOUPPER,AUTOLOWER,FULL};
+
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+
+/* ---------------------------------------------------------------------- */
+
+FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
+  Fix (lmp, narg, arg)
+{
+  if (narg < 7) error->all ("Illegal fix ave/correlate command");
+
+  MPI_Comm_rank(world,&me);
+
+  nevery = atoi(arg[3]);
+  nrepeat = atoi(arg[4]);
+  nfreq = atoi(arg[5]);
+
+  global_freq = nfreq;
+  time_depend = 1;
+
+  // parse values until one isn't recognized
+
+  which = new int[narg-6];
+  argindex = new int[narg-6];
+  ids = new char*[narg-6];
+  value2index = new int[narg-6];
+  nvalues = 0;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    if (strncmp(arg[iarg],"c_",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];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Illegal fix ave/correlate command");
+	argindex[nvalues] = atoi(ptr+1);
+	*ptr = '\0';
+      } else argindex[nvalues] = 0;
+
+      n = strlen(suffix) + 1;
+      ids[nvalues] = new char[n];
+      strcpy(ids[nvalues],suffix);
+      delete [] suffix;
+
+      nvalues++;
+      iarg++;
+    } else break;
+  }
+
+  // optional args
+
+  type = AUTO;
+  ave = ONE;
+  startstep = 0;
+  prefactor = 1.0;
+  fp = NULL;
+  char *title1 = NULL;
+  char *title2 = NULL;
+  char *title3 = NULL;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"type") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
+      if (strcmp(arg[iarg+1],"auto") == 0) type = AUTO;
+      else if (strcmp(arg[iarg+1],"upper") == 0) type = UPPER;
+      else if (strcmp(arg[iarg+1],"lower") == 0) type = LOWER;
+      else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER;
+      else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER;
+      else if (strcmp(arg[iarg+1],"full") == 0) type = FULL;
+      else error->all("Illegal fix ave/correlate command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"ave") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
+      if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
+      else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
+      else error->all("Illegal fix ave/correlate command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"start") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
+      startstep = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"prefactor") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
+      prefactor = atof(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"file") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
+      if (me == 0) {
+	fp = fopen(arg[iarg+1],"w");
+	if (fp == NULL) {
+	  char str[128];
+	  sprintf(str,"Cannot open fix ave/correlate file %s",arg[iarg+1]);
+	  error->one(str);
+	}
+      }
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"title1") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/correlate 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/correlate 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/correlate 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/correlate command");
+  }
+
+  // setup and error check
+  // for fix inputs, check that fix frequency is acceptable
+
+  if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
+    error->all("Illegal fix ave/correlate command");
+  if (nfreq % nevery || (nrepeat-1)*nevery > nfreq)
+    error->all("Illegal fix ave/correlate command");
+
+  for (int i = 0; i < nvalues; i++) {
+    if (which[i] == COMPUTE) {
+      int icompute = modify->find_compute(ids[i]);
+      if (icompute < 0)
+	error->all ("Compute ID for fix ave/correlate does not exist");
+      if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0)
+	error->all ("Fix ave/correlate compute does not calculate a scalar");
+      if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
+	error->all ("Fix ave/correlate compute does not calculate a vector");
+      if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
+	error->all ("Fix ave/correlate compute vector "
+		    "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/correlate does not exist");
+      if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0)
+	error->all ("Fix ave/correlate fix does not calculate a scalar");
+      if (argindex[i] && modify->fix[ifix]->vector_flag == 0)
+	error->all ("Fix ave/correlate fix does not calculate a vector");
+      if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
+	error->all ("Fix ave/correlate fix vector is accessed out-of-range");
+      if (nevery % modify->fix[ifix]->global_freq)
+	error->all("Fix for fix ave/correlate "
+		   "not computed at compatible time");
+
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+	error->all ("Variable name for fix ave/correlate does not exist");
+      if (input->variable->equalstyle(ivariable) == 0)
+	error->all ("Fix ave/correlate variable is not equal-style variable");
+    }
+  }
+
+  // npair = # of correlation pairs to calculate
+
+  if (type == AUTO) npair = nvalues;
+  if (type == UPPER || type == LOWER) npair = nvalues*(nvalues-1)/2;
+  if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues*(nvalues+1)/2;
+  if (type == FULL) npair = nvalues*nvalues;
+
+  // print file comment lines
+
+  if (fp && me == 0) {
+    if (title1) fprintf(fp,"%s\n",title1);
+    else fprintf(fp,"# Time-correlated data for fix %s\n",id);
+    if (title2) fprintf(fp,"%s\n",title2);
+    else fprintf(fp,"# Timestep Number-of-time-windows\n");
+    if (title3) fprintf(fp,"%s\n",title3);
+    else {
+      fprintf(fp,"# Index TimeDelta Ncount");
+      if (type == AUTO)
+	for (int i = 0; i < nvalues; i++)
+	  fprintf(fp," %s*%s",arg[6+i],arg[6+i]);
+      else if (type == UPPER)
+	for (int i = 0; i < nvalues; i++)
+	  for (int j = i+1; j < nvalues; j++)
+	    fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
+      else if (type == LOWER) 
+	for (int i = 0; i < nvalues; i++)
+	  for (int j = 0; j < i-1; j++)
+	    fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
+      else if (type == AUTOUPPER)
+	for (int i = 0; i < nvalues; i++)
+	  for (int j = i; j < nvalues; j++)
+	    fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
+      else if (type == AUTOLOWER) 
+	for (int i = 0; i < nvalues; i++)
+	  for (int j = 0; j < i; j++)
+	    fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
+      else if (type == FULL)
+	for (int i = 0; i < nvalues; i++)
+	  for (int j = 0; j < nvalues; j++)
+	    fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
+      fprintf(fp,"\n");
+    }
+  }
+
+  delete [] title1;
+  delete [] title2;
+  delete [] title3;
+  
+  // allocate and initialize memory for averaging
+  // set count and corr to zero since they accumulate
+  // also set save versions to zero in case accessed via compute_array()
+
+  values = memory->create_2d_double_array(nrepeat,nvalues,
+					  "ave/correlate:values");
+  count = (int *) memory->smalloc(nrepeat*sizeof(int),
+				  "ave/correlate:count");
+  save_count = (int *) memory->smalloc(nrepeat*sizeof(int),
+				       "ave/correlate:save_count");
+  corr = memory->create_2d_double_array(nrepeat,npair,
+					"ave/correlate:corr");
+  save_corr = memory->create_2d_double_array(nrepeat,npair,
+					     "ave/correlate:save_corr");
+
+  int i,j;
+  for (i = 0; i < nrepeat; i++) {
+    save_count[i] = count[i] = 0;
+    for (j = 0; j < npair; j++)
+      save_corr[i][j] = corr[i][j] = 0.0;
+  }
+
+  lastindex = -1;
+  firstindex = 0;
+  nsample = 0;
+
+  // this fix produces a global array
+
+  array_flag = 1;
+  size_array_rows = nrepeat;
+  size_array_cols = npair+2;
+  extarray = 0;
+
+  // nvalid = next step on which end_of_step does something
+  // this step if multiple of nevery, else next multiple
+  // startstep is lower bound
+
+  nvalid = update->ntimestep;
+  if (startstep > nvalid) nvalid = startstep;
+  if (nvalid % nevery) nvalid = (nvalid/nevery)*nevery + nevery;
+
+  // add nvalid to all computes that store invocation times
+  // since don't know a priori which are invoked by this fix
+  // once in end_of_step() can set timestep for ones actually invoked
+
+  modify->addstep_compute_all(nvalid);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAveCorrelate::~FixAveCorrelate()
+{
+  delete [] which;
+  delete [] argindex;
+  delete [] value2index;
+  for (int i = 0; i < nvalues; i++) delete [] ids[i];
+  delete [] ids;
+
+  memory->destroy_2d_double_array(values);
+  memory->sfree(count);
+  memory->sfree(save_count);
+  memory->destroy_2d_double_array(corr);
+  memory->destroy_2d_double_array(save_corr);
+
+  if (fp && me == 0) fclose(fp);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAveCorrelate::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveCorrelate::init()
+{
+  // set current indices for all computes,fixes,variables
+
+  for (int i = 0; i < nvalues; i++) {
+    if (which[i] == COMPUTE) {
+      int icompute = modify->find_compute(ids[i]);
+      if (icompute < 0)
+	error->all("Compute ID for fix ave/correlate does not exist");
+      value2index[i] = icompute;
+
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0) 
+	error->all("Fix ID for fix ave/correlate does not exist");
+      value2index[i] = ifix;
+
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0) 
+	error->all("Variable name for fix ave/correlate does not exist");
+      value2index[i] = ivariable;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   only does something if nvalid = current timestep
+------------------------------------------------------------------------- */
+
+void FixAveCorrelate::setup(int vflag)
+{
+  end_of_step();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveCorrelate::end_of_step()
+{
+  int i,j,k,m;
+  double scalar;
+
+  // skip if not step which requires doing something
+
+  int ntimestep = update->ntimestep;
+  if (ntimestep != nvalid) return;
+
+  // accumulate results of computes,fixes,variables to origin
+  // compute/fix/variable may invoke computes so wrap with clear/add
+
+  modify->clearstep_compute();
+
+  // lastindex = index in values ring of latest time sample
+
+  lastindex++;
+  if (lastindex == nrepeat) lastindex = 0;
+
+  for (i = 0; i < nvalues; i++) {
+    m = value2index[i];
+    
+    // invoke compute if not previously invoked
+    
+    if (which[i] == COMPUTE) {
+      Compute *compute = modify->compute[m];
+      
+      if (argindex[i] == 0) {
+	if (!(compute->invoked_flag & INVOKED_SCALAR)) {
+	  compute->compute_scalar();
+	  compute->invoked_flag |= INVOKED_SCALAR;
+	}
+	scalar = compute->scalar;
+      } else {
+	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+	  compute->compute_vector();
+	  compute->invoked_flag |= INVOKED_VECTOR;
+	}
+	scalar = compute->vector[argindex[i]-1];
+      }
+      
+    // access fix fields, guaranteed to be ready
+      
+    } else if (which[i] == FIX) {
+      if (argindex[i] == 0) 
+	scalar = modify->fix[m]->compute_scalar();
+      else
+	scalar = modify->fix[m]->compute_vector(argindex[i]-1);
+      
+    // evaluate equal-style variable
+      
+    } else if (which[i] == VARIABLE)
+      scalar = input->variable->compute_equal(m);
+    
+    values[lastindex][i] = scalar;
+  }
+
+  // fistindex = index in values ring of earliest time sample
+  // nsample = number of time samples in values ring
+
+  if (nsample < nrepeat) nsample++;
+  else {
+    firstindex++;
+    if (firstindex == nrepeat) firstindex = 0;
+  }
+
+  nvalid += nevery;
+  modify->addstep_compute(nvalid);
+
+  // calculate all Cij() enabled by latest values
+
+  accumulate();
+  if (ntimestep % nfreq) return;
+
+  // save results in save_count and save_corr
+
+  for (i = 0; i < nrepeat; i++) {
+    save_count[i] = count[i];
+    if (count[i])
+      for (j = 0; j < npair; j++)
+	save_corr[i][j] = prefactor*corr[i][j]/count[i];
+    else
+      for (j = 0; j < npair; j++)
+	save_corr[i][j] = 0.0;
+  }
+
+  // output to file
+
+  if (fp && me == 0) {
+    fprintf(fp,"%d %d\n",ntimestep,nrepeat);
+    for (i = 0; i < nrepeat; i++) {
+      fprintf(fp,"%d %d %d",i+1,i*nevery,count[i]);
+      if (count[i])
+	for (j = 0; j < npair; j++)
+	  fprintf(fp," %g",prefactor*corr[i][j]/count[i]);
+      else
+	for (j = 0; j < npair; j++)
+	  fprintf(fp," 0.0");
+      fprintf(fp,"\n");
+    }
+    fflush(fp);
+  }
+
+  // zero accumulation if requested
+  // recalculate Cij(0)
+
+  if (ave == ONE) {
+    for (i = 0; i < nrepeat; i++) {
+      count[i] = 0;
+      for (j = 0; j < npair; j++)
+	corr[i][j] = 0.0;
+    }
+    nsample = 1;
+    accumulate();
+  }
+}
+
+/* ----------------------------------------------------------------------
+   accumulate correlation data using more recently added values
+------------------------------------------------------------------------- */
+
+void FixAveCorrelate::accumulate()
+{
+  int i,j,k,m,n,ipair;
+
+  for (k = 0; k < nsample; k++) count[k]++;
+
+  if (type == AUTO) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++) {
+	corr[k][ipair++] += values[m][i]*values[n][i];
+      }
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  } else if (type == UPPER) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++)
+	for (j = i+1; j < nvalues; j++)
+	  corr[k][ipair++] += values[m][i]*values[n][j];
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  } else if (type == LOWER) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++)
+	for (j = 0; j < i-1; j++)
+	  corr[k][ipair++] += values[m][i]*values[n][j];
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  } else if (type == AUTOUPPER) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++)
+	for (j = i; j < nvalues; j++)
+	  corr[k][ipair++] += values[m][i]*values[n][j];
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  } else if (type == AUTOLOWER) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++)
+	for (j = 0; j < i; j++)
+	  corr[k][ipair++] += values[m][i]*values[n][j];
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  } else if (type == FULL) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++)
+	for (j = 0; j < nvalues; j++)
+	  corr[k][ipair++] += values[m][i]*values[n][j];
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   return I,J array value
+------------------------------------------------------------------------- */
+
+double FixAveCorrelate::compute_array(int i, int j)
+{
+  if (j == 0) return 1.0*i*nevery;
+  else if (j == 1) return 1.0*save_count[i];
+  else if (save_count[i]) return save_corr[i][j-2];
+  return 0.0;
+}
diff --git a/src/fix_ave_correlate.h b/src/fix_ave_correlate.h
new file mode 100644
index 0000000000..f2e95b756a
--- /dev/null
+++ b/src/fix_ave_correlate.h
@@ -0,0 +1,66 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(ave/correlate,FixAveCorrelate)
+
+#else
+
+#ifndef LMP_FIX_AVE_CORRELATE_H
+#define LMP_FIX_AVE_CORRELATE_H
+
+#include "stdio.h"
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixAveCorrelate : public Fix {
+ public:
+  FixAveCorrelate(class LAMMPS *, int, char **);
+  ~FixAveCorrelate();
+  int setmask();
+  void init();
+  void setup(int);
+  void end_of_step();
+  double compute_array(int,int);
+
+ private:
+  int me,nvalues;
+  int nrepeat,nfreq,nvalid;
+  int *which,*argindex,*value2index;
+  char **ids;
+  FILE *fp;
+
+  int type,ave,startstep;
+  double prefactor;
+  char *title1,*title2,*title3;
+
+  int firstindex;      // index in values ring of earliest time sample
+  int lastindex;       // index in values ring of latest time sample
+  int nsample;         // number of time samples in values ring
+
+  int npair;           // number of correlation pairs to calculate
+  int *count;
+  double **values,**corr;
+
+  int *save_count;     // saved values at Nfreq for output via compute_array()
+  double **save_corr;
+    
+  void accumulate();
+};
+
+}
+
+#endif
+#endif
diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp
index 16073a4136..747e9a3b3c 100644
--- a/src/fix_ave_histo.cpp
+++ b/src/fix_ave_histo.cpp
@@ -238,10 +238,10 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
   // kind = inputs are all global, or all per-atom, or all local
   // for fix inputs, check that fix frequency is acceptable
 
-  if (nevery <= 0) error->all("Illegal fix ave/histo command");
-  if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
+  if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
+    error->all("Illegal fix ave/histo command");
+  if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all("Illegal fix ave/histo command");
-
   if (lo >= hi) error->all("Illegal fix ave/histo command");
   if (nbins <= 0) error->all("Illegal fix ave/histo command");
 
diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp
index eb05df32d4..7497e865ce 100644
--- a/src/fix_ave_spatial.cpp
+++ b/src/fix_ave_spatial.cpp
@@ -225,10 +225,10 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
 
   // setup and error check
 
-  if (nevery <= 0) error->all("Illegal fix ave/spatial command");
-  if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
+  if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
+    error->all("Illegal fix ave/spatial command");
+  if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all("Illegal fix ave/spatial command");
-
   if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
 
   for (int i = 0; i < nvalues; i++) {
diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp
index d07f52cab3..94dd8f56c5 100644
--- a/src/fix_ave_time.cpp
+++ b/src/fix_ave_time.cpp
@@ -165,8 +165,9 @@ 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)
+  if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
+    error->all("Illegal fix ave/time command");
+  if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all("Illegal fix ave/time command");
 
   for (int i = 0; i < nvalues; i++) {
@@ -284,7 +285,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   delete [] title2;
   delete [] title3;
 
-  // allocate and initialize memory for averaging
+  // allocate memory for averaging
 
   vector = vector_total = NULL;
   vector_list = NULL;
diff --git a/src/math_extra.h b/src/math_extra.h
index 8f6bb8ed60..5bbd9046a7 100755
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -83,7 +83,7 @@ namespace MathExtra {
 
   inline void multiply_shape_shape(const double *one, const double *two,
                                    double *ans);
-};
+}
 
 /* ----------------------------------------------------------------------
    normalize a vector
diff --git a/src/variable.cpp b/src/variable.cpp
index 9bd7db45cd..2c72d6ef5b 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -35,6 +35,9 @@ using namespace LAMMPS_NS;
 #define VARDELTA 4
 #define MAXLEVEL 4
 
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
 enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
@@ -42,14 +45,17 @@ enum{ARG,OP};
 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
        EQ,NE,LT,LE,GT,GE,AND,OR,
        SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,
-       CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,
+       CEIL,FLOOR,ROUND,
        VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
+enum{SUM,XMIN,XMAX,AVE,TRAP};
 
 #define INVOKED_SCALAR 1
 #define INVOKED_VECTOR 2
 #define INVOKED_ARRAY 4
 #define INVOKED_PERATOM 8
 
+#define BIG 1.0e20
+
 /* ---------------------------------------------------------------------- */
 
 Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
@@ -758,10 +764,10 @@ double Variable::evaluate(char *str, Tree **tree)
 	} else if (nbracket == 2 && compute->array_flag) {
 
 	  if (index1 > compute->size_array_rows)
-	    error->all("Variable formula compute vector "
+	    error->all("Variable formula compute array "
 		       "is accessed out-of-range");
 	  if (index2 > compute->size_array_cols)
-	    error->all("Variable formula compute vector "
+	    error->all("Variable formula compute array "
 		       "is accessed out-of-range");
 	  if (update->whichflag == 0) {
 	    if (compute->invoked_array != update->ntimestep)
@@ -804,7 +810,7 @@ double Variable::evaluate(char *str, Tree **tree)
 		   compute->size_peratom_cols > 0) {
 
 	  if (index2 > compute->size_peratom_cols)
-	    error->all("Variable formula compute vector "
+	    error->all("Variable formula compute array "
 		       "is accessed out-of-range");
 	  if (update->whichflag == 0) {
 	    if (compute->invoked_peratom != update->ntimestep)
@@ -850,7 +856,7 @@ double Variable::evaluate(char *str, Tree **tree)
 	  if (tree == NULL)
 	    error->all("Per-atom compute in equal-style variable formula");
 	  if (index1 > compute->size_peratom_cols)
-	    error->all("Variable formula compute vector "
+	    error->all("Variable formula compute array "
 		       "is accessed out-of-range");
 	  if (update->whichflag == 0) {
 	    if (compute->invoked_peratom != update->ntimestep)
@@ -946,9 +952,9 @@ double Variable::evaluate(char *str, Tree **tree)
 	} else if (nbracket == 2 && fix->array_flag) {
 
 	  if (index1 > fix->size_array_rows)
-	    error->all("Variable formula fix vector is accessed out-of-range");
+	    error->all("Variable formula fix array is accessed out-of-range");
 	  if (index2 > fix->size_array_cols)
-	    error->all("Variable formula fix vector is accessed out-of-range");
+	    error->all("Variable formula fix array is accessed out-of-range");
 	  if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
 	    error->all("Fix in variable not computed at compatible time");
 
@@ -979,7 +985,7 @@ double Variable::evaluate(char *str, Tree **tree)
 		   fix->size_peratom_cols > 0) {
 
 	  if (index2 > fix->size_peratom_cols)
-	    error->all("Variable formula fix vector is accessed out-of-range");
+	    error->all("Variable formula fix array is accessed out-of-range");
 	  if (update->whichflag > 0 && 
 	      update->ntimestep % fix->peratom_freq)
 	    error->all("Fix in variable not computed at compatible time");
@@ -1014,7 +1020,7 @@ double Variable::evaluate(char *str, Tree **tree)
 	  if (tree == NULL)
 	    error->all("Per-atom fix in equal-style variable formula");
 	  if (index1 > fix->size_peratom_cols)
-	    error->all("Variable formula fix vector is accessed out-of-range");
+	    error->all("Variable formula fix array is accessed out-of-range");
 	  if (update->whichflag > 0 && 
 	      update->ntimestep % fix->peratom_freq)
 	    error->all("Fix in variable not computed at compatible time");
@@ -1098,13 +1104,13 @@ double Variable::evaluate(char *str, Tree **tree)
 	delete [] id;
 
       // ----------------
-      // math/group function or atom value/vector or thermo keyword
+      // math/group/special function or atom value/vector or thermo keyword
       // ----------------
 
       } else {
 
 	// ----------------
-	// math or group function
+	// math or group or special function
 	// ----------------
 
 	if (str[i] == '(') {
@@ -1116,7 +1122,10 @@ double Variable::evaluate(char *str, Tree **tree)
 			    treestack,ntreestack,argstack,nargstack));
 	  else if (group_function(word,contents,tree,
 				  treestack,ntreestack,argstack,nargstack));
-	  else error->all("Invalid math or group function in variable formula");
+	  else if (special_function(word,contents,tree,
+				    treestack,ntreestack,argstack,nargstack));
+	  else error->all("Invalid math/group/special function "
+			  "in variable formula");
 	  delete [] contents;
 
 	// ----------------
@@ -1311,7 +1320,7 @@ double Variable::evaluate(char *str, Tree **tree)
    process an evaulation tree
    customize by adding a math function:
      sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan()
-     ceil(),floor(),round(),ramp(),stagger(),logfreq()
+     ceil(),floor(),round()
 ---------------------------------------------------------------------- */
 
 double Variable::eval_tree(Tree *tree, int i)
@@ -1424,43 +1433,6 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == ROUND)
     return MYROUND(eval_tree(tree->left,i));
 
-  if (tree->type == RAMP) {
-    if (update->whichflag == 0)
-      error->all("Cannot use ramp in variable formula between runs");
-    double arg1 = eval_tree(tree->left,i);
-    double arg2 = eval_tree(tree->right,i);
-    double delta = update->ntimestep - update->beginstep;
-    delta /= update->endstep - update->beginstep;
-    return arg1 + delta*(arg2-arg1);
-  }
-
-  if (tree->type == STAGGER) {
-    int arg1 = static_cast<int> (eval_tree(tree->left,i));
-    int arg2 = static_cast<int> (eval_tree(tree->right,i));
-    if (arg1 <= 0 || arg2 <= 0 || arg1 <= arg2)
-      error->all("Invalid math function in variable formula");
-    int lower = update->ntimestep/arg1 * arg1;
-    int delta = update->ntimestep - lower;
-    if (delta < arg2) return 1.0*(lower+arg2);
-    else return 1.0*(lower+arg1);
-  }
-
-  if (tree->type == LOGFREQ) {
-    int arg1 = static_cast<int> (eval_tree(tree->left,i));
-    int arg2 = static_cast<int> (eval_tree(tree->middle,i));
-    int arg3 = static_cast<int> (eval_tree(tree->right,i));
-    if (arg1 <= 0 || arg2 <= 0 || arg3 <= 0 || arg2 >= arg3)
-      error->all("Invalid math function in variable formula");
-    if (update->ntimestep < arg1) return 1.0*arg1;
-    else {
-      int lower = arg1;
-      while (update->ntimestep >= arg3*lower) lower *= arg3;
-      int multiple = update->ntimestep/lower;
-      if (multiple < arg2) return 1.0*(multiple+1)*lower;
-      else return 1.0*(lower*arg3);
-    }
-  }
-
   return 0.0;
 }
 
@@ -1518,19 +1490,19 @@ int Variable::int_between_brackets(char *&ptr)
 
   while (*ptr && *ptr != ']') {
     if (!isdigit(*ptr)) 
-      error->all("Non digit character between brackets in input command");
+      error->all("Non digit character between brackets in variable");
     ptr++;
   }
 
-  if (*ptr != ']') error->all("Mismatched brackets in input command");
-  if (ptr == start) error->all("Empty brackets in input command");
+  if (*ptr != ']') error->all("Mismatched brackets in variable");
+  if (ptr == start) error->all("Empty brackets in variable");
 
   *ptr = '\0';
   int index = atoi(start);
   *ptr = ']';
 
   if (index == 0) 
-    error->all("Index between input command brackets must be positive");
+    error->all("Index between variable brackets must be positive");
   return index;
 }
 
@@ -1542,7 +1514,7 @@ int Variable::int_between_brackets(char *&ptr)
    return 0 if not a match, 1 if successfully processed
    customize by adding a math function in 2 places:
      sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan()
-     ceil(),floor(),round(),ramp(),stagger(),logfreq()
+     ceil(),floor(),round()
 ------------------------------------------------------------------------- */
 
 int Variable::math_function(char *word, char *contents, Tree **tree,
@@ -1557,8 +1529,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       strcmp(word,"tan") && strcmp(word,"asin") &&
       strcmp(word,"acos") && strcmp(word,"atan") &&
       strcmp(word,"ceil") && strcmp(word,"floor") && 
-      strcmp(word,"round") && strcmp(word,"ramp") &&
-      strcmp(word,"stagger") && strcmp(word,"logfreq"))
+      strcmp(word,"round"))
     return 0;
 
   // parse contents for arg1,arg2,arg3 separated by commas
@@ -1709,50 +1680,6 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     if (narg != 1) error->all("Invalid math function in variable formula");
     if (tree) newtree->type = ROUND;
     else argstack[nargstack++] = MYROUND(value1);
-
-  } else if (strcmp(word,"ramp") == 0) {
-    if (narg != 2) error->all("Invalid math function in variable formula");
-    if (update->whichflag == 0)
-      error->all("Cannot use ramp in variable formula between runs");
-    if (tree) newtree->type = RAMP;
-    else {
-      double delta = update->ntimestep - update->beginstep;
-      delta /= update->endstep - update->beginstep;
-      argstack[nargstack++] = value1 + delta*(value2-value1);
-    }
-
-  } else if (strcmp(word,"stagger") == 0) {
-    if (narg != 2) error->all("Invalid math function in variable formula");
-    if (tree) newtree->type = STAGGER;
-    else {
-      int ivalue1 = static_cast<int> (value1);
-      int ivalue2 = static_cast<int> (value2);
-      if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
-	error->all("Invalid math function in variable formula");
-      int lower = update->ntimestep/ivalue1 * ivalue1;
-      int delta = update->ntimestep - lower;
-      if (delta < value2) argstack[nargstack++] = lower+ivalue2;
-      else argstack[nargstack++] = lower+ivalue1;
-    }
-
-  } else if (strcmp(word,"logfreq") == 0) {
-    if (narg != 3) error->all("Invalid math function in variable formula");
-    if (tree) newtree->type = LOGFREQ;
-    else {
-      int ivalue1 = static_cast<int> (value1);
-      int ivalue2 = static_cast<int> (value2);
-      int ivalue3 = static_cast<int> (value3);
-      if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
-	error->all("Invalid math function in variable formula");
-      if (update->ntimestep < ivalue1) argstack[nargstack++] = ivalue1;
-      else {
-	int lower = ivalue1;
-	while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
-	int multiple = update->ntimestep/lower;
-	if (multiple < ivalue2) argstack[nargstack++] = (multiple+1)*lower;
-	else argstack[nargstack++] = lower*ivalue3;
-      }
-    }
   }
 
   delete [] arg1;
@@ -2007,6 +1934,236 @@ int Variable::region_function(char *id)
   return iregion;
 }
 
+/* ----------------------------------------------------------------------
+   process a special function in formula
+   push result onto tree or arg stack
+   word = special function
+   contents = str between parentheses with one,two,three args
+   return 0 if not a match, 1 if successfully processed
+   customize by adding a special function in 2 places
+     ramp(x,y),stagger(x,y),logfreq(x,y,z),sum(x),min(x),max(x),
+     ave(x),trap(x)
+------------------------------------------------------------------------- */
+
+int Variable::special_function(char *word, char *contents, Tree **tree,
+			       Tree **treestack, int &ntreestack,
+			       double *argstack, int &nargstack)
+{
+  // word not a match to any special function
+
+  if (strcmp(word,"ramp") && strcmp(word,"stagger") && 
+      strcmp(word,"logfreq") && strcmp(word,"sum") &&
+      strcmp(word,"min") && strcmp(word,"max") &&
+      strcmp(word,"ave") && strcmp(word,"trap"))
+    return 0;
+
+  // parse contents for arg1,arg2,arg3 separated by commas
+  // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
+
+  char *arg1,*arg2,*arg3;
+  char *ptr1,*ptr2;
+
+  ptr1 = strchr(contents,',');
+  if (ptr1) {
+    *ptr1 = '\0';
+    ptr2 = strchr(ptr1+1,',');
+    if (ptr2) *ptr2 = '\0';
+  } else ptr2 = NULL;
+
+  int n = strlen(contents) + 1;
+  arg1 = new char[n];
+  strcpy(arg1,contents);
+  int narg = 1;
+  if (ptr1) {
+    n = strlen(ptr1+1) + 1;
+    arg2 = new char[n];
+    strcpy(arg2,ptr1+1);
+    narg = 2;
+  } else arg2 = NULL;
+  if (ptr2) {
+    n = strlen(ptr2+1) + 1;
+    arg3 = new char[n];
+    strcpy(arg3,ptr2+1);
+    narg = 3;
+  } else arg3 = NULL;
+
+  // match word to special function
+
+  double value;
+
+  if (strcmp(word,"ramp") == 0) {
+    if (narg != 2) error->all("Invalid special function in variable formula");
+    if (update->whichflag == 0)
+      error->all("Cannot use ramp in variable formula between runs");
+    double value1 = numeric(arg1);
+    double value2 = numeric(arg2);
+    double delta = update->ntimestep - update->beginstep;
+    delta /= update->endstep - update->beginstep;
+    value = value1 + delta*(value2-value1);
+
+  } else if (strcmp(word,"stagger") == 0) {
+    if (narg != 2) error->all("Invalid special function in variable formula");
+    int ivalue1 = inumeric(arg1);
+    int ivalue2 = inumeric(arg2);
+    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
+      error->all("Invalid special function in variable formula");
+    int lower = update->ntimestep/ivalue1 * ivalue1;
+    int delta = update->ntimestep - lower;
+    if (delta < ivalue2) value = lower+ivalue2;
+    else value = lower+ivalue1;
+
+  } else if (strcmp(word,"logfreq") == 0) {
+    if (narg != 3) error->all("Invalid special function in variable formula");
+    int ivalue1 = inumeric(arg1);
+    int ivalue2 = inumeric(arg2);
+    int ivalue3 = inumeric(arg3);
+    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+      error->all("Invalid special function in variable formula");
+    if (update->ntimestep < ivalue1) value = ivalue1;
+    else {
+      int lower = ivalue1;
+      while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
+      int multiple = update->ntimestep/lower;
+      if (multiple < ivalue2) value = (multiple+1)*lower;
+      else value = lower*ivalue3;
+    }
+
+  } else {
+    if (narg != 1) error->all("Invalid special function in variable formula");
+
+    int method;
+    if (strcmp(word,"sum") == 0) method = SUM;
+    else if (strcmp(word,"min") == 0) method = XMIN;
+    else if (strcmp(word,"max") == 0) method = XMAX;
+    else if (strcmp(word,"ave") == 0) method = AVE;
+    else if (strcmp(word,"trap") == 0) method = TRAP;
+
+    Compute *compute = NULL;
+    Fix *fix = NULL;
+    int index,nvec,nstride;
+
+    if (strstr(arg1,"c_") == arg1) {
+      ptr1 = strchr(arg1,'[');
+      if (ptr1) {
+	ptr2 = ptr1;
+	index = int_between_brackets(ptr2);
+	*ptr1 = '\0';
+      } else index = 0;
+
+      int icompute = modify->find_compute(&arg1[2]);
+      if (icompute < 0) error->all("Invalid compute ID in variable formula");
+      compute = modify->compute[icompute];
+      if (index == 0 && compute->vector_flag) {
+	if (update->whichflag == 0) {
+	  if (compute->invoked_vector != update->ntimestep)
+	    error->all("Compute used in variable between runs is not current");
+	} else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+	  compute->compute_vector();
+	  compute->invoked_flag |= INVOKED_VECTOR;
+	}
+	nvec = compute->size_vector;
+	nstride = 1;
+      } else if (index && compute->array_flag) {
+	if (index > compute->size_array_cols)
+	  error->all("Variable formula compute array "
+		     "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;
+	}
+	nvec = compute->size_array_rows;
+	nstride = compute->size_array_cols;
+      } else error->all("Mismatched compute in variable formula");
+
+    } else if (strstr(arg1,"f_") == arg1) {
+      ptr1 = strchr(arg1,'[');
+      if (ptr1) {
+	ptr2 = ptr1;
+	index = int_between_brackets(ptr2);
+	*ptr1 = '\0';
+      } else index = 0;
+
+      int ifix = modify->find_fix(&arg1[2]);
+      if (ifix < 0) error->all("Invalid fix ID in variable formula");
+      fix = modify->fix[ifix];
+      if (index == 0 && fix->vector_flag) {
+	if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+	  error->all("Fix in variable not computed at compatible time");
+	nvec = fix->size_vector;
+	nstride = 1;
+      } else if (index && fix->array_flag) {
+	if (index > fix->size_array_cols)
+	  error->all("Variable formula fix array is accessed out-of-range");
+	if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+	  error->all("Fix in variable not computed at compatible time");
+	nvec = fix->size_array_rows;
+	nstride = fix->size_array_cols;
+      } else error->all("Mismatched fix in variable formula");
+
+    } else error->all("Invalid special function in variable formula");
+    
+    value = 0.0;
+    if (method == XMIN) value = BIG;
+    if (method == XMAX) value = -BIG;
+
+    if (compute) {
+      double *vec;
+      if (index) vec = &compute->array[0][index];
+      else vec = compute->vector; 
+
+      int j = 0;
+      for (int i = 0; i < nvec; i++) {
+	if (method == SUM) value += vec[j];
+	else if (method == XMIN) value = MIN(value,vec[j]);
+	else if (method == XMAX) value = MAX(value,vec[j]);
+	else if (method == AVE) value += vec[j];
+	else if (method == TRAP) {
+	  if (i > 0 && i < nvec-1) value += vec[j];
+	  else value += 0.5*vec[j];
+	}
+	j += nstride;
+      }
+    }
+
+    if (fix) {
+      double one;
+      for (int i = 0; i < nvec; i++) {
+	if (index) one = fix->compute_array(i,index-1);
+	else one = fix->compute_vector(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) {
+	  if (i > 1 && i < nvec) value += one;
+	  else value += 0.5*one;
+	}
+      }
+    }
+
+    if (method == AVE) value /= nvec;
+  }
+    
+  delete [] arg1;
+  delete [] arg2;
+  delete [] arg3;
+
+  // save value in tree or on argstack
+
+  if (tree) {
+    Tree *newtree = new Tree();
+    newtree->type = VALUE;
+    newtree->value = value;
+    newtree->left = newtree->middle = newtree->right = NULL;
+    treestack[ntreestack++] = newtree;
+  } else argstack[nargstack++] = value;
+
+  return 1;
+}
+
 /* ----------------------------------------------------------------------
    extract a global value from a per-atom quantity in a formula
    flag = 0 -> word is an atom vector
@@ -2127,3 +2284,50 @@ void Variable::atom_vector(char *word, Tree **tree,
   else if (strcmp(word,"fy") == 0) newtree->array = &atom->f[0][1];
   else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
 }
+
+/* ----------------------------------------------------------------------
+   read a floating point value from a string
+   generate an error if not a legitimate floating point value
+------------------------------------------------------------------------- */
+
+double Variable::numeric(char *str)
+{
+  int n = strlen(str);
+  for (int i = 0; i < n; i++) {
+    if (isdigit(str[i])) continue;
+    if (str[i] == '-' || str[i] == '+' || str[i] == '.') continue;
+    if (str[i] == 'e' || str[i] == 'E') continue;
+    error->all("Expected floating point parameter in variable definition");
+  }
+
+  return atof(str);
+}
+
+/* ----------------------------------------------------------------------
+   read an integer value from a string
+   generate an error if not a legitimate integer value
+------------------------------------------------------------------------- */
+
+int Variable::inumeric(char *str)
+{
+  int n = strlen(str);
+  for (int i = 0; i < n; i++) {
+    if (isdigit(str[i]) || str[i] == '-' || str[i] == '+') continue;
+    error->all("Expected integer parameter in variable definition");
+  }
+
+  return atoi(str);
+}
+
+/* ----------------------------------------------------------------------
+   debug routine for printing formula tree recursively
+------------------------------------------------------------------------- */
+
+void Variable::print_tree(Tree *tree, int level)
+{
+  printf("TREE %d: %d %g\n",level,tree->type,tree->value);
+  if (tree->left) print_tree(tree->left,level+1);
+  if (tree->middle) print_tree(tree->middle,level+1);
+  if (tree->right) print_tree(tree->right,level+1);
+  return;
+}
diff --git a/src/variable.h b/src/variable.h
index 75279348c1..963fca5609 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -63,10 +63,15 @@ class Variable : protected Pointers {
   int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);
   int group_function(char *, char *, Tree **, Tree **, int &, double *, int &);
   int region_function(char *);
+  int special_function(char *, char *, Tree **, Tree **, 
+		       int &, double *, int &);
   void peratom2global(int, char *, double *, int, int,
 		      Tree **, Tree **, int &, double *, int &);
   int is_atom_vector(char *);
   void atom_vector(char *, Tree **, Tree **, int &);
+  double numeric(char *);
+  int inumeric(char *);
+  void print_tree(Tree *, int);
 };
 
 }