/* ----------------------------------------------------------------------
   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 author: Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */

#include "stdlib.h"
#include "string.h"
#include "fix_ave_cyl.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "lattice.h"
#include "modify.h"
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "math.h"

using namespace LAMMPS_NS;

enum{CENTER,ORIGIN};
enum{X,V,F,DENSITY_NUMBER,DENSITY_MASS,COMPUTE,FIX,VARIABLE};
enum{SAMPLE,ALL};
enum{BOX,LATTICE,REDUCED};
enum{ONE,RUNNING,WINDOW};

#define INVOKED_PERATOM 8

#define BIG 1000000000

/* ---------------------------------------------------------------------- */

FixAveCyl::FixAveCyl(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{
  if (narg < 10) error->all("Illegal fix ave/spatial command");

  MPI_Comm_rank(world,&me);

  nevery = atoi(arg[3]);
  nrepeat = atoi(arg[4]);
  nfreq = atoi(arg[5]);

  global_freq = nfreq;
  no_change_box = 1;
  time_depend = 1;

  if (strcmp(arg[6],"x") == 0) dim = 0;
  else if (strcmp(arg[6],"y") == 0) dim = 1;
  else if (strcmp(arg[6],"z") == 0) dim = 2;
  else error->all("Illegal fix ave/spatial command");

  dim1 = (dim + 1) % 3;
  dim2 = (dim + 2) % 3;

  if (strcmp(arg[7],"center") == 0) originflag = CENTER;
  else if (strcmp(arg[7],"origin") == 0) originflag = ORIGIN;
  else error->all("Illegal fix ave/spatial command");

  delta = atof(arg[8]);

  // parse values until one isn't recognized

  which = new int[narg-9];
  argindex = new int[narg-9];
  ids = new char*[narg-9];
  value2index = new int[narg-9];
  nvalues = 0;

  int iarg = 9;
  while (iarg < narg) {
    ids[nvalues] = NULL;

    if (strcmp(arg[iarg],"x") == 0) {
      which[nvalues] = X;
      argindex[nvalues++] = 0;
    } else if (strcmp(arg[iarg],"y") == 0) {
      which[nvalues] = X;
      argindex[nvalues++] = 1;
    } else if (strcmp(arg[iarg],"z") == 0) {
      which[nvalues] = X;
      argindex[nvalues++] = 2;

    } else if (strcmp(arg[iarg],"vx") == 0) {
      which[nvalues] = V;
      argindex[nvalues++] = 0;
    } else if (strcmp(arg[iarg],"vy") == 0) {
      which[nvalues] = V;
      argindex[nvalues++] = 1;
    } else if (strcmp(arg[iarg],"vz") == 0) {
      which[nvalues] = V;
      argindex[nvalues++] = 2;

    } else if (strcmp(arg[iarg],"fx") == 0) {
      which[nvalues] = F;
      argindex[nvalues++] = 0;
    } else if (strcmp(arg[iarg],"fy") == 0) {
      which[nvalues] = F;
      argindex[nvalues++] = 1;
    } else if (strcmp(arg[iarg],"fz") == 0) {
      which[nvalues] = F;
      argindex[nvalues++] = 2;

    } else if (strcmp(arg[iarg],"density/number") == 0) {
      which[nvalues] = DENSITY_NUMBER;
      argindex[nvalues++] = 0;
    } else if (strcmp(arg[iarg],"density/mass") == 0) {
      which[nvalues] = DENSITY_MASS;
      argindex[nvalues++] = 0;

    } else 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/spatial 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);
      nvalues++;
      delete [] suffix;

    } else break;

    iarg++;
  }

  // optional args

  normflag = ALL;
  scaleflag = LATTICE;
  fp = NULL;
  ave = ONE;
  char *title1 = NULL;
  char *title2 = NULL;
  char *title3 = NULL;

  while (iarg < narg) {
    if (strcmp(arg[iarg],"norm") == 0) {
      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
      if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL;
      else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE;
      else error->all("Illegal fix ave/spatial command");
      iarg += 2;
    } else if (strcmp(arg[iarg],"units") == 0) {
      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
      else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED;
      else error->all("Illegal fix ave/spatial command");
      iarg += 2;
    } else if (strcmp(arg[iarg],"file") == 0) {
      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
      if (me == 0) {
	fp = fopen(arg[iarg+1],"w");
	if (fp == NULL) {
	  char str[128];
	  sprintf(str,"Cannot open fix ave/spatial file %s",arg[iarg+1]);
	  error->one(str);
	}
      }
      iarg += 2;
    } else if (strcmp(arg[iarg],"ave") == 0) {
      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
      if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
      else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
      else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW;
      else error->all("Illegal fix ave/spatial command");
      if (ave == WINDOW) {
	if (iarg+3 > narg) error->all("Illegal fix ave/spatial command");
	nwindow = atoi(arg[iarg+2]);
	if (nwindow <= 0) error->all("Illegal fix ave/spatial command");
      }
      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");
  }

  // setup and error check

  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++) {
    if (which[i] == COMPUTE) {
      int icompute = modify->find_compute(ids[i]);
      if (icompute < 0)
	error->all("Compute ID for fix ave/spatial does not exist");
      if (modify->compute[icompute]->peratom_flag == 0)
	error->all("Fix ave/spatial compute does not calculate per-atom values");
      if (argindex[i] == 0 && 
	  modify->compute[icompute]->size_peratom_cols != 0)
	error->all("Fix ave/spatial compute does not "
		   "calculate a per-atom vector");
      if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
	error->all("Fix ave/spatial compute does not "
		   "calculate a per-atom array");
      if (argindex[i] && 
	  argindex[i] > modify->compute[icompute]->size_peratom_cols)
	error->all("Fix ave/spatial 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/spatial does not exist");
      if (modify->fix[ifix]->peratom_flag == 0)
	error->all("Fix ave/spatial fix does not calculate per-atom values");
      if (argindex[i] && modify->fix[ifix]->size_peratom_cols != 0)
	error->all("Fix ave/spatial fix does not calculate a per-atom vector");
      if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
	error->all("Fix ave/spatial fix does not calculate a per-atom array");
      if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
	error->all("Fix ave/spatial fix vector is accessed out-of-range");
    } else if (which[i] == VARIABLE) {
      int ivariable = input->variable->find(ids[i]);
      if (ivariable < 0)
	error->all("Variable name for fix ave/spatial does not exist");
      if (input->variable->atomstyle(ivariable) == 0)
	error->all("Fix ave/spatial variable is not atom-style variable");
    }
  }

  // print file comment lines

  if (fp && me == 0) {
    if (title1) fprintf(fp,"%s\n",title1);
    else fprintf(fp,"# Spatial-averaged data for fix %s and group %s\n",
		 id,arg[1]);
    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++) fprintf(fp," %s",arg[9+i]);
      fprintf(fp,"\n");
    }
  }

  delete [] title1;
  delete [] title2;
  delete [] title3;

  // this fix produces a global array

  array_flag = 1;
  size_local_rows = BIG;
  size_local_cols = nvalues+2;
  extarray = 0;

  // setup scaling

  int triclinic = domain->triclinic;
  if (triclinic == 1 && scaleflag != REDUCED)
    error->all("Fix ave/spatial for triclinic boxes requires units reduced");

  if (scaleflag == LATTICE && domain->lattice == NULL)
    error->all("Use of fix ave/spatial with undefined lattice");

  if (scaleflag == LATTICE) {
    xscale = domain->lattice->xlattice;
    yscale = domain->lattice->ylattice;
    zscale = domain->lattice->zlattice;
  }
  else xscale = yscale = zscale = 1.0;

  // apply scaling factors

  double scale;
  if (dim == 0) scale = xscale;
  if (dim == 1) scale = yscale;
  if (dim == 2) scale = zscale;
  delta *= scale;

  invdelta = 1.0/delta;

  // initializations

  irepeat = 0;
  iwindow = window_limit = 0;
  norm = 0;
  nlayers = maxlayer = 0;
  coord = NULL;
  layer_volume = NULL;
  count_one = count_many = count_sum = count_total = NULL;
  count_list = NULL;
  values_one = values_many = values_sum = values_total = NULL;
  values_list = NULL;

  maxatomvar = 0;
  varatom = NULL;

  maxatomlayer = 0;
  layer = NULL;

  // nvalid = next step on which end_of_step does something
  // can be this timestep if multiple of nfreq and nrepeat = 1
  // else backup from next multiple of nfreq

  nvalid = (update->ntimestep/nfreq)*nfreq + nfreq;
  if (nvalid-nfreq == update->ntimestep && nrepeat == 1)
    nvalid = update->ntimestep;
  else
    nvalid -= (nrepeat-1)*nevery;
  if (nvalid < update->ntimestep) nvalid += nfreq;

  // 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);
}

/* ---------------------------------------------------------------------- */

FixAveCyl::~FixAveCyl()
{
  delete [] which;
  delete [] argindex;
  for (int i = 0; i < nvalues; i++) delete [] ids[i];
  delete [] ids;
  delete [] value2index;

  if (fp && me == 0) fclose(fp);

  memory->sfree(coord);
  memory->sfree(layer_volume);
  memory->sfree(count_one);
  memory->sfree(count_many);
  memory->sfree(count_sum);
  memory->sfree(count_total);
  memory->destroy_2d_double_array(count_list);
  memory->destroy_2d_double_array(values_one);
  memory->destroy_2d_double_array(values_many);
  memory->destroy_2d_double_array(values_sum);
  memory->destroy_2d_double_array(values_total);
  memory->destroy_3d_double_array(values_list);

  memory->sfree(varatom);
  memory->sfree(layer);
}

/* ---------------------------------------------------------------------- */

int FixAveCyl::setmask()
{
  int mask = 0;
  mask |= END_OF_STEP;
  return mask;
}

/* ---------------------------------------------------------------------- */

void FixAveCyl::init()
{
  // # of layers cannot vary for ave = RUNNING or WINDOW

  if (ave == RUNNING || ave == WINDOW) {
    if (scaleflag != REDUCED && domain->box_change)
      error->all("Fix ave/spatial settings invalid with changing box");
  }

  // 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) {
      int icompute = modify->find_compute(ids[m]);
      if (icompute < 0)
	error->all("Compute ID for fix ave/spatial does not exist");
      value2index[m] = icompute;
      
    } else if (which[m] == FIX) {
      int ifix = modify->find_fix(ids[m]);
      if (ifix < 0) 
	error->all("Fix ID for fix ave/spatial does not exist");
      value2index[m] = ifix;

      if (nevery % modify->fix[ifix]->peratom_freq)
	error->all("Fix for fix ave/spatial not computed at compatible time");

    } else if (which[m] == VARIABLE) {
      int ivariable = input->variable->find(ids[m]);
      if (ivariable < 0) 
	error->all("Variable name for fix ave/spatial does not exist");
      value2index[m] = ivariable;

    } else value2index[m] = -1;
  }
}

/* ----------------------------------------------------------------------
   only does something if nvalid = current timestep
------------------------------------------------------------------------- */

void FixAveCyl::setup(int vflag)
{
  end_of_step();
}

/* ---------------------------------------------------------------------- */

void FixAveCyl::end_of_step()
{
  int i,j,m,n,ilayer;

  // skip if not step which requires doing something

  int ntimestep = update->ntimestep;
  if (ntimestep != nvalid) return;

  // if computing the first sample, setup layers
  // compute current origin = boundary for some layer
  // lo = layer boundary immediately below boxlo
  // hi = layer boundary immediately above boxhi
  // allocate and initialize arrays based on new layer count

  if (irepeat == 0) {
    double *boxlo,*boxhi,*prd;
    double prdmin;
    if (scaleflag == REDUCED) {
      boxlo = domain->boxlo_lamda;
      boxhi = domain->boxhi_lamda;
      prd = domain->prd_lamda;
    } else {
      boxlo = domain->boxlo;
      boxhi = domain->boxhi;
      prd = domain->prd;
    }
    
    origin1 = 0.5 * (boxlo[dim1] + boxhi[dim1]);
    origin2 = 0.5 * (boxlo[dim2] + boxhi[dim2]);
    if (originflag == ORIGIN) {
      origin1 = 0.0;
      origin2 = 0.0;
    }

    prdmin = (prd[dim1] < prd[dim2]) ? prd[dim1] : prd[dim2];
    nlayers = static_cast<int> ( prdmin * 0.5 * invdelta );

    if (nlayers > maxlayer) {
      maxlayer = nlayers;
      coord = (double *) memory->srealloc(coord,nlayers*sizeof(double),
					  "ave/spatial:coord");
      layer_volume = (double *) memory->srealloc(layer_volume,nlayers*sizeof(double),
					  "ave/spatial:coord");
      count_one = (double *) 
	memory->srealloc(count_one,nlayers*sizeof(double),
			 "ave/spatial:count_one");
      count_many = (double *) 
	memory->srealloc(count_many,nlayers*sizeof(double),
			 "ave/spatial:count_many");
      count_sum = (double *) 
	memory->srealloc(count_sum,nlayers*sizeof(double),
			 "ave/spatial:count_sum");
      count_total = (double *) 
	memory->srealloc(count_total,nlayers*sizeof(double),
			 "ave/spatial:count_total");

      values_one = memory->grow_2d_double_array(values_one,nlayers,nvalues,
						"ave/spatial:values_one");
      values_many = memory->grow_2d_double_array(values_many,nlayers,nvalues,
						 "ave/spatial:values_many");
      values_sum = memory->grow_2d_double_array(values_sum,nlayers,nvalues,
						"ave/spatial:values_sum");
      values_total = memory->grow_2d_double_array(values_total,nlayers,nvalues,
						  "ave/spatial:values_total");

      // initialize count and values total to zero since they accumulate

      for (m = 0; m < nlayers; m++) {
	for (i = 0; i < nvalues; i++) values_total[m][i] = 0.0;
	count_total[m] = 0.0;
      }

      // only allocate count and values list for ave = WINDOW
      // only happens once since nlayers never changes for these ave settings
      
      if (ave == WINDOW) {
	count_list =
	  memory->create_2d_double_array(nwindow,nlayers,
					 "ave/spatial:count_list");
	values_list =
	  memory->create_3d_double_array(nwindow,nlayers,nvalues,
					 "ave/spatial:values_list");
      }
    }

    const double PI = 4.0 * atan(1.0);
    for (m = 0; m < nlayers; m++) {
      coord[m] = (m+0.5)*delta;
      layer_volume[m] = PI*delta*delta*prd[dim]*((m+1.0)*(m+1.0)-m*m);
      count_many[m] = count_sum[m] = 0.0;
      for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0;
    }
  }
  
  // zero out arrays for one sample

  for (m = 0; m < nlayers; m++) {
    count_one[m] = 0.0;
    for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0;
  }

  // assign each atom to a layer
  // remap each atom's relevant coord back into box via PBC if necessary
  // if scaleflag = REDUCED, box coords -> lamda coords

  double **x = atom->x;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;


  if (nlocal > maxatomlayer) {
    maxatomlayer = atom->nmax;
    memory->sfree(layer);
    layer = (int *) 
      memory->smalloc(maxatomlayer*sizeof(int),"ave/spatial:layer");
  }

  double *boxlo,*boxhi,*prd;
  double x1remap,x2remap;
  int periodicity = domain->periodicity[dim];

  if (periodicity) {
    if (scaleflag == REDUCED) {
      boxlo = domain->boxlo_lamda;
      boxhi = domain->boxhi_lamda;
      prd = domain->prd_lamda;
    } else {
      boxlo = domain->boxlo;
      boxhi = domain->boxhi;
      prd = domain->prd;
    }
  }

  if (scaleflag == REDUCED) domain->x2lamda(nlocal);
  
  for (i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      x1remap = x[i][dim1];
      if (periodicity) {
	if (x1remap < boxlo[dim1]) x1remap += prd[0];
	if (x1remap >= boxhi[dim1]) x1remap -= prd[0];
      }
      x2remap = x[i][dim2];
      if (periodicity) {
	if (x2remap < boxlo[dim2]) x2remap += prd[1];
	if (x2remap >= boxhi[dim2]) x2remap -= prd[1];
      }
      
      ilayer = static_cast<int> (sqrt((x1remap-origin1)*(x1remap-origin1)+(x2remap-origin2)*(x2remap-origin2)) * invdelta);
      if (ilayer < 0) ilayer = 0;
      if (ilayer >= nlayers) ilayer = nlayers-1;
      layer[i] = ilayer;
      count_one[ilayer] += 1.0;
    }

  if (scaleflag == REDUCED) domain->lamda2x(nlocal);

  // perform the computation for one sample
  // accumulate results of attributes,computes,fixes,variables to local copy
  // sum within each layer, only include atoms in fix group
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (m = 0; m < nvalues; m++) {
    n = value2index[m];
    j = argindex[m];

    // X,V,F adds coords,velocities,forces to values

    if (which[m] == X || which[m] == V || which[m] == F) {
      double **attribute;
      if (which[m] == X) attribute = x;
      else if (which[m] == V) attribute = atom->v;
      else attribute = atom->f;

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit)
	  values_one[layer[i]][m] += attribute[i][j];

    // DENSITY_NUMBER adds 1 to values

    } else if (which[m] == DENSITY_NUMBER) {

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit)
	  values_one[layer[i]][m] += 1.0;

    // DENSITY_MASS adds mass to values

    } else if (which[m] == DENSITY_MASS) {
      int *type = atom->type;
      double *mass = atom->mass;
      double *rmass = atom->rmass;

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit)
	  if (rmass) values_one[layer[i]][m] += rmass[i];
	  else values_one[layer[i]][m] += mass[type[i]];

    // COMPUTE adds its scalar or vector component to values
    // invoke compute if not previously invoked

    } else if (which[m] == COMPUTE) {
      Compute *compute = modify->compute[n];
      if (!(compute->invoked_flag & INVOKED_PERATOM)) {
	compute->compute_peratom();
	compute->invoked_flag |= INVOKED_PERATOM;
      }
      double *vector = compute->vector_atom;
      double **array = compute->array_atom;
      int jm1 = j - 1;

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit)
	  if (j == 0) values_one[layer[i]][m] += vector[i];
	  else values_one[layer[i]][m] += array[i][jm1];
      
    // FIX adds its scalar or vector component to values
    // access fix fields, guaranteed to be ready

    } else if (which[m] == FIX) {
      double *vector = modify->fix[n]->vector_atom;
      double **array = modify->fix[n]->array_atom;
      int jm1 = j - 1;

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit) {
	  if (j == 0) values_one[layer[i]][m] += vector[i];
	  else values_one[layer[i]][m] += array[i][jm1];
	}

    // VARIABLE adds its per-atom quantities to values
    // evaluate atom-style variable

    } else if (which[m] == VARIABLE) {
      if (nlocal > maxatomvar) {
	maxatomvar = atom->nmax;
	memory->sfree(varatom);
	varatom = (double *) 
	  memory->smalloc(maxatomvar*sizeof(double),"ave/spatial:varatom");
      }

      input->variable->compute_atom(n,igroup,varatom,1,0);

      for (i = 0; i < nlocal; i++)
	if (mask[i] & groupbit)
	  values_one[layer[i]][m] += varatom[i];
    }
  }

  // process a single sample
  // if normflag = ALL, accumulate values,count separately to many
  // if normflag = SAMPLE, one = value/count, accumulate one to many
  // exception is SAMPLE density: no normalization by atom count

  if (normflag == ALL) {
    for (m = 0; m < nlayers; m++) {
      count_many[m] += count_one[m];
      for (j = 0; j < nvalues; j++)
	values_many[m][j] += values_one[m][j];
    }
  } else {
    MPI_Allreduce(count_one,count_many,nlayers,MPI_DOUBLE,MPI_SUM,world);
    for (m = 0; m < nlayers; m++) {
      if (count_many[m] > 0.0)
	for (j = 0; j < nvalues; j++) {
	  if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
	    values_many[m][j] += values_one[m][j];
	  else values_many[m][j] += values_one[m][j]/count_many[m];
	}
      count_sum[m] += count_many[m];
    }
  }

  // done if irepeat < nrepeat
  // else reset irepeat and nvalid

  irepeat++;
  if (irepeat < nrepeat) {
    nvalid += nevery;
    modify->addstep_compute(nvalid);
    return;
  }

  irepeat = 0;
  nvalid = ntimestep+nfreq - (nrepeat-1)*nevery;
  modify->addstep_compute(nvalid);

  // time average across samples
  // if normflag = ALL, final is total value / total count
  // if normflag = SAMPLE, final is sum of ave / repeat
  // exception is density: normalized by repeat, not total count

  double repeat = nrepeat;
  double mv2d = force->mv2d;

  if (normflag == ALL) {
    MPI_Allreduce(count_many,count_sum,nlayers,MPI_DOUBLE,MPI_SUM,world);
    MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nlayers*nvalues,
		  MPI_DOUBLE,MPI_SUM,world);
    for (m = 0; m < nlayers; m++) {
      if (count_sum[m] > 0.0)
	for (j = 0; j < nvalues; j++)
          if (which[j] == DENSITY_NUMBER) values_sum[m][j] /= repeat;
          else if (which[j] == DENSITY_MASS) values_sum[m][j] *= mv2d/repeat;
	  else values_sum[m][j] /= count_sum[m];
      count_sum[m] /= repeat;
    }
  } else {
    MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nlayers*nvalues,
		  MPI_DOUBLE,MPI_SUM,world);
    for (m = 0; m < nlayers; m++) {
      for (j = 0; j < nvalues; j++)
	values_sum[m][j] /= repeat;
      count_sum[m] /= repeat;
    }
  }

  // density is additionally normalized by layer volume

  for (j = 0; j < nvalues; j++)
    if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
      for (m = 0; m < nlayers; m++)
	values_sum[m][j] /= layer_volume[m];

  // if ave = ONE, only single Nfreq timestep value is needed
  // if ave = RUNNING, combine with all previous Nfreq timestep values
  // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values

  if (ave == ONE) {
    for (m = 0; m < nlayers; m++) {
      for (i = 0; i < nvalues; i++) 
	values_total[m][i] = values_sum[m][i];
      count_total[m] = count_sum[m];
    }
    norm = 1;

  } else if (ave == RUNNING) {
    for (m = 0; m < nlayers; m++) {
      for (i = 0; i < nvalues; i++)
	values_total[m][i] += values_sum[m][i];
      count_total[m] += count_sum[m];
    }
    norm++;

  } else if (ave == WINDOW) {
    for (m = 0; m < nlayers; m++) {
      for (i = 0; i < nvalues; i++) {
	values_total[m][i] += values_sum[m][i];
	if (window_limit) values_total[m][i] -= values_list[iwindow][m][i];
	values_list[iwindow][m][i] = values_sum[m][i];
      }
      count_total[m] += count_sum[m];
      if (window_limit) count_total[m] -= count_list[iwindow][m];
      count_list[iwindow][m] = count_sum[m];
    }

    iwindow++;
    if (iwindow == nwindow) {
      iwindow = 0;
      window_limit = 1;
    }
    if (window_limit) norm = nwindow;
    else norm = iwindow;
  }

  // output result to file
  
  if (fp && me == 0) {
    fprintf(fp,"%d %d\n",ntimestep,nlayers);
    for (m = 0; m < nlayers; m++) {
      fprintf(fp,"  %d %g %g",m+1,coord[m],count_total[m]/norm);
      for (i = 0; i < nvalues; i++) fprintf(fp," %g",values_total[m][i]/norm);
      fprintf(fp,"\n");
    }
    fflush(fp);
  }
}

/* ----------------------------------------------------------------------
   return I,J array value
   if I exceeds current layers, return 0.0 instead of generating an error
   1st column = bin coord, 2nd column = count, remaining columns = Nvalues
------------------------------------------------------------------------- */

double FixAveCyl::compute_array(int i, int j)
{
  if (values_total == NULL) return 0.0;
  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;
}

/* ----------------------------------------------------------------------
   memory usage of varatom and layer
------------------------------------------------------------------------- */

double FixAveCyl::memory_usage()
{
  double bytes = maxatomvar * sizeof(double);
  bytes += maxatomlayer * sizeof(int);
  return bytes;
}