diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 4e4b4224e..547dcfaa6 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -1,2409 +1,2410 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "string.h"
 #include "dump_custom.h"
 #include "atom.h"
 #include "force.h"
 #include "domain.h"
 #include "region.h"
 #include "group.h"
 #include "input.h"
 #include "variable.h"
 #include "update.h"
 #include "modify.h"
 #include "compute.h"
 #include "fix.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 // customize by adding keyword
 // also customize compute_atom_property.cpp
 
 enum{ID,MOL,TYPE,ELEMENT,MASS,
      X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,
      XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI,
      IX,IY,IZ,
      VX,VY,VZ,FX,FY,FZ,
      Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
      OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
      TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
      COMPUTE,FIX,VARIABLE};
 enum{LT,LE,GT,GE,EQ,NEQ};
 enum{INT,DOUBLE,STRING};
 
 #define INVOKED_PERATOM 8
 
 /* ---------------------------------------------------------------------- */
 
 DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
   Dump(lmp, narg, arg)
 {
   if (narg == 5) error->all(FLERR,"No dump custom arguments specified");
 
   clearstep = 1;
 
   nevery = force->inumeric(FLERR,arg[3]);
 
   // size_one may be shrunk below if additional optional args exist
 
   size_one = nfield = narg - 5;
   pack_choice = new FnPtrPack[nfield];
   vtype = new int[nfield];
 
   iregion = -1;
   idregion = NULL;
   nthresh = 0;
   thresh_array = NULL;
   thresh_op = NULL;
   thresh_value = NULL;
 
   // computes, fixes, variables which the dump accesses
 
   memory->create(field2index,nfield,"dump:field2index");
   memory->create(argindex,nfield,"dump:argindex");
 
   ncompute = 0;
   id_compute = NULL;
   compute = NULL;
 
   nfix = 0;
   id_fix = NULL;
   fix = NULL;
 
   nvariable = 0;
   id_variable = NULL;
   variable = NULL;
   vbuf = NULL;
 
   // process attributes
   // ioptional = start of additional optional args
-  // only dump image style processes optional args
+  // only dump image and dump movie styles processes optional args
 
   ioptional = parse_fields(narg,arg);
 
-  if (ioptional < narg && strcmp(style,"image") != 0)
+  if ((ioptional < narg) &&
+      (strcmp(style,"image") != 0) && (strcmp(style,"movie") != 0))
     error->all(FLERR,"Invalid attribute in dump custom command");
   size_one = nfield = ioptional - 5;
 
   // atom selection arrays
 
   maxlocal = 0;
   choose = NULL;
   dchoose = NULL;
   clist = NULL;
 
   // element names
 
   ntypes = atom->ntypes;
   typenames = NULL;
 
   // setup format strings
 
   vformat = new char*[size_one];
 
   format_default = new char[3*size_one+1];
   format_default[0] = '\0';
 
   for (int i = 0; i < size_one; i++) {
     if (vtype[i] == INT) strcat(format_default,"%d ");
     else if (vtype[i] == DOUBLE) strcat(format_default,"%g ");
     else if (vtype[i] == STRING) strcat(format_default,"%s ");
     vformat[i] = NULL;
   }
 
   // setup column string
 
   int n = 0;
   for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2;
   columns = new char[n];
   columns[0] = '\0';
   for (int iarg = 5; iarg < narg; iarg++) {
     strcat(columns,arg[iarg]);
     strcat(columns," ");
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 DumpCustom::~DumpCustom()
 {
   delete [] pack_choice;
   delete [] vtype;
   memory->destroy(field2index);
   memory->destroy(argindex);
 
   delete [] idregion;
   memory->destroy(thresh_array);
   memory->destroy(thresh_op);
   memory->destroy(thresh_value);
 
   for (int i = 0; i < ncompute; i++) delete [] id_compute[i];
   memory->sfree(id_compute);
   delete [] compute;
 
   for (int i = 0; i < nfix; i++) delete [] id_fix[i];
   memory->sfree(id_fix);
   delete [] fix;
 
   for (int i = 0; i < nvariable; i++) delete [] id_variable[i];
   memory->sfree(id_variable);
   delete [] variable;
   for (int i = 0; i < nvariable; i++) memory->destroy(vbuf[i]);
   delete [] vbuf;
 
   memory->destroy(choose);
   memory->destroy(dchoose);
   memory->destroy(clist);
 
   if (typenames) {
     for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
     delete [] typenames;
   }
 
   for (int i = 0; i < size_one; i++) delete [] vformat[i];
   delete [] vformat;
 
   delete [] columns;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::init_style()
 {
   delete [] format;
   char *str;
   if (format_user) str = format_user;
   else str = format_default;
 
   int n = strlen(str) + 1;
   format = new char[n];
   strcpy(format,str);
 
   // default for element names = C
 
   if (typenames == NULL) {
     typenames = new char*[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) {
       typenames[itype] = new char[2];
       strcpy(typenames[itype],"C");
     }
   }
 
   // tokenize the format string and add space at end of each format element
 
   char *ptr;
   for (int i = 0; i < size_one; i++) {
     if (i == 0) ptr = strtok(format," \0");
     else ptr = strtok(NULL," \0");
     if (ptr == NULL) error->all(FLERR,"Dump_modify format string is too short");
     delete [] vformat[i];
     vformat[i] = new char[strlen(ptr) + 2];
     strcpy(vformat[i],ptr);
     vformat[i] = strcat(vformat[i]," ");
   }
 
   // setup boundary string
 
   domain->boundary_string(boundstr);
 
   // setup function ptrs
 
   if (binary && domain->triclinic == 0)
     header_choice = &DumpCustom::header_binary;
   else if (binary && domain->triclinic == 1)
     header_choice = &DumpCustom::header_binary_triclinic;
   else if (!binary && domain->triclinic == 0)
     header_choice = &DumpCustom::header_item;
   else if (!binary && domain->triclinic == 1)
     header_choice = &DumpCustom::header_item_triclinic;
 
   if (binary) write_choice = &DumpCustom::write_binary;
   else write_choice = &DumpCustom::write_text;
 
   // find current ptr for each compute,fix,variable
   // check that fix frequency is acceptable
 
   int icompute;
   for (int i = 0; i < ncompute; i++) {
     icompute = modify->find_compute(id_compute[i]);
     if (icompute < 0) error->all(FLERR,"Could not find dump custom compute ID");
     compute[i] = modify->compute[icompute];
   }
 
   int ifix;
   for (int i = 0; i < nfix; i++) {
     ifix = modify->find_fix(id_fix[i]);
     if (ifix < 0) error->all(FLERR,"Could not find dump custom fix ID");
     fix[i] = modify->fix[ifix];
     if (nevery % modify->fix[ifix]->peratom_freq)
       error->all(FLERR,"Dump custom and fix not computed at compatible times");
   }
 
   int ivariable;
   for (int i = 0; i < nvariable; i++) {
     ivariable = input->variable->find(id_variable[i]);
     if (ivariable < 0)
       error->all(FLERR,"Could not find dump custom variable name");
     variable[i] = ivariable;
   }
 
   // set index and check validity of region
 
   if (iregion >= 0) {
     iregion = domain->find_region(idregion);
     if (iregion == -1)
       error->all(FLERR,"Region ID for dump custom does not exist");
   }
 
   // open single file, one time only
 
   if (multifile == 0) openfile();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::write_header(bigint ndump)
 {
   if (multiproc) (this->*header_choice)(ndump);
   else if (me == 0) (this->*header_choice)(ndump);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::header_binary(bigint ndump)
 {
   fwrite(&update->ntimestep,sizeof(bigint),1,fp);
   fwrite(&ndump,sizeof(bigint),1,fp);
   fwrite(&domain->triclinic,sizeof(int),1,fp);
   fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp);
   fwrite(&boxxlo,sizeof(double),1,fp);
   fwrite(&boxxhi,sizeof(double),1,fp);
   fwrite(&boxylo,sizeof(double),1,fp);
   fwrite(&boxyhi,sizeof(double),1,fp);
   fwrite(&boxzlo,sizeof(double),1,fp);
   fwrite(&boxzhi,sizeof(double),1,fp);
   fwrite(&size_one,sizeof(int),1,fp);
   if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp);
   else fwrite(&nprocs,sizeof(int),1,fp);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::header_binary_triclinic(bigint ndump)
 {
   fwrite(&update->ntimestep,sizeof(bigint),1,fp);
   fwrite(&ndump,sizeof(bigint),1,fp);
   fwrite(&domain->triclinic,sizeof(int),1,fp);
   fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp);
   fwrite(&boxxlo,sizeof(double),1,fp);
   fwrite(&boxxhi,sizeof(double),1,fp);
   fwrite(&boxylo,sizeof(double),1,fp);
   fwrite(&boxyhi,sizeof(double),1,fp);
   fwrite(&boxzlo,sizeof(double),1,fp);
   fwrite(&boxzhi,sizeof(double),1,fp);
   fwrite(&boxxy,sizeof(double),1,fp);
   fwrite(&boxxz,sizeof(double),1,fp);
   fwrite(&boxyz,sizeof(double),1,fp);
   fwrite(&size_one,sizeof(int),1,fp);
   if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp);
   else fwrite(&nprocs,sizeof(int),1,fp);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::header_item(bigint ndump)
 {
   fprintf(fp,"ITEM: TIMESTEP\n");
   fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep);
   fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
   fprintf(fp,BIGINT_FORMAT "\n",ndump);
   fprintf(fp,"ITEM: BOX BOUNDS %s\n",boundstr);
   fprintf(fp,"%g %g\n",boxxlo,boxxhi);
   fprintf(fp,"%g %g\n",boxylo,boxyhi);
   fprintf(fp,"%g %g\n",boxzlo,boxzhi);
   fprintf(fp,"ITEM: ATOMS %s\n",columns);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::header_item_triclinic(bigint ndump)
 {
   fprintf(fp,"ITEM: TIMESTEP\n");
   fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep);
   fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
   fprintf(fp,BIGINT_FORMAT "\n",ndump);
   fprintf(fp,"ITEM: BOX BOUNDS xy xz yz %s\n",boundstr);
   fprintf(fp,"%g %g %g\n",boxxlo,boxxhi,boxxy);
   fprintf(fp,"%g %g %g\n",boxylo,boxyhi,boxxz);
   fprintf(fp,"%g %g %g\n",boxzlo,boxzhi,boxyz);
   fprintf(fp,"ITEM: ATOMS %s\n",columns);
 }
 
 /* ---------------------------------------------------------------------- */
 
 int DumpCustom::count()
 {
   int i;
 
   // grow choose and variable vbuf arrays if needed
 
   int nlocal = atom->nlocal;
   if (nlocal > maxlocal) {
     maxlocal = atom->nmax;
 
     memory->destroy(choose);
     memory->destroy(dchoose);
     memory->destroy(clist);
     memory->create(choose,maxlocal,"dump:choose");
     memory->create(dchoose,maxlocal,"dump:dchoose");
     memory->create(clist,maxlocal,"dump:clist");
 
     for (i = 0; i < nvariable; i++) {
       memory->destroy(vbuf[i]);
       memory->create(vbuf[i],maxlocal,"dump:vbuf");
     }
   }
 
   // invoke Computes for per-atom quantities
 
   if (ncompute) {
     for (i = 0; i < ncompute; i++)
       if (!(compute[i]->invoked_flag & INVOKED_PERATOM)) {
         compute[i]->compute_peratom();
         compute[i]->invoked_flag |= INVOKED_PERATOM;
       }
   }
 
   // evaluate atom-style Variables for per-atom quantities
 
   if (nvariable)
     for (i = 0; i < nvariable; i++)
       input->variable->compute_atom(variable[i],igroup,vbuf[i],1,0);
 
   // choose all local atoms for output
 
   for (i = 0; i < nlocal; i++) choose[i] = 1;
 
   // un-choose if not in group
 
   if (igroup) {
     int *mask = atom->mask;
     for (i = 0; i < nlocal; i++)
       if (!(mask[i] & groupbit))
         choose[i] = 0;
   }
 
   // un-choose if not in region
 
   if (iregion >= 0) {
     Region *region = domain->regions[iregion];
     double **x = atom->x;
     for (i = 0; i < nlocal; i++)
       if (choose[i] && region->match(x[i][0],x[i][1],x[i][2]) == 0)
         choose[i] = 0;
   }
 
   // un-choose if any threshhold criterion isn't met
 
   if (nthresh) {
     double *ptr;
     double value;
     int nstride;
     int nlocal = atom->nlocal;
 
     for (int ithresh = 0; ithresh < nthresh; ithresh++) {
 
       // customize by adding to if statement
 
       if (thresh_array[ithresh] == ID) {
         int *tag = atom->tag;
         for (i = 0; i < nlocal; i++) dchoose[i] = tag[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == MOL) {
         if (!atom->molecule_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         int *molecule = atom->molecule;
         for (i = 0; i < nlocal; i++) dchoose[i] = molecule[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == TYPE) {
         int *type = atom->type;
         for (i = 0; i < nlocal; i++) dchoose[i] = type[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ELEMENT) {
         int *type = atom->type;
         for (i = 0; i < nlocal; i++) dchoose[i] = type[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == MASS) {
         if (atom->rmass) {
           ptr = atom->rmass;
           nstride = 1;
         } else {
           double *mass = atom->mass;
           int *type = atom->type;
           for (i = 0; i < nlocal; i++) dchoose[i] = mass[type[i]];
           ptr = dchoose;
           nstride = 1;
         }
 
       } else if (thresh_array[ithresh] == X) {
         ptr = &atom->x[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == Y) {
         ptr = &atom->x[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == Z) {
         ptr = &atom->x[0][2];
         nstride = 3;
 
       } else if (thresh_array[ithresh] == XS) {
         double **x = atom->x;
         double boxxlo = domain->boxlo[0];
         double invxprd = 1.0/domain->xprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (x[i][0] - boxxlo) * invxprd;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == YS) {
         double **x = atom->x;
         double boxylo = domain->boxlo[1];
         double invyprd = 1.0/domain->yprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (x[i][1] - boxylo) * invyprd;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ZS) {
         double **x = atom->x;
         double boxzlo = domain->boxlo[2];
         double invzprd = 1.0/domain->zprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (x[i][2] - boxzlo) * invzprd;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == XSTRI) {
         double **x = atom->x;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[0]*(x[i][0]-boxlo[0]) +
             h_inv[5]*(x[i][1]-boxlo[1]) + h_inv[4]*(x[i][2]-boxlo[2]);
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == YSTRI) {
         double **x = atom->x;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[1]*(x[i][1]-boxlo[1]) +
             h_inv[3]*(x[i][2]-boxlo[2]);
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ZSTRI) {
         double **x = atom->x;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[2]*(x[i][2]-boxlo[2]);
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == XU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double xprd = domain->xprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = x[i][0] + ((image[i] & IMGMASK) - IMGMAX) * xprd;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == YU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double yprd = domain->yprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = x[i][1] + 
             ((image[i] >> IMGBITS & IMGMASK) - IMGMAX) * yprd;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ZU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double zprd = domain->zprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = x[i][2] + ((image[i] >> IMG2BITS) - IMGMAX) * zprd;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == XUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *h = domain->h;
         int xbox,ybox,zbox;
         for (i = 0; i < nlocal; i++) {
           xbox = (image[i] & IMGMASK) - IMGMAX;
           ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
           zbox = (image[i] >> IMG2BITS) - IMGMAX;
           dchoose[i] = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
         }
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == YUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *h = domain->h;
         int ybox,zbox;
         for (i = 0; i < nlocal; i++) {
           ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
           zbox = (image[i] >> IMG2BITS) - IMGMAX;
           dchoose[i] = x[i][1] + h[1]*ybox + h[3]*zbox;
         }
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ZUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *h = domain->h;
         int zbox;
         for (i = 0; i < nlocal; i++) {
           zbox = (image[i] >> IMG2BITS) - IMGMAX;
           dchoose[i] = x[i][2] + h[2]*zbox;
         }
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == XSU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double boxxlo = domain->boxlo[0];
         double invxprd = 1.0/domain->xprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (x[i][0] - boxxlo) * invxprd + (image[i] & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == YSU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double boxylo = domain->boxlo[1];
         double invyprd = 1.0/domain->yprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] =
             (x[i][1] - boxylo) * invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == ZSU) {
         double **x = atom->x;
         tagint *image = atom->image;
         double boxzlo = domain->boxlo[2];
         double invzprd = 1.0/domain->zprd;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (x[i][2] - boxzlo) * invzprd + (image[i] >> IMG2BITS) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == XSUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[0]*(x[i][0]-boxlo[0]) +
             h_inv[5]*(x[i][1]-boxlo[1]) +
             h_inv[4]*(x[i][2]-boxlo[2]) +
             (image[i] & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == YSUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[1]*(x[i][1]-boxlo[1]) +
             h_inv[3]*(x[i][2]-boxlo[2]) +
             (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ZSUTRI) {
         double **x = atom->x;
         tagint *image = atom->image;
         double *boxlo = domain->boxlo;
         double *h_inv = domain->h_inv;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = h_inv[2]*(x[i][2]-boxlo[2]) +
             (image[i] >> IMG2BITS) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == IX) {
         tagint *image = atom->image;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (image[i] & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == IY) {
         tagint *image = atom->image;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == IZ) {
         tagint *image = atom->image;
         for (i = 0; i < nlocal; i++)
           dchoose[i] = (image[i] >> IMG2BITS) - IMGMAX;
         ptr = dchoose;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == VX) {
         ptr = &atom->v[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == VY) {
         ptr = &atom->v[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == VZ) {
         ptr = &atom->v[0][2];
         nstride = 3;
       } else if (thresh_array[ithresh] == FX) {
         ptr = &atom->f[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == FY) {
         ptr = &atom->f[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == FZ) {
         ptr = &atom->f[0][2];
         nstride = 3;
 
       } else if (thresh_array[ithresh] == Q) {
         if (!atom->q_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = atom->q;
         nstride = 1;
       } else if (thresh_array[ithresh] == MUX) {
         if (!atom->mu_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->mu[0][0];
         nstride = 4;
       } else if (thresh_array[ithresh] == MUY) {
         if (!atom->mu_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->mu[0][1];
         nstride = 4;
       } else if (thresh_array[ithresh] == MUZ) {
         if (!atom->mu_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->mu[0][2];
         nstride = 4;
       } else if (thresh_array[ithresh] == MU) {
         if (!atom->mu_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->mu[0][3];
         nstride = 4;
 
       } else if (thresh_array[ithresh] == RADIUS) {
         if (!atom->radius_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = atom->radius;
         nstride = 1;
       } else if (thresh_array[ithresh] == DIAMETER) {
         if (!atom->radius_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         double *radius = atom->radius;
         for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == OMEGAX) {
         if (!atom->omega_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->omega[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == OMEGAY) {
         if (!atom->omega_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->omega[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == OMEGAZ) {
         if (!atom->omega_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->omega[0][2];
         nstride = 3;
       } else if (thresh_array[ithresh] == ANGMOMX) {
         if (!atom->angmom_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->angmom[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == ANGMOMY) {
         if (!atom->angmom_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->angmom[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == ANGMOMZ) {
         if (!atom->angmom_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->angmom[0][2];
         nstride = 3;
       } else if (thresh_array[ithresh] == TQX) {
         if (!atom->torque_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->torque[0][0];
         nstride = 3;
       } else if (thresh_array[ithresh] == TQY) {
         if (!atom->torque_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->torque[0][1];
         nstride = 3;
       } else if (thresh_array[ithresh] == TQZ) {
         if (!atom->torque_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = &atom->torque[0][2];
         nstride = 3;
 
       } else if (thresh_array[ithresh] == SPIN) {
         if (!atom->spin_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         int *spin = atom->spin;
         for (i = 0; i < nlocal; i++) dchoose[i] = spin[i];
         ptr = dchoose;
         nstride = 1;
       } else if (thresh_array[ithresh] == ERADIUS) {
         if (!atom->eradius_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = atom->eradius;
         nstride = 1;
       } else if (thresh_array[ithresh] == ERVEL) {
         if (!atom->ervel_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = atom->ervel;
         nstride = 1;
       } else if (thresh_array[ithresh] == ERFORCE) {
         if (!atom->erforce_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
         ptr = atom->erforce;
         nstride = 1;
 
       } else if (thresh_array[ithresh] == COMPUTE) {
         i = nfield + ithresh;
         if (argindex[i] == 0) {
           ptr = compute[field2index[i]]->vector_atom;
           nstride = 1;
         } else {
           ptr = &compute[field2index[i]]->array_atom[0][argindex[i]-1];
           nstride = compute[field2index[i]]->size_peratom_cols;
         }
 
       } else if (thresh_array[ithresh] == FIX) {
         i = nfield + ithresh;
         if (argindex[i] == 0) {
           ptr = fix[field2index[i]]->vector_atom;
           nstride = 1;
         } else {
           ptr = &fix[field2index[i]]->array_atom[0][argindex[i]-1];
           nstride = fix[field2index[i]]->size_peratom_cols;
         }
 
       } else if (thresh_array[ithresh] == VARIABLE) {
         i = nfield + ithresh;
         ptr = vbuf[field2index[i]];
         nstride = 1;
       }
 
       // unselect atoms that don't meet threshhold criterion
 
       value = thresh_value[ithresh];
 
       if (thresh_op[ithresh] == LT) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr >= value) choose[i] = 0;
       } else if (thresh_op[ithresh] == LE) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr > value) choose[i] = 0;
       } else if (thresh_op[ithresh] == GT) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr <= value) choose[i] = 0;
       } else if (thresh_op[ithresh] == GE) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr < value) choose[i] = 0;
       } else if (thresh_op[ithresh] == EQ) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr != value) choose[i] = 0;
       } else if (thresh_op[ithresh] == NEQ) {
         for (i = 0; i < nlocal; i++, ptr += nstride)
           if (choose[i] && *ptr == value) choose[i] = 0;
       }
     }
   }
 
   // compress choose flags into clist
   // nchoose = # of selected atoms
   // clist[i] = local index of each selected atom
 
   nchoose = 0;
   for (i = 0; i < nlocal; i++)
     if (choose[i]) clist[nchoose++] = i;
 
   return nchoose;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack(int *ids)
 {
   for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n);
   if (ids) {
     int *tag = atom->tag;
     for (int i = 0; i < nchoose; i++)
       ids[i] = tag[clist[i]];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::write_data(int n, double *mybuf)
 {
   (this->*write_choice)(n,mybuf);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::write_binary(int n, double *mybuf)
 {
   n *= size_one;
   fwrite(&n,sizeof(int),1,fp);
   fwrite(mybuf,sizeof(double),n,fp);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::write_text(int n, double *mybuf)
 {
   int i,j;
 
   int m = 0;
   for (i = 0; i < n; i++) {
     for (j = 0; j < size_one; j++) {
       if (vtype[j] == INT) fprintf(fp,vformat[j],static_cast<int> (mybuf[m]));
       else if (vtype[j] == DOUBLE) fprintf(fp,vformat[j],mybuf[m]);
       else if (vtype[j] == STRING)
         fprintf(fp,vformat[j],typenames[(int) mybuf[m]]);
       m++;
     }
     fprintf(fp,"\n");
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int DumpCustom::parse_fields(int narg, char **arg)
 {
   // customize by adding to if statement
 
   int i;
   for (int iarg = 5; iarg < narg; iarg++) {
     i = iarg-5;
 
     if (strcmp(arg[iarg],"id") == 0) {
       pack_choice[i] = &DumpCustom::pack_id;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"mol") == 0) {
       if (!atom->molecule_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_molecule;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"type") == 0) {
       pack_choice[i] = &DumpCustom::pack_type;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"element") == 0) {
       pack_choice[i] = &DumpCustom::pack_type;
       vtype[i] = STRING;
     } else if (strcmp(arg[iarg],"mass") == 0) {
       pack_choice[i] = &DumpCustom::pack_mass;
       vtype[i] = DOUBLE;
 
     } else if (strcmp(arg[iarg],"x") == 0) {
       pack_choice[i] = &DumpCustom::pack_x;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"y") == 0) {
       pack_choice[i] = &DumpCustom::pack_y;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"z") == 0) {
       pack_choice[i] = &DumpCustom::pack_z;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"xs") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xs_triclinic;
       else pack_choice[i] = &DumpCustom::pack_xs;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"ys") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_ys_triclinic;
       else pack_choice[i] = &DumpCustom::pack_ys;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"zs") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zs_triclinic;
       else pack_choice[i] = &DumpCustom::pack_zs;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"xu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_xu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"yu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_yu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_yu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"zu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_zu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"xsu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xsu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_xsu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"ysu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_ysu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_ysu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"zsu") == 0) {
       if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zsu_triclinic;
       else pack_choice[i] = &DumpCustom::pack_zsu;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"ix") == 0) {
       pack_choice[i] = &DumpCustom::pack_ix;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"iy") == 0) {
       pack_choice[i] = &DumpCustom::pack_iy;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"iz") == 0) {
       pack_choice[i] = &DumpCustom::pack_iz;
       vtype[i] = INT;
 
     } else if (strcmp(arg[iarg],"vx") == 0) {
       pack_choice[i] = &DumpCustom::pack_vx;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"vy") == 0) {
       pack_choice[i] = &DumpCustom::pack_vy;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"vz") == 0) {
       pack_choice[i] = &DumpCustom::pack_vz;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"fx") == 0) {
       pack_choice[i] = &DumpCustom::pack_fx;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"fy") == 0) {
       pack_choice[i] = &DumpCustom::pack_fy;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"fz") == 0) {
       pack_choice[i] = &DumpCustom::pack_fz;
       vtype[i] = DOUBLE;
 
     } else if (strcmp(arg[iarg],"q") == 0) {
       if (!atom->q_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_q;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"mux") == 0) {
       if (!atom->mu_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_mux;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"muy") == 0) {
       if (!atom->mu_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_muy;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"muz") == 0) {
       if (!atom->mu_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_muz;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"mu") == 0) {
       if (!atom->mu_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_mu;
       vtype[i] = DOUBLE;
 
     } else if (strcmp(arg[iarg],"radius") == 0) {
       if (!atom->radius_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_radius;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"diameter") == 0) {
       if (!atom->radius_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_diameter;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"omegax") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_omegax;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"omegay") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_omegay;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"omegaz") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_omegaz;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"angmomx") == 0) {
       if (!atom->angmom_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_angmomx;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"angmomy") == 0) {
       if (!atom->angmom_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_angmomy;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"angmomz") == 0) {
       if (!atom->angmom_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_angmomz;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"tqx") == 0) {
       if (!atom->torque_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_tqx;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"tqy") == 0) {
       if (!atom->torque_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_tqy;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"tqz") == 0) {
       if (!atom->torque_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_tqz;
       vtype[i] = DOUBLE;
 
     } else if (strcmp(arg[iarg],"spin") == 0) {
       if (!atom->spin_flag)
         error->all(FLERR,"Dumping an atom quantity that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_spin;
       vtype[i] = INT;
     } else if (strcmp(arg[iarg],"eradius") == 0) {
       if (!atom->eradius_flag)
         error->all(FLERR,"Dumping an atom quantity that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_eradius;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"ervel") == 0) {
       if (!atom->ervel_flag)
         error->all(FLERR,"Dumping an atom quantity that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_ervel;
       vtype[i] = DOUBLE;
     } else if (strcmp(arg[iarg],"erforce") == 0) {
       if (!atom->erforce_flag)
         error->all(FLERR,"Dumping an atom quantity that isn't allocated");
       pack_choice[i] = &DumpCustom::pack_erforce;
       vtype[i] = DOUBLE;
 
     // compute value = c_ID
     // 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;
       vtype[i] = DOUBLE;
 
       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(FLERR,"Invalid attribute in dump custom command");
         argindex[i] = atoi(ptr+1);
         *ptr = '\0';
       } else argindex[i] = 0;
 
       n = modify->find_compute(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump custom compute ID");
       if (modify->compute[n]->peratom_flag == 0)
         error->all(FLERR,"Dump custom compute does not compute per-atom info");
       if (argindex[i] == 0 && modify->compute[n]->size_peratom_cols > 0)
         error->all(FLERR,"Dump custom compute does not calculate per-atom vector");
       if (argindex[i] > 0 && modify->compute[n]->size_peratom_cols == 0)
         error->all(FLERR,"Dump custom compute does not calculate per-atom array");
       if (argindex[i] > 0 &&
           argindex[i] > modify->compute[n]->size_peratom_cols)
         error->all(FLERR,"Dump custom compute vector is accessed out-of-range");
 
       field2index[i] = add_compute(suffix);
       delete [] suffix;
 
     // fix value = f_ID
     // if no trailing [], then arg is set to 0, else arg is between []
 
     } else if (strncmp(arg[iarg],"f_",2) == 0) {
       pack_choice[i] = &DumpCustom::pack_fix;
       vtype[i] = DOUBLE;
 
       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(FLERR,"Invalid attribute in dump custom command");
         argindex[i] = atoi(ptr+1);
         *ptr = '\0';
       } else argindex[i] = 0;
 
       n = modify->find_fix(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump custom fix ID");
       if (modify->fix[n]->peratom_flag == 0)
         error->all(FLERR,"Dump custom fix does not compute per-atom info");
       if (argindex[i] == 0 && modify->fix[n]->size_peratom_cols > 0)
         error->all(FLERR,"Dump custom fix does not compute per-atom vector");
       if (argindex[i] > 0 && modify->fix[n]->size_peratom_cols == 0)
         error->all(FLERR,"Dump custom fix does not compute per-atom array");
       if (argindex[i] > 0 &&
           argindex[i] > modify->fix[n]->size_peratom_cols)
         error->all(FLERR,"Dump custom fix vector is accessed out-of-range");
 
       field2index[i] = add_fix(suffix);
       delete [] suffix;
 
     // variable value = v_name
 
     } else if (strncmp(arg[iarg],"v_",2) == 0) {
       pack_choice[i] = &DumpCustom::pack_variable;
       vtype[i] = DOUBLE;
 
       int n = strlen(arg[iarg]);
       char *suffix = new char[n];
       strcpy(suffix,&arg[iarg][2]);
 
       argindex[i] = 0;
 
       n = input->variable->find(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump custom variable name");
       if (input->variable->atomstyle(n) == 0)
         error->all(FLERR,"Dump custom variable is not atom-style variable");
 
       field2index[i] = add_variable(suffix);
       delete [] suffix;
 
     } else return iarg;
   }
 
   return narg;
 }
 
 /* ----------------------------------------------------------------------
    add Compute to list of Compute objects used by dump
    return index of where this Compute is in list
    if already in list, do not add, just return index, else add to list
 ------------------------------------------------------------------------- */
 
 int DumpCustom::add_compute(char *id)
 {
   int icompute;
   for (icompute = 0; icompute < ncompute; icompute++)
     if (strcmp(id,id_compute[icompute]) == 0) break;
   if (icompute < ncompute) return icompute;
 
   id_compute = (char **)
     memory->srealloc(id_compute,(ncompute+1)*sizeof(char *),"dump:id_compute");
   delete [] compute;
   compute = new Compute*[ncompute+1];
 
   int n = strlen(id) + 1;
   id_compute[ncompute] = new char[n];
   strcpy(id_compute[ncompute],id);
   ncompute++;
   return ncompute-1;
 }
 
 /* ----------------------------------------------------------------------
    add Fix to list of Fix objects used by dump
    return index of where this Fix is in list
    if already in list, do not add, just return index, else add to list
 ------------------------------------------------------------------------- */
 
 int DumpCustom::add_fix(char *id)
 {
   int ifix;
   for (ifix = 0; ifix < nfix; ifix++)
     if (strcmp(id,id_fix[ifix]) == 0) break;
   if (ifix < nfix) return ifix;
 
   id_fix = (char **)
     memory->srealloc(id_fix,(nfix+1)*sizeof(char *),"dump:id_fix");
   delete [] fix;
   fix = new Fix*[nfix+1];
 
   int n = strlen(id) + 1;
   id_fix[nfix] = new char[n];
   strcpy(id_fix[nfix],id);
   nfix++;
   return nfix-1;
 }
 
 /* ----------------------------------------------------------------------
    add Variable to list of Variables used by dump
    return index of where this Variable is in list
    if already in list, do not add, just return index, else add to list
 ------------------------------------------------------------------------- */
 
 int DumpCustom::add_variable(char *id)
 {
   int ivariable;
   for (ivariable = 0; ivariable < nvariable; ivariable++)
     if (strcmp(id,id_variable[ivariable]) == 0) break;
   if (ivariable < nvariable) return ivariable;
 
   id_variable = (char **)
     memory->srealloc(id_variable,(nvariable+1)*sizeof(char *),
                      "dump:id_variable");
   delete [] variable;
   variable = new int[nvariable+1];
   delete [] vbuf;
   vbuf = new double*[nvariable+1];
   for (int i = 0; i <= nvariable; i++) vbuf[i] = NULL;
 
   int n = strlen(id) + 1;
   id_variable[nvariable] = new char[n];
   strcpy(id_variable[nvariable],id);
   nvariable++;
   return nvariable-1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int DumpCustom::modify_param(int narg, char **arg)
 {
   if (strcmp(arg[0],"region") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     if (strcmp(arg[1],"none") == 0) iregion = -1;
     else {
       iregion = domain->find_region(arg[1]);
       if (iregion == -1)
         error->all(FLERR,"Dump_modify region ID does not exist");
       int n = strlen(arg[1]) + 1;
       idregion = new char[n];
       strcpy(idregion,arg[1]);
     }
     return 2;
   }
 
   if (strcmp(arg[0],"element") == 0) {
     if (narg < ntypes+1)
       error->all(FLERR,"Dump modify element names do not match atom types");
 
     if (typenames) {
       for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
       delete [] typenames;
       typenames = NULL;
     }
 
     typenames = new char*[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) {
       int n = strlen(arg[itype]) + 1;
       typenames[itype] = new char[n];
       strcpy(typenames[itype],arg[itype]);
     }
     return ntypes+1;
   }
 
   if (strcmp(arg[0],"thresh") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     if (strcmp(arg[1],"none") == 0) {
       if (nthresh) {
         memory->destroy(thresh_array);
         memory->destroy(thresh_op);
         memory->destroy(thresh_value);
         thresh_array = NULL;
         thresh_op = NULL;
         thresh_value = NULL;
       }
       nthresh = 0;
       return 2;
     }
 
     if (narg < 4) error->all(FLERR,"Illegal dump_modify command");
 
     // grow threshhold arrays
 
     memory->grow(thresh_array,nthresh+1,"dump:thresh_array");
     memory->grow(thresh_op,(nthresh+1),"dump:thresh_op");
     memory->grow(thresh_value,(nthresh+1),"dump:thresh_value");
 
     // set attribute type of threshhold
     // customize by adding to if statement
 
     if (strcmp(arg[1],"id") == 0) thresh_array[nthresh] = ID;
     else if (strcmp(arg[1],"mol") == 0) thresh_array[nthresh] = MOL;
     else if (strcmp(arg[1],"type") == 0) thresh_array[nthresh] = TYPE;
     else if (strcmp(arg[1],"mass") == 0) thresh_array[nthresh] = MASS;
 
     else if (strcmp(arg[1],"x") == 0) thresh_array[nthresh] = X;
     else if (strcmp(arg[1],"y") == 0) thresh_array[nthresh] = Y;
     else if (strcmp(arg[1],"z") == 0) thresh_array[nthresh] = Z;
 
     else if (strcmp(arg[1],"xs") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = XS;
     else if (strcmp(arg[1],"xs") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = XSTRI;
     else if (strcmp(arg[1],"ys") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = YS;
     else if (strcmp(arg[1],"ys") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = YSTRI;
     else if (strcmp(arg[1],"zs") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = ZS;
     else if (strcmp(arg[1],"zs") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = ZSTRI;
 
     else if (strcmp(arg[1],"xu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = XU;
     else if (strcmp(arg[1],"xu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = XUTRI;
     else if (strcmp(arg[1],"yu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = YU;
     else if (strcmp(arg[1],"yu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = YUTRI;
     else if (strcmp(arg[1],"zu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = ZU;
     else if (strcmp(arg[1],"zu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = ZUTRI;
 
     else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = XSU;
     else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = XSUTRI;
     else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = YSU;
     else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = YSUTRI;
     else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 0)
       thresh_array[nthresh] = ZSU;
     else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 1)
       thresh_array[nthresh] = ZSUTRI;
 
     else if (strcmp(arg[1],"ix") == 0) thresh_array[nthresh] = IX;
     else if (strcmp(arg[1],"iy") == 0) thresh_array[nthresh] = IY;
     else if (strcmp(arg[1],"iz") == 0) thresh_array[nthresh] = IZ;
     else if (strcmp(arg[1],"vx") == 0) thresh_array[nthresh] = VX;
     else if (strcmp(arg[1],"vy") == 0) thresh_array[nthresh] = VY;
     else if (strcmp(arg[1],"vz") == 0) thresh_array[nthresh] = VZ;
     else if (strcmp(arg[1],"fx") == 0) thresh_array[nthresh] = FX;
     else if (strcmp(arg[1],"fy") == 0) thresh_array[nthresh] = FY;
     else if (strcmp(arg[1],"fz") == 0) thresh_array[nthresh] = FZ;
 
     else if (strcmp(arg[1],"q") == 0) thresh_array[nthresh] = Q;
     else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX;
     else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY;
     else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ;
     else if (strcmp(arg[1],"mu") == 0) thresh_array[nthresh] = MU;
 
     else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
     else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER;
     else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
     else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;
     else if (strcmp(arg[1],"omegaz") == 0) thresh_array[nthresh] = OMEGAZ;
     else if (strcmp(arg[1],"angmomx") == 0) thresh_array[nthresh] = ANGMOMX;
     else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY;
     else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ;
     else if (strcmp(arg[1],"tqx") == 0) thresh_array[nthresh] = TQX;
     else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY;
     else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ;
 
     else if (strcmp(arg[1],"spin") == 0) thresh_array[nthresh] = SPIN;
     else if (strcmp(arg[1],"eradius") == 0) thresh_array[nthresh] = ERADIUS;
     else if (strcmp(arg[1],"ervel") == 0) thresh_array[nthresh] = ERVEL;
     else if (strcmp(arg[1],"erforce") == 0) thresh_array[nthresh] = ERFORCE;
 
     // compute value = c_ID
     // if no trailing [], then arg is set to 0, else arg is between []
     // must grow field2index and argindex arrays, since access is beyond nfield
 
     else if (strncmp(arg[1],"c_",2) == 0) {
       thresh_array[nthresh] = COMPUTE;
       memory->grow(field2index,nfield+nthresh+1,"dump:field2index");
       memory->grow(argindex,nfield+nthresh+1,"dump:argindex");
       int n = strlen(arg[1]);
       char *suffix = new char[n];
       strcpy(suffix,&arg[1][2]);
 
       char *ptr = strchr(suffix,'[');
       if (ptr) {
         if (suffix[strlen(suffix)-1] != ']')
           error->all(FLERR,"Invalid attribute in dump modify command");
         argindex[nfield+nthresh] = atoi(ptr+1);
         *ptr = '\0';
       } else argindex[nfield+nthresh] = 0;
 
       n = modify->find_compute(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump modify compute ID");
 
       if (modify->compute[n]->peratom_flag == 0)
         error->all(FLERR,
                    "Dump modify compute ID does not compute per-atom info");
       if (argindex[nfield+nthresh] == 0 &&
           modify->compute[n]->size_peratom_cols > 0)
         error->all(FLERR,
                    "Dump modify compute ID does not compute per-atom vector");
       if (argindex[nfield+nthresh] > 0 &&
           modify->compute[n]->size_peratom_cols == 0)
         error->all(FLERR,
                    "Dump modify compute ID does not compute per-atom array");
       if (argindex[nfield+nthresh] > 0 &&
           argindex[nfield+nthresh] > modify->compute[n]->size_peratom_cols)
         error->all(FLERR,"Dump modify compute ID vector is not large enough");
 
       field2index[nfield+nthresh] = add_compute(suffix);
       delete [] suffix;
 
     // fix value = f_ID
     // if no trailing [], then arg is set to 0, else arg is between []
     // must grow field2index and argindex arrays, since access is beyond nfield
 
     } else if (strncmp(arg[1],"f_",2) == 0) {
       thresh_array[nthresh] = FIX;
       memory->grow(field2index,nfield+nthresh+1,"dump:field2index");
       memory->grow(argindex,nfield+nthresh+1,"dump:argindex");
       int n = strlen(arg[1]);
       char *suffix = new char[n];
       strcpy(suffix,&arg[1][2]);
 
       char *ptr = strchr(suffix,'[');
       if (ptr) {
         if (suffix[strlen(suffix)-1] != ']')
           error->all(FLERR,"Invalid attribute in dump modify command");
         argindex[nfield+nthresh] = atoi(ptr+1);
         *ptr = '\0';
       } else argindex[nfield+nthresh] = 0;
 
       n = modify->find_fix(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump modify fix ID");
 
       if (modify->fix[n]->peratom_flag == 0)
         error->all(FLERR,"Dump modify fix ID does not compute per-atom info");
       if (argindex[nfield+nthresh] == 0 &&
           modify->fix[n]->size_peratom_cols > 0)
         error->all(FLERR,"Dump modify fix ID does not compute per-atom vector");
       if (argindex[nfield+nthresh] > 0 &&
           modify->fix[n]->size_peratom_cols == 0)
         error->all(FLERR,"Dump modify fix ID does not compute per-atom array");
       if (argindex[nfield+nthresh] > 0 &&
           argindex[nfield+nthresh] > modify->fix[n]->size_peratom_cols)
         error->all(FLERR,"Dump modify fix ID vector is not large enough");
 
       field2index[nfield+nthresh] = add_fix(suffix);
       delete [] suffix;
 
     // variable value = v_ID
     // must grow field2index and argindex arrays, since access is beyond nfield
 
     } else if (strncmp(arg[1],"v_",2) == 0) {
       thresh_array[nthresh] = VARIABLE;
       memory->grow(field2index,nfield+nthresh+1,"dump:field2index");
       memory->grow(argindex,nfield+nthresh+1,"dump:argindex");
       int n = strlen(arg[1]);
       char *suffix = new char[n];
       strcpy(suffix,&arg[1][2]);
 
       argindex[nfield+nthresh] = 0;
 
       n = input->variable->find(suffix);
       if (n < 0) error->all(FLERR,"Could not find dump modify variable name");
       if (input->variable->atomstyle(n) == 0)
         error->all(FLERR,"Dump modify variable is not atom-style variable");
 
       field2index[nfield+nthresh] = add_variable(suffix);
       delete [] suffix;
 
     } else error->all(FLERR,"Invalid dump_modify threshhold operator");
 
     // set operation type of threshhold
 
     if (strcmp(arg[2],"<") == 0) thresh_op[nthresh] = LT;
     else if (strcmp(arg[2],"<=") == 0) thresh_op[nthresh] = LE;
     else if (strcmp(arg[2],">") == 0) thresh_op[nthresh] = GT;
     else if (strcmp(arg[2],">=") == 0) thresh_op[nthresh] = GE;
     else if (strcmp(arg[2],"==") == 0) thresh_op[nthresh] = EQ;
     else if (strcmp(arg[2],"!=") == 0) thresh_op[nthresh] = NEQ;
     else error->all(FLERR,"Invalid dump_modify threshhold operator");
 
     // set threshhold value
 
     thresh_value[nthresh] = force->numeric(FLERR,arg[3]);
 
     nthresh++;
     return 4;
   }
 
   return 0;
 }
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory in buf, choose, variable arrays
 ------------------------------------------------------------------------- */
 
 bigint DumpCustom::memory_usage()
 {
   bigint bytes = Dump::memory_usage();
   bytes += memory->usage(choose,maxlocal);
   bytes += memory->usage(dchoose,maxlocal);
   bytes += memory->usage(clist,maxlocal);
   bytes += memory->usage(vbuf,nvariable,maxlocal);
   return bytes;
 }
 
 /* ----------------------------------------------------------------------
    extraction of Compute, Fix, Variable results
 ------------------------------------------------------------------------- */
 
 void DumpCustom::pack_compute(int n)
 {
   double *vector = compute[field2index[n]]->vector_atom;
   double **array = compute[field2index[n]]->array_atom;
   int index = argindex[n];
 
   if (index == 0) {
     for (int i = 0; i < nchoose; i++) {
       buf[n] = vector[clist[i]];
       n += size_one;
     }
   } else {
     index--;
     for (int i = 0; i < nchoose; i++) {
       buf[n] = array[clist[i]][index];
       n += size_one;
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_fix(int n)
 {
   double *vector = fix[field2index[n]]->vector_atom;
   double **array = fix[field2index[n]]->array_atom;
   int index = argindex[n];
 
   if (index == 0) {
     for (int i = 0; i < nchoose; i++) {
       buf[n] = vector[clist[i]];
       n += size_one;
     }
   } else {
     index--;
     for (int i = 0; i < nchoose; i++) {
       buf[n] = array[clist[i]][index];
       n += size_one;
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_variable(int n)
 {
   double *vector = vbuf[field2index[n]];
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = vector[clist[i]];
     n += size_one;
   }
 }
 
 /* ----------------------------------------------------------------------
    one method for every attribute dump custom can output
    the atom property is packed into buf starting at n with stride size_one
    customize a new attribute by adding a method
 ------------------------------------------------------------------------- */
 
 void DumpCustom::pack_id(int n)
 {
   int *tag = atom->tag;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = tag[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_molecule(int n)
 {
   int *molecule = atom->molecule;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = molecule[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_type(int n)
 {
   int *type = atom->type;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = type[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_mass(int n)
 {
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
 
   if (rmass) {
     for (int i = 0; i < nchoose; i++) {
       buf[n] = rmass[clist[i]];
       n += size_one;
     }
   } else {
     for (int i = 0; i < nchoose; i++) {
       buf[n] = mass[type[clist[i]]];
       n += size_one;
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_x(int n)
 {
   double **x = atom->x;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = x[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_y(int n)
 {
   double **x = atom->x;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = x[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_z(int n)
 {
   double **x = atom->x;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = x[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xs(int n)
 {
   double **x = atom->x;
 
   double boxxlo = domain->boxlo[0];
   double invxprd = 1.0/domain->xprd;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (x[clist[i]][0] - boxxlo) * invxprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ys(int n)
 {
   double **x = atom->x;
 
   double boxylo = domain->boxlo[1];
   double invyprd = 1.0/domain->yprd;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (x[clist[i]][1] - boxylo) * invyprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zs(int n)
 {
   double **x = atom->x;
 
   double boxzlo = domain->boxlo[2];
   double invzprd = 1.0/domain->zprd;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (x[clist[i]][2] - boxzlo) * invzprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xs_triclinic(int n)
 {
   int j;
   double **x = atom->x;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = h_inv[0]*(x[j][0]-boxlo[0]) + h_inv[5]*(x[j][1]-boxlo[1]) +
       h_inv[4]*(x[j][2]-boxlo[2]);
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ys_triclinic(int n)
 {
   int j;
   double **x = atom->x;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = h_inv[1]*(x[j][1]-boxlo[1]) + h_inv[3]*(x[j][2]-boxlo[2]);
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zs_triclinic(int n)
 {
   double **x = atom->x;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = h_inv[2]*(x[clist[i]][2]-boxlo[2]);
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double xprd = domain->xprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = x[j][0] + ((image[j] & IMGMASK) - IMGMAX) * xprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_yu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double yprd = domain->yprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = x[j][1] + ((image[j] >> IMGBITS & IMGMASK) - IMGMAX) * yprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double zprd = domain->zprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = x[j][2] + ((image[j] >> IMG2BITS) - IMGMAX) * zprd;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *h = domain->h;
   int xbox,ybox,zbox;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     xbox = (image[j] & IMGMASK) - IMGMAX;
     ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
     zbox = (image[j] >> IMG2BITS) - IMGMAX;
     buf[n] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_yu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *h = domain->h;
   int ybox,zbox;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
     zbox = (image[j] >> IMG2BITS) - IMGMAX;
     buf[n] = x[j][1] + h[1]*ybox + h[3]*zbox;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *h = domain->h;
   int zbox;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     zbox = (image[j] >> IMG2BITS) - IMGMAX;
     buf[n] = x[j][2] + h[2]*zbox;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xsu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double boxxlo = domain->boxlo[0];
   double invxprd = 1.0/domain->xprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = (x[j][0] - boxxlo) * invxprd + (image[j] & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ysu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double boxylo = domain->boxlo[1];
   double invyprd = 1.0/domain->yprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = (x[j][1] - boxylo) * invyprd + (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zsu(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double boxzlo = domain->boxlo[2];
   double invzprd = 1.0/domain->zprd;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = (x[j][2] - boxzlo) * invzprd + (image[j] >> IMG2BITS) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_xsu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = h_inv[0]*(x[j][0]-boxlo[0]) + h_inv[5]*(x[j][1]-boxlo[1]) +
       h_inv[4]*(x[j][2]-boxlo[2]) + (image[j] & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ysu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = h_inv[1]*(x[j][1]-boxlo[1]) + h_inv[3]*(x[j][2]-boxlo[2]) +
       (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_zsu_triclinic(int n)
 {
   int j;
   double **x = atom->x;
   tagint *image = atom->image;
 
   double *boxlo = domain->boxlo;
   double *h_inv = domain->h_inv;
 
   for (int i = 0; i < nchoose; i++) {
     j = clist[i];
     buf[n] = h_inv[2]*(x[j][2]-boxlo[2]) + (image[j] >> IMG2BITS) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ix(int n)
 {
   tagint *image = atom->image;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (image[clist[i]] & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_iy(int n)
 {
   tagint *image = atom->image;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (image[clist[i]] >> IMGBITS & IMGMASK) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_iz(int n)
 {
   tagint *image = atom->image;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = (image[clist[i]] >> IMG2BITS) - IMGMAX;
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_vx(int n)
 {
   double **v = atom->v;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = v[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_vy(int n)
 {
   double **v = atom->v;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = v[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_vz(int n)
 {
   double **v = atom->v;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = v[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_fx(int n)
 {
   double **f = atom->f;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = f[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_fy(int n)
 {
   double **f = atom->f;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = f[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_fz(int n)
 {
   double **f = atom->f;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = f[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_q(int n)
 {
   double *q = atom->q;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = q[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_mux(int n)
 {
   double **mu = atom->mu;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = mu[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_muy(int n)
 {
   double **mu = atom->mu;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = mu[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_muz(int n)
 {
   double **mu = atom->mu;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = mu[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_mu(int n)
 {
   double **mu = atom->mu;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = mu[clist[i]][3];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_radius(int n)
 {
   double *radius = atom->radius;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = radius[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_diameter(int n)
 {
   double *radius = atom->radius;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = 2.0*radius[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_omegax(int n)
 {
   double **omega = atom->omega;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = omega[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_omegay(int n)
 {
   double **omega = atom->omega;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = omega[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_omegaz(int n)
 {
   double **omega = atom->omega;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = omega[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_angmomx(int n)
 {
   double **angmom = atom->angmom;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = angmom[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_angmomy(int n)
 {
   double **angmom = atom->angmom;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = angmom[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_angmomz(int n)
 {
   double **angmom = atom->angmom;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = angmom[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_tqx(int n)
 {
   double **torque = atom->torque;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = torque[clist[i]][0];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_tqy(int n)
 {
   double **torque = atom->torque;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = torque[clist[i]][1];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_tqz(int n)
 {
   double **torque = atom->torque;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = torque[clist[i]][2];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_spin(int n)
 {
   int *spin = atom->spin;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = spin[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_eradius(int n)
 {
   double *eradius = atom->eradius;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = eradius[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_ervel(int n)
 {
   double *ervel = atom->ervel;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = ervel[clist[i]];
     n += size_one;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpCustom::pack_erforce(int n)
 {
   double *erforce = atom->erforce;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = erforce[clist[i]];
     n += size_one;
   }
 }
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index a0063ff30..51c77158e 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -1,1068 +1,1070 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "ctype.h"
 #include "stdlib.h"
 #include "string.h"
 #include "dump_image.h"
 #include "image.h"
 #include "atom.h"
 #include "domain.h"
 #include "group.h"
 #include "force.h"
 #include "comm.h"
 #include "input.h"
 #include "variable.h"
 #include "math_const.h"
 #include "error.h"
 #include "memory.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
 #define BIG 1.0e20
 
-enum{PPM,JPG,PNG};
 enum{NUMERIC,ATOM,TYPE,ELEMENT,ATTRIBUTE};
 enum{STATIC,DYNAMIC};
 enum{NO,YES};
 
 /* ---------------------------------------------------------------------- */
 
 DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
   DumpCustom(lmp, narg, arg)
 {
   if (binary || multiproc) error->all(FLERR,"Invalid dump image filename");
 
   // force binary flag on to avoid corrupted output on Windows
 
   binary = 1;
 
   // set filetype based on filename suffix
 
   int n = strlen(filename);
   if (strlen(filename) > 4 && strcmp(&filename[n-4],".jpg") == 0)
     filetype = JPG;
   else if (strlen(filename) > 4 && strcmp(&filename[n-4],".JPG") == 0)
     filetype = JPG;
   else if (strlen(filename) > 5 && strcmp(&filename[n-5],".jpeg") == 0)
     filetype = JPG;
   else if (strlen(filename) > 5 && strcmp(&filename[n-5],".JPEG") == 0)
     filetype = JPG;
   else if (strlen(filename) > 4 && strcmp(&filename[n-4],".png") == 0)
     filetype = PNG;
   else if (strlen(filename) > 4 && strcmp(&filename[n-4],".PNG") == 0)
     filetype = PNG;
   else filetype = PPM;
 
 #ifndef LAMMPS_JPEG
   if (filetype == JPG)
     error->all(FLERR,"Support for writing images in JPEG format not included");
 #endif
 #ifndef LAMMPS_PNG
   if (filetype == PNG)
     error->all(FLERR,"Support for writing images in PNG format not included");
 #endif
 
   // atom color,diameter settings
 
   if (nfield != 2) error->all(FLERR,"Illegal dump image command");
 
   acolor = ATTRIBUTE;
   if (strcmp(arg[5],"type") == 0) acolor = TYPE;
   else if (strcmp(arg[5],"element") == 0) acolor = ELEMENT;
 
   adiam = ATTRIBUTE;
   if (strcmp(arg[6],"type") == 0) adiam = TYPE;
   else if (strcmp(arg[6],"element") == 0) adiam = ELEMENT;
 
   // create Image class with single colormap for atoms
   // change defaults for 2d
 
   image = new Image(lmp,1);
 
   if (domain->dimension == 2) {
     image->theta = 0.0;
     image->phi = 0.0;
     image->up[0] = 0.0; image->up[1] = 1.0; image->up[2] = 0.0;
   }
 
   // set defaults for optional args
 
   atomflag = YES;
   if (atom->nbondtypes == 0) bondflag = NO;
   else {
     bondflag = YES;
     bcolor = ATOM;
     bdiam = NUMERIC;
     bdiamvalue = 0.5;
   }
 
   thetastr = phistr = NULL;
   cflag = STATIC;
   cx = cy = cz = 0.5;
   cxstr = cystr = czstr = NULL;
 
   upxstr = upystr = upzstr = NULL;
   zoomstr = NULL;
   perspstr = NULL;
   boxflag = YES;
   boxdiam = 0.02;
   axesflag = NO;
 
   // parse optional args
 
   int iarg = ioptional;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"adiam") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       adiam = NUMERIC;
       adiamvalue = force->numeric(FLERR,arg[iarg+1]);
       if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"atom") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES;
       else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO;
       else error->all(FLERR,"Illegal dump image command");
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"bond") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       if (atom->nbondtypes == 0)
         error->all(FLERR,"Dump image bond not allowed with no bond types");
       bondflag = YES;
       if (strcmp(arg[iarg+1],"none") == 0) bondflag = NO;
       else if (strcmp(arg[iarg+1],"atom") == 0) bcolor = ATOM;
       else if (strcmp(arg[iarg+1],"type") == 0) bcolor = TYPE;
       else error->all(FLERR,"Illegal dump image command");
       if (!islower(arg[iarg+2][0])) {
           bdiam = NUMERIC;
           bdiamvalue = force->numeric(FLERR,arg[iarg+2]);
           if (bdiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
       } else if (strcmp(arg[iarg+2],"atom") == 0) bdiam = ATOM;
       else if (strcmp(arg[iarg+2],"type") == 0) bdiam = TYPE;
       else if (strcmp(arg[iarg+2],"none") == 0) bondflag = NO;
       else error->all(FLERR,"Illegal dump image command");
       iarg += 3;
 
     } else if (strcmp(arg[iarg],"size") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       int width = force->inumeric(FLERR,arg[iarg+1]);
       int height = force->inumeric(FLERR,arg[iarg+2]);
       if (width <= 0 || height <= 0)
         error->all(FLERR,"Illegal dump image command");
       image->width = width;
       image->height = height;
       iarg += 3;
 
     } else if (strcmp(arg[iarg],"view") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         int n = strlen(&arg[iarg+1][2]) + 1;
         thetastr = new char[n];
         strcpy(thetastr,&arg[iarg+1][2]);
       } else {
         double theta = force->numeric(FLERR,arg[iarg+1]);
         if (theta < 0.0 || theta > 180.0)
           error->all(FLERR,"Invalid dump image theta value");
         theta *= MY_PI/180.0;
         image->theta = theta;
       }
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
         int n = strlen(&arg[iarg+2][2]) + 1;
         phistr = new char[n];
         strcpy(phistr,&arg[iarg+2][2]);
       } else {
         double phi = force->numeric(FLERR,arg[iarg+2]);
         phi *= MY_PI/180.0;
         image->phi = phi;
       }
       iarg += 3;
 
     } else if (strcmp(arg[iarg],"center") == 0) {
       if (iarg+5 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"s") == 0) cflag = STATIC;
       else if (strcmp(arg[iarg+1],"d") == 0) cflag = DYNAMIC;
       else error->all(FLERR,"Illegal dump image command");
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
         int n = strlen(&arg[iarg+2][2]) + 1;
         cxstr = new char[n];
         strcpy(cxstr,&arg[iarg+2][2]);
         cflag = DYNAMIC;
       } else cx = force->numeric(FLERR,arg[iarg+2]);
       if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
         int n = strlen(&arg[iarg+3][2]) + 1;
         cystr = new char[n];
         strcpy(cystr,&arg[iarg+3][2]);
         cflag = DYNAMIC;
       } else cy = force->numeric(FLERR,arg[iarg+3]);
       if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) {
         int n = strlen(&arg[iarg+4][2]) + 1;
         czstr = new char[n];
         strcpy(czstr,&arg[iarg+4][2]);
         cflag = DYNAMIC;
       } else cz = force->numeric(FLERR,arg[iarg+4]);
       iarg += 5;
 
     } else if (strcmp(arg[iarg],"up") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         int n = strlen(&arg[iarg+1][2]) + 1;
         upxstr = new char[n];
         strcpy(upxstr,&arg[iarg+1][2]);
       } else image->up[0] = force->numeric(FLERR,arg[iarg+1]);
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
         int n = strlen(&arg[iarg+2][2]) + 1;
         upystr = new char[n];
         strcpy(upystr,&arg[iarg+2][2]);
       } else image->up[1] = force->numeric(FLERR,arg[iarg+2]);
       if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
         int n = strlen(&arg[iarg+3][2]) + 1;
         upzstr = new char[n];
         strcpy(upzstr,&arg[iarg+3][2]);
       } else image->up[2] = force->numeric(FLERR,arg[iarg+3]);
       iarg += 4;
 
     } else if (strcmp(arg[iarg],"zoom") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         int n = strlen(&arg[iarg+1][2]) + 1;
         zoomstr = new char[n];
         strcpy(zoomstr,&arg[iarg+1][2]);
       } else {
         double zoom = force->numeric(FLERR,arg[iarg+1]);
         if (zoom <= 0.0) error->all(FLERR,"Illegal dump image command");
         image->zoom = zoom;
       }
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"persp") == 0) {
       error->all(FLERR,"Dump image persp option is not yet supported");
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         int n = strlen(&arg[iarg+1][2]) + 1;
         perspstr = new char[n];
         strcpy(perspstr,&arg[iarg+1][2]);
       } else {
         double persp = force->numeric(FLERR,arg[iarg+1]);
         if (persp < 0.0) error->all(FLERR,"Illegal dump image command");
         image->persp = persp;
       }
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"box") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"yes") == 0) boxflag = YES;
       else if (strcmp(arg[iarg+1],"no") == 0) boxflag = NO;
       else error->all(FLERR,"Illegal dump image command");
       boxdiam = force->numeric(FLERR,arg[iarg+2]);
       if (boxdiam < 0.0) error->all(FLERR,"Illegal dump image command");
       iarg += 3;
 
     } else if (strcmp(arg[iarg],"axes") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"yes") == 0) axesflag = YES;
       else if (strcmp(arg[iarg+1],"no") == 0) axesflag = NO;
       else error->all(FLERR,"Illegal dump image command");
       axeslen = force->numeric(FLERR,arg[iarg+2]);
       axesdiam = force->numeric(FLERR,arg[iarg+3]);
       if (axeslen < 0.0 || axesdiam < 0.0)
         error->all(FLERR,"Illegal dump image command");
       iarg += 4;
 
     } else if (strcmp(arg[iarg],"shiny") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       double shiny = force->numeric(FLERR,arg[iarg+1]);
       if (shiny < 0.0 || shiny > 1.0)
         error->all(FLERR,"Illegal dump image command");
       image->shiny = shiny;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"ssao") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"yes") == 0) image->ssao = YES;
       else if (strcmp(arg[iarg+1],"no") == 0) image->ssao = NO;
       else error->all(FLERR,"Illegal dump image command");
       int seed = force->inumeric(FLERR,arg[iarg+2]);
       if (seed <= 0) error->all(FLERR,"Illegal dump image command");
       image->seed = seed;
       double ssaoint = force->numeric(FLERR,arg[iarg+3]);
       if (ssaoint < 0.0 || ssaoint > 1.0)
         error->all(FLERR,"Illegal dump image command");
       image->ssaoint = ssaoint;
       iarg += 4;
 
     } else error->all(FLERR,"Illegal dump image command");
   }
 
   // allocate image buffer now that image size is known
 
   image->buffers();
 
   // communication neede for bonds colored by atoms
 
   if (bondflag) {
     if (bcolor == ATOM || bdiam == ATOM) comm_forward = 3;
     else comm_forward = 1;
   }
 
   // additional defaults for dump_modify options
 
   diamtype = new double[ntypes+1];
   diamelement = new double[ntypes+1];
   colortype = new double*[ntypes+1];
   colorelement = new double*[ntypes+1];
 
   for (int i = 1; i <= ntypes; i++) {
     diamtype[i] = 1.0;
     if (i % 6 == 1) colortype[i] = image->color2rgb("red");
     else if (i % 6 == 2) colortype[i] = image->color2rgb("green");
     else if (i % 6 == 3) colortype[i] = image->color2rgb("blue");
     else if (i % 6 == 4) colortype[i] = image->color2rgb("yellow");
     else if (i % 6 == 5) colortype[i] = image->color2rgb("aqua");
     else if (i % 6 == 0) colortype[i] = image->color2rgb("cyan");
   }
 
   if (bondflag) {
     bdiamtype = new double[atom->nbondtypes+1];
     bcolortype = new double*[atom->nbondtypes+1];
     for (int i = 1; i <= atom->nbondtypes; i++) {
       bdiamtype[i] = 0.5;
       if (i % 6 == 1) bcolortype[i] = image->color2rgb("red");
       else if (i % 6 == 2) bcolortype[i] = image->color2rgb("green");
       else if (i % 6 == 3) bcolortype[i] = image->color2rgb("blue");
       else if (i % 6 == 4) bcolortype[i] = image->color2rgb("yellow");
       else if (i % 6 == 5) bcolortype[i] = image->color2rgb("aqua");
       else if (i % 6 == 0) bcolortype[i] = image->color2rgb("cyan");
     }
   } else {
     bdiamtype = NULL;
     bcolortype = NULL;
   }
 
   // viewflag = DYNAMIC if any view parameter is dynamic
 
   viewflag = STATIC;
   if (thetastr || phistr || cflag == DYNAMIC ||
       upxstr || upystr || upzstr || zoomstr || perspstr) viewflag = DYNAMIC;
 
   box_bounds();
   if (cflag == STATIC) box_center();
   if (viewflag == STATIC) view_params();
 
   // local data
 
   maxbufcopy = 0;
   chooseghost = NULL;
   bufcopy = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 DumpImage::~DumpImage()
 {
   delete image;
 
   delete [] diamtype;
   delete [] diamelement;
   delete [] colortype;
   delete [] colorelement;
   delete [] bdiamtype;
   delete [] bcolortype;
   memory->destroy(chooseghost);
   memory->destroy(bufcopy);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpImage::init_style()
 {
   if (multifile == 0)
     error->all(FLERR,"Dump image requires one snapshot per file");
   if (sort_flag) error->all(FLERR,"Dump image cannot perform sorting");
 
   DumpCustom::init_style();
 
   // check variables
 
   if (thetastr) {
     thetavar = input->variable->find(thetastr);
     if (thetavar < 0)
       error->all(FLERR,"Variable name for dump image theta does not exist");
     if (!input->variable->equalstyle(thetavar))
       error->all(FLERR,"Variable for dump image theta is invalid style");
   }
   if (phistr) {
     phivar = input->variable->find(phistr);
     if (phivar < 0)
       error->all(FLERR,"Variable name for dump image phi does not exist");
     if (!input->variable->equalstyle(phivar))
       error->all(FLERR,"Variable for dump image phi is invalid style");
   }
   if (cxstr) {
     cxvar = input->variable->find(cxstr);
     if (cxvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(cxvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (cystr) {
     cyvar = input->variable->find(cystr);
     if (cyvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(cyvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (czstr) {
     czvar = input->variable->find(czstr);
     if (czvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(czvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (upxstr) {
     upxvar = input->variable->find(upxstr);
     if (upxvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(upxvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (upystr) {
     upyvar = input->variable->find(upystr);
     if (upyvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(upyvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (upzstr) {
     upzvar = input->variable->find(upzstr);
     if (upzvar < 0)
       error->all(FLERR,"Variable name for dump image center does not exist");
     if (!input->variable->equalstyle(upzvar))
       error->all(FLERR,"Variable for dump image center is invalid style");
   }
   if (zoomstr) {
     zoomvar = input->variable->find(zoomstr);
     if (zoomvar < 0)
       error->all(FLERR,"Variable name for dump image zoom does not exist");
     if (!input->variable->equalstyle(zoomvar))
       error->all(FLERR,"Variable for dump image zoom is invalid style");
   }
   if (perspstr) {
     perspvar = input->variable->find(perspstr);
     if (perspvar < 0)
       error->all(FLERR,"Variable name for dump image persp does not exist");
     if (!input->variable->equalstyle(perspvar))
       error->all(FLERR,"Variable for dump image persp is invalid style");
   }
 
   // set up type -> element mapping
 
   if (atomflag && acolor == ELEMENT) {
     for (int i = 1; i <= ntypes; i++) {
       colorelement[i] = image->element2color(typenames[i]);
       if (colorelement[i] == NULL)
         error->all(FLERR,"Invalid dump image element name");
     }
   }
 
   if (atomflag && adiam == ELEMENT) {
     for (int i = 1; i <= ntypes; i++) {
       diamelement[i] = image->element2diam(typenames[i]);
       if (diamelement[i] == 0.0)
         error->all(FLERR,"Invalid dump image element name");
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpImage::write()
 {
   // open new file
 
   openfile();
 
   // reset box center and view parameters if dynamic
 
   box_bounds();
   if (cflag == DYNAMIC) box_center();
   if (viewflag == DYNAMIC) view_params();
 
   // nme = # of atoms this proc will contribute to dump
 
   nme = count();
 
   if (nme > maxbuf) {
     maxbuf = nme;
     memory->destroy(buf);
     memory->create(buf,maxbuf*size_one,"dump:buf");
   }
 
   // pack buf with color & diameter
 
   pack(NULL);
 
   // set minmax color range if using dynamic atom color map
 
   if (acolor == ATTRIBUTE && image->map_dynamic(0)) {
     double two[2],twoall[2];
     double lo = BIG;
     double hi = -BIG;
     int m = 0;
     for (int i = 0; i < nchoose; i++) {
       lo = MIN(lo,buf[m]);
       hi = MAX(hi,buf[m]);
       m += size_one;
     }
     two[0] = -lo;
     two[1] = hi;
     MPI_Allreduce(two,twoall,2,MPI_DOUBLE,MPI_MAX,world);
     int flag = image->map_minmax(0,-twoall[0],twoall[1]);
     if (flag) error->all(FLERR,"Invalid color map min/max values");
   }
 
   // create image on each proc, then merge them
 
   image->clear();
   create_image();
   image->merge();
 
   // write image file
 
   if (me == 0) {
     if (filetype == JPG) image->write_JPG(fp);
     else if (filetype == PNG) image->write_PNG(fp);
     else image->write_PPM(fp);
-    fclose(fp);
+    if (multifile) {
+      fclose(fp);
+      fp = NULL;
+    }
   }
 }
 
 /* ----------------------------------------------------------------------
    simulation box bounds
 ------------------------------------------------------------------------- */
 
 void DumpImage::box_bounds()
 {
   if (domain->triclinic == 0) {
     boxxlo = domain->boxlo[0];
     boxxhi = domain->boxhi[0];
     boxylo = domain->boxlo[1];
     boxyhi = domain->boxhi[1];
     boxzlo = domain->boxlo[2];
     boxzhi = domain->boxhi[2];
   } else {
     boxxlo = domain->boxlo_bound[0];
     boxxhi = domain->boxhi_bound[0];
     boxylo = domain->boxlo_bound[1];
     boxyhi = domain->boxhi_bound[1];
     boxzlo = domain->boxlo_bound[2];
     boxzhi = domain->boxhi_bound[2];
     boxxy = domain->xy;
     boxxz = domain->xz;
     boxyz = domain->yz;
   }
 }
 
 /* ----------------------------------------------------------------------
    reset view parameters
    called once from constructor if view is STATIC
    called every snapshot from write() if view is DYNAMIC
 ------------------------------------------------------------------------- */
 
 void DumpImage::box_center()
 {
   if (cxstr) cx = input->variable->compute_equal(cxvar);
   if (cystr) cy = input->variable->compute_equal(cyvar);
   if (czstr) cz = input->variable->compute_equal(czvar);
 
   image->xctr = boxxlo + cx*(boxxhi-boxxlo);
   image->yctr = boxylo + cy*(boxyhi-boxylo);
   image->zctr = boxzlo + cz*(boxzhi-boxzlo);
 }
 
 /* ----------------------------------------------------------------------
    reset view parameters in Image class
    called once from constructor if view is STATIC
    called every snapshot from write() if view is DYNAMIC
 ------------------------------------------------------------------------- */
 
 void DumpImage::view_params()
 {
   // view direction theta and phi
 
   if (thetastr) {
     double theta = input->variable->compute_equal(thetavar);
     if (theta < 0.0 || theta > 180.0)
       error->all(FLERR,"Invalid dump image theta value");
     theta *= MY_PI/180.0;
     image->theta = theta;
   }
 
   if (phistr) {
     double phi = input->variable->compute_equal(phivar);
     phi *= MY_PI/180.0;
     image->phi = phi;
   }
 
   // up vector
 
   if (upxstr) image->up[0] = input->variable->compute_equal(upxvar);
   if (upystr) image->up[1] = input->variable->compute_equal(upyvar);
   if (upzstr) image->up[2] = input->variable->compute_equal(upzvar);
 
   // zoom and perspective
 
   if (zoomstr) image->zoom = input->variable->compute_equal(zoomvar);
   if (image->zoom <= 0.0) error->all(FLERR,"Invalid dump image zoom value");
   if (perspstr) image->persp = input->variable->compute_equal(perspvar);
   if (image->persp < 0.0) error->all(FLERR,"Invalid dump image persp value");
 
   // remainder of view setup is internal to Image class
 
   image->view_params(boxxlo,boxxhi,boxylo,boxyhi,boxzlo,boxzhi);
 }
 
 /* ----------------------------------------------------------------------
    create image for atoms on this proc
    every pixel has depth
 ------------------------------------------------------------------------- */
 
 void DumpImage::create_image()
 {
   int i,j,m,itype,atom1,atom2;
   double diameter,delx,dely,delz;
   double *color,*color1,*color2;
   double xmid[3];
 
   // render my atoms
 
   if (atomflag) {
     double **x = atom->x;
 
     m = 0;
     for (i = 0; i < nchoose; i++) {
       j = clist[i];
 
       if (acolor == TYPE) {
         itype = static_cast<int> (buf[m]);
         color = colortype[itype];
       } else if (acolor == ELEMENT) {
         itype = static_cast<int> (buf[m]);
         color = colorelement[itype];
       } else if (acolor == ATTRIBUTE) {
         color = image->map_value2color(0,buf[m]);
       }
 
       if (adiam == NUMERIC) {
         diameter = adiamvalue;
       } else if (adiam == TYPE) {
         itype = static_cast<int> (buf[m+1]);
         diameter = diamtype[itype];
       } else if (adiam == ELEMENT) {
         itype = static_cast<int> (buf[m+1]);
         diameter = diamelement[itype];
       } else if (adiam == ATTRIBUTE) {
         diameter = buf[m+1];
       }
 
       image->draw_sphere(x[j],color,diameter);
       m += size_one;
     }
   }
 
   // render bonds for my atoms
   // both atoms in bond must be selected for bond to be rendered
   // if newton_bond is off, only render bond once
   // render bond in 2 pieces if crosses periodic boundary
   // if bond is deleted (type = 0), do not render
   // if bond is turned off (type < 0), still render
 
   if (bondflag) {
     double **x = atom->x;
     int *tag = atom->tag;
     int **bond_atom = atom->bond_atom;
     int **bond_type = atom->bond_type;
     int *num_bond = atom->num_bond;
     int *type = atom->type;
     int nlocal = atom->nlocal;
     int nall = atom->nlocal + atom->nghost;
     int newton_bond = force->newton_bond;
 
     // communicate choose flag for ghost atoms to know if they are selected
     // if bcolor/bdiam = ATOM, setup bufcopy to comm atom color/diam attributes
 
     if (nall > maxbufcopy) {
       maxbufcopy = atom->nmax;
       memory->destroy(chooseghost);
       memory->create(chooseghost,maxbufcopy,"dump:chooseghost");
       if (comm_forward == 3) {
         memory->destroy(bufcopy);
         memory->create(bufcopy,maxbufcopy,2,"dump:bufcopy");
       }
     }
 
     for (i = 0; i < nlocal; i++) chooseghost[i] = choose[i];
 
     if (comm_forward == 3) {
       for (i = 0; i < nlocal; i++) bufcopy[i][0] = bufcopy[i][1] = 0.0;
       m = 0;
       for (i = 0; i < nchoose; i++) {
         j = clist[i];
         bufcopy[j][0] = buf[m];
         bufcopy[j][1] = buf[m+1];
         m += size_one;
       }
     }
 
     comm->forward_comm_dump(this);
 
     for (i = 0; i < nchoose; i++) {
       atom1 = clist[i];
       for (m = 0; m < num_bond[atom1]; m++) {
         atom2 = atom->map(bond_atom[atom1][m]);
         if (atom2 < 0 || !chooseghost[atom2]) continue;
         if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
         if (bond_type[atom1][m] == 0) continue;
 
         if (bcolor == ATOM) {
           if (acolor == TYPE) {
             color1 = colortype[type[atom1]];
             color2 = colortype[type[atom2]];
           } else if (acolor == ELEMENT) {
             color1 = colorelement[type[atom1]];
             color2 = colorelement[type[atom2]];
           } else if (acolor == ATTRIBUTE) {
             color1 = image->map_value2color(0,bufcopy[atom1][0]);
             color2 = image->map_value2color(0,bufcopy[atom2][0]);
           }
         } else if (bcolor == TYPE) {
           itype = bond_type[atom1][m];
           if (itype < 0) itype = -itype;
           color = bcolortype[itype];
         }
 
         if (bdiam == NUMERIC) {
           diameter = bdiamvalue;
         } else if (bdiam == ATOM) {
           if (adiam == NUMERIC) {
             diameter = adiamvalue;
           } else if (adiam == TYPE) {
             diameter = MIN(diamtype[type[atom1]],diamtype[type[atom1]]);
           } else if (adiam == ELEMENT) {
             diameter = MIN(diamelement[type[atom1]],diamelement[type[atom1]]);
           } else if (adiam == ATTRIBUTE) {
             diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]);
           }
         } else if (bdiam == TYPE) {
           itype = bond_type[atom1][m];
           if (itype < 0) itype = -itype;
           diameter = bdiamtype[itype];
         }
 
         // draw cylinder in 2 pieces if bcolor = ATOM
         // or bond crosses periodic boundary
 
         delx = x[atom2][0] - x[atom1][0];
         dely = x[atom2][1] - x[atom1][1];
         delz = x[atom2][2] - x[atom1][2];
 
         if (bcolor == ATOM || domain->minimum_image_check(delx,dely,delz)) {
           domain->minimum_image(delx,dely,delz);
           xmid[0] = x[atom1][0] + 0.5*delx;
           xmid[1] = x[atom1][1] + 0.5*dely;
           xmid[2] = x[atom1][2] + 0.5*delz;
           if (bcolor == ATOM)
             image->draw_cylinder(x[atom1],xmid,color1,diameter,3);
           else image->draw_cylinder(x[atom1],xmid,color,diameter,3);
           xmid[0] = x[atom2][0] - 0.5*delx;
           xmid[1] = x[atom2][1] - 0.5*dely;
           xmid[2] = x[atom2][2] - 0.5*delz;
           if (bcolor == ATOM)
             image->draw_cylinder(xmid,x[atom2],color2,diameter,3);
           else image->draw_cylinder(xmid,x[atom2],color,diameter,3);
 
         } else image->draw_cylinder(x[atom1],x[atom2],color,diameter,3);
       }
     }
   }
 
   // render outline of simulation box, orthogonal or triclinic
 
   if (boxflag) {
     double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo);
     if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo);
     diameter *= boxdiam;
 
     double (*boxcorners)[3];
     double box[8][3];
     if (domain->triclinic == 0) {
       box[0][0] = boxxlo; box[0][1] = boxylo; box[0][2] = boxzlo;
       box[1][0] = boxxhi; box[1][1] = boxylo; box[1][2] = boxzlo;
       box[2][0] = boxxlo; box[2][1] = boxyhi; box[2][2] = boxzlo;
       box[3][0] = boxxhi; box[3][1] = boxyhi; box[3][2] = boxzlo;
       box[4][0] = boxxlo; box[4][1] = boxylo; box[4][2] = boxzhi;
       box[5][0] = boxxhi; box[5][1] = boxylo; box[5][2] = boxzhi;
       box[6][0] = boxxlo; box[6][1] = boxyhi; box[6][2] = boxzhi;
       box[7][0] = boxxhi; box[7][1] = boxyhi; box[7][2] = boxzhi;
       boxcorners = box;
     } else {
       domain->box_corners();
       boxcorners = domain->corners;
     }
 
     image->draw_box(boxcorners,diameter);
   }
 
   // render XYZ axes in red/green/blue
   // offset by 10% of box size and scale by axeslen
 
   if (axesflag) {
     double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo);
     if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo);
     diameter *= axesdiam;
 
     double (*boxcorners)[3];
     double axes[4][3];
     if (domain->triclinic == 0) {
       axes[0][0] = boxxlo; axes[0][1] = boxylo; axes[0][2] = boxzlo;
       axes[1][0] = boxxhi; axes[1][1] = boxylo; axes[1][2] = boxzlo;
       axes[2][0] = boxxlo; axes[2][1] = boxyhi; axes[2][2] = boxzlo;
       axes[3][0] = boxxlo; axes[3][1] = boxylo; axes[3][2] = boxzhi;
     } else {
       domain->box_corners();
       boxcorners = domain->corners;
       axes[0][0] = boxcorners[0][0];
       axes[0][1] = boxcorners[0][1];
       axes[0][2] = boxcorners[0][2];
       axes[1][0] = boxcorners[1][0];
       axes[1][1] = boxcorners[1][1];
       axes[1][2] = boxcorners[1][2];
       axes[2][0] = boxcorners[2][0];
       axes[2][1] = boxcorners[2][1];
       axes[2][2] = boxcorners[2][2];
       axes[3][0] = boxcorners[4][0];
       axes[3][1] = boxcorners[4][1];
       axes[3][2] = boxcorners[4][2];
     }
 
     double offset = MAX(boxxhi-boxxlo,boxyhi-boxylo);
     if (domain->dimension == 3) offset = MAX(offset,boxzhi-boxzlo);
     offset *= 0.1;
     axes[0][0] -= offset; axes[0][1] -= offset; axes[0][2] -= offset;
     axes[1][0] -= offset; axes[1][1] -= offset; axes[1][2] -= offset;
     axes[2][0] -= offset; axes[2][1] -= offset; axes[2][2] -= offset;
     axes[3][0] -= offset; axes[3][1] -= offset; axes[3][2] -= offset;
 
     axes[1][0] = axes[0][0] + axeslen*(axes[1][0]-axes[0][0]);
     axes[1][1] = axes[0][1] + axeslen*(axes[1][1]-axes[0][1]);
     axes[1][2] = axes[0][2] + axeslen*(axes[1][2]-axes[0][2]);
     axes[2][0] = axes[0][0] + axeslen*(axes[2][0]-axes[0][0]);
     axes[2][1] = axes[0][1] + axeslen*(axes[2][1]-axes[0][1]);
     axes[2][2] = axes[0][2] + axeslen*(axes[2][2]-axes[0][2]);
     axes[3][0] = axes[0][0] + axeslen*(axes[3][0]-axes[0][0]);
     axes[3][1] = axes[0][1] + axeslen*(axes[3][1]-axes[0][1]);
     axes[3][2] = axes[0][2] + axeslen*(axes[3][2]-axes[0][2]);
 
     image->draw_axes(axes,diameter);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int DumpImage::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
 {
   int i,j,m;
 
   m = 0;
 
   if (comm_forward == 1) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = chooseghost[j];
     }
   } else {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = chooseghost[j];
       buf[m++] = bufcopy[j][0];
       buf[m++] = bufcopy[j][1];
     }
   }
 
   return comm_forward;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DumpImage::unpack_comm(int n, int first, double *buf)
 {
   int i,m,last;
 
   m = 0;
   last = first + n;
 
   if (comm_forward == 1)
     for (i = first; i < last; i++) chooseghost[i] = static_cast<int> (buf[m++]);
   else {
     for (i = first; i < last; i++) {
       chooseghost[i] = static_cast<int> (buf[m++]);
       bufcopy[i][0] = buf[m++];
       bufcopy[i][1] = buf[m++];
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int DumpImage::modify_param(int narg, char **arg)
 {
   int n = DumpCustom::modify_param(narg,arg);
   if (n) return n;
 
   if (strcmp(arg[0],"acolor") == 0) {
     if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
     int nlo,nhi;
     force->bounds(arg[1],atom->ntypes,nlo,nhi);
 
     // ptrs = list of ncount colornames separated by '/'
 
     int ncount = 1;
     char *nextptr;
     char *ptr = arg[2];
     while (nextptr = strchr(ptr,'/')) {
       ptr = nextptr + 1;
       ncount++;
     }
     char **ptrs = new char*[ncount+1];
     ncount = 0;
     ptrs[ncount++] = strtok(arg[2],"/");
     while (ptrs[ncount++] = strtok(NULL,"/"));
     ncount--;
 
     // assign each of ncount colors in round-robin fashion to types
 
     int m = 0;
     for (int i = nlo; i <= nhi; i++) {
       colortype[i] = image->color2rgb(ptrs[m%ncount]);
       if (colortype[i] == NULL)
         error->all(FLERR,"Invalid color in dump_modify command");
       m++;
     }
 
     delete [] ptrs;
     return 3;
   }
 
   if (strcmp(arg[0],"adiam") == 0) {
     if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
     int nlo,nhi;
     force->bounds(arg[1],atom->ntypes,nlo,nhi);
     double diam = force->numeric(FLERR,arg[2]);
     if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command");
     for (int i = nlo; i <= nhi; i++) diamtype[i] = diam;
     return 3;
   }
 
   if (strcmp(arg[0],"amap") == 0) {
     if (narg < 6) error->all(FLERR,"Illegal dump_modify command");
     if (strlen(arg[3]) != 2) error->all(FLERR,"Illegal dump_modify command");
     int factor = 2;
     if (arg[3][0] == 's') factor = 1;
     int nentry = force->inumeric(FLERR,arg[5]);
     if (nentry < 1) error->all(FLERR,"Illegal dump_modify command");
     int n = 6 + factor*nentry;
     if (narg < n) error->all(FLERR,"Illegal dump_modify command");
     int flag = image->map_reset(0,n-1,&arg[1]);
     if (flag) error->all(FLERR,"Illegal dump_modify command");
     return n;
   }
 
   if (strcmp(arg[0],"bcolor") == 0) {
     if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
     if (atom->nbondtypes == 0)
       error->all(FLERR,"Dump modify bcolor not allowed with no bond types");
     int nlo,nhi;
     force->bounds(arg[1],atom->nbondtypes,nlo,nhi);
 
     // ptrs = list of ncount colornames separated by '/'
 
     int ncount = 1;
     char *nextptr;
     char *ptr = arg[2];
     while (nextptr = strchr(ptr,'/')) {
       ptr = nextptr + 1;
       ncount++;
     }
     char **ptrs = new char*[ncount+1];
     ncount = 0;
     ptrs[ncount++] = strtok(arg[2],"/");
     while (ptrs[ncount++] = strtok(NULL,"/"));
     ncount--;
 
     // assign each of ncount colors in round-robin fashion to types
 
     int m = 0;
     for (int i = nlo; i <= nhi; i++) {
       bcolortype[i] = image->color2rgb(ptrs[m%ncount]);
       if (bcolortype[i] == NULL)
         error->all(FLERR,"Invalid color in dump_modify command");
       m++;
     }
 
     delete [] ptrs;
     return 3;
   }
 
   if (strcmp(arg[0],"bdiam") == 0) {
     if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
     if (atom->nbondtypes == 0)
       error->all(FLERR,"Dump modify bdiam not allowed with no bond types");
     int nlo,nhi;
     force->bounds(arg[1],atom->ntypes,nlo,nhi);
     double diam = force->numeric(FLERR,arg[2]);
     if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command");
     for (int i = nlo; i <= nhi; i++) bdiamtype[i] = diam;
     return 3;
   }
 
   if (strcmp(arg[0],"backcolor") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     double *color = image->color2rgb(arg[1]);
     if (color == NULL) error->all(FLERR,"Invalid color in dump_modify command");
     image->background[0] = static_cast<int> (color[0]*255.0);
     image->background[1] = static_cast<int> (color[1]*255.0);
     image->background[2] = static_cast<int> (color[2]*255.0);
     return 2;
   }
 
   if (strcmp(arg[0],"boxcolor") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     image->boxcolor = image->color2rgb(arg[1]);
     if (image->boxcolor == NULL)
       error->all(FLERR,"Invalid color in dump_modify command");
     return 2;
   }
 
   if (strcmp(arg[0],"color") == 0) {
     if (narg < 5) error->all(FLERR,"Illegal dump_modify command");
     int flag = image->addcolor(arg[1],force->numeric(FLERR,arg[2]),force->numeric(FLERR,arg[3]),force->numeric(FLERR,arg[4]));
     if (flag) error->all(FLERR,"Illegal dump_modify command");
     return 5;
   }
 
   return 0;
 }
diff --git a/src/dump_image.h b/src/dump_image.h
index ab97d7f4a..9c42c584c 100644
--- a/src/dump_image.h
+++ b/src/dump_image.h
@@ -1,188 +1,190 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef DUMP_CLASS
 
 DumpStyle(image,DumpImage)
 
 #else
 
 #ifndef LMP_DUMP_IMAGE_H
 #define LMP_DUMP_IMAGE_H
 
 #include "dump_custom.h"
 
 namespace LAMMPS_NS {
 
 class DumpImage : public DumpCustom {
  public:
   DumpImage(class LAMMPS *, int, char**);
-  ~DumpImage();
+  virtual ~DumpImage();
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
 
- private:
+ protected:
   int filetype;
+  enum{PPM,JPG,PNG};
+
   int acolor,adiam;                // what determines color/diam of atoms
   double adiamvalue;               // atom diameter value
   int atomflag,bondflag;           // 0/1 for draw atoms,bonds
   int bcolor,bdiam;                // what determines color/diam of bonds
   double bdiamvalue;               // bond diameter value
   char *thetastr,*phistr;          // variables for view theta,phi
   int thetavar,phivar;             // index to theta,phi vars
   int cflag;                       // static/dynamic box center
   double cx,cy,cz;                 // fractional box center
   char *cxstr,*cystr,*czstr;       // variables for box center
   int cxvar,cyvar,czvar;           // index to box center vars
   char *upxstr,*upystr,*upzstr;    // view up vector variables
   int upxvar,upyvar,upzvar;        // index to up vector vars
   char *zoomstr,*perspstr;         // view zoom and perspective variables
   int zoomvar,perspvar;            // index to zoom,persp vars
   int boxflag,axesflag;            // 0/1 for draw box and axes
   double boxdiam,axeslen,axesdiam; // params for drawing box and axes
 
   int viewflag;                    // overall view is static or dynamic
 
   double *diamtype,*diamelement,*bdiamtype;         // per-type diameters
   double **colortype,**colorelement,**bcolortype;   // per-type colors
 
   class Image *image;              // class that renders each image
 
   int *chooseghost;                // extended choose array for comm
   double **bufcopy;                // buffer for communicating bond/atom info
   int maxbufcopy;
 
-  void init_style();
+  virtual void init_style();
   int modify_param(int, char **);
   void write();
 
   void box_center();
   void view_params();
   void box_bounds();
 
   void create_image();
 };
 
 }
 
 #endif
 #endif
 
 /* ERROR/WARNING messages:
 
 E: Invalid dump image filename
 
 The file produced by dump image cannot be binary and must
 be for a single processor.
 
 E: Support for writing images in JPEG format not included
 
 LAMMPS was not built with the -DLAMMPS_JPEG switch in the Makefile.
 
 E: Support for writing images in PNG format not included
 
 LAMMPS was not built with the -DLAMMPS_PNG switch in the Makefile.
 
 E: Illegal ... command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
 E: Dump image bond not allowed with no bond types
 
 Self-explanatory.
 
 E: Invalid dump image theta value
 
 Theta must be between 0.0 and 180.0 inclusive.
 
 E: Dump image persp option is not yet supported
 
 Self-explanatory.
 
 E: Dump image requires one snapshot per file
 
 Use a "*" in the filename.
 
 E: Dump image cannot perform sorting
 
 Self-explanatory.
 
 E: Variable name for dump image theta does not exist
 
 Self-explanatory.
 
 E: Variable for dump image theta is invalid style
 
 Must be an equal-style variable.
 
 E: Variable name for dump image phi does not exist
 
 Self-explanatory.
 
 E: Variable for dump image phi is invalid style
 
 Must be an equal-style variable.
 
 E: Variable name for dump image center does not exist
 
 Self-explanatory.
 
 E: Variable for dump image center is invalid style
 
 Must be an equal-style variable.
 
 E: Variable name for dump image zoom does not exist
 
 Self-explanatory.
 
 E: Variable for dump image zoom is invalid style
 
 Must be an equal-style variable.
 
 E: Variable name for dump image persp does not exist
 
 Self-explanatory.
 
 E: Variable for dump image persp is invalid style
 
 Must be an equal-style variable.
 
 E: Invalid dump image element name
 
 The specified element name was not in the standard list of elements.
 See the dump_modify doc page.
 
 E: Invalid dump image zoom value
 
 Zoom value must be > 0.0.
 
 E: Invalid dump image persp value
 
 Persp value must be >= 0.0.
 
 E: Invalid color in dump_modify command
 
 The specified color name was not in the list of recognized colors.
 See the dump_modify doc page.
 
 E: Dump modify bcolor not allowed with no bond types
 
 Self-explanatory.
 
 E: Dump modify bdiam not allowed with no bond types
 
 Self-explanatory.
 
 */
diff --git a/src/dump_movie.cpp b/src/dump_movie.cpp
new file mode 100644
index 000000000..0e6fc1952
--- /dev/null
+++ b/src/dump_movie.cpp
@@ -0,0 +1,59 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "dump_movie.h"
+#include "comm.h"
+#include "error.h"
+#include "memory.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+DumpMovie::DumpMovie(LAMMPS *lmp, int narg, char **arg) :
+  DumpImage(lmp, narg, arg)
+{
+  if (multiproc || compressed || multifile)
+    error->all(FLERR,"Invalid dump movie filename");
+
+  filetype = PPM;
+  bitrate = 2000;
+  framerate = 24;
+  fp = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpMovie::openfile()
+{
+  char moviecmd[1024];
+
+  if ((comm->me == 0) && (fp == NULL)) {
+    sprintf(moviecmd,"ffmpeg -y -r %d -f image2pipe -c:v ppm -i - "
+            "-r 24 -b:v %dk %s 2> /dev/null ", framerate, bitrate, filename);
+    fp = popen(moviecmd,"w");
+  }
+}
+/* ---------------------------------------------------------------------- */
+
+void DumpMovie::init_style()
+{
+  // initialize image style circumventing multifile check
+  multifile = 1;
+  DumpImage::init_style();
+  multifile = 0;
+}
diff --git a/src/dump_movie.h b/src/dump_movie.h
new file mode 100644
index 000000000..0ed5c577b
--- /dev/null
+++ b/src/dump_movie.h
@@ -0,0 +1,47 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef DUMP_CLASS
+
+DumpStyle(movie,DumpMovie)
+
+#else
+
+#ifndef LMP_DUMP_MOVIE_H
+#define LMP_DUMP_MOVIE_H
+
+#include "dump_image.h"
+
+namespace LAMMPS_NS {
+
+class DumpMovie : public DumpImage {
+ public:
+  DumpMovie(LAMMPS *, int, char**);
+
+  virtual void openfile();
+  virtual void init_style();
+
+ protected:
+  int bitrate;                  // bitrate of video file in kbps
+  int framerate;                // frame rate of animation
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+
+*/