diff --git a/src/MAKE/Makefile.g++ b/src/MAKE/Makefile.g++
index df0332892..90771bb1e 100755
--- a/src/MAKE/Makefile.g++
+++ b/src/MAKE/Makefile.g++
@@ -1,108 +1,111 @@
 # g++ = RedHat Linux box, g++4, gfortran, MPICH2, FFTW
 
 SHELL = /bin/sh
 
 # ---------------------------------------------------------------------
 # compiler/linker settings
 # specify flags and libraries needed for your compiler
 
 CC =		g++4
 CCFLAGS =	-g -O
 DEPFLAGS =	-M
 LINK =		g++4
 LINKFLAGS =	-g -O
 LIB =           
 ARCHIVE =	ar
 ARFLAGS =	-rc
 SIZE =		size
 
 # ---------------------------------------------------------------------
 # LAMMPS-specific settings
 # specify settings for LAMMPS features you will use
 
 # LAMMPS ifdef options, see doc/Section_start.html
 
 LMP_INC =	-DLAMMPS_GZIP -DLAMMPS_JPEG
 
 # MPI library, can be src/STUBS dummy lib
 # INC = path for mpi.h, MPI compiler settings
 # PATH = path for MPI library
 # LIB = name of MPI library
 
 MPI_INC =       -DMPICH_SKIP_MPICXX
 MPI_PATH = 
 MPI_LIB =	-lmpich -lpthread
 
 # FFT library, can be -DFFT_NONE if not using PPPM from KSPACE package
 # INC = -DFFT_FFTW, -DFFT_INTEL, -DFFT_NONE, etc, FFT compiler settings
 # PATH = path for FFT library
 # LIB = name of FFT library
 
 FFT_INC =       -DFFT_FFTW 
 FFT_PATH = 
 FFT_LIB =	-lfftw
 
 # JPEG library, only needed if -DLAMMPS_JPEG listed with LMP_INC
 # INC = path for jpeglib.h
 # PATH = path for JPEG library
 # LIB = name of JPEG library
 
 JPG_INC =       
 JPG_PATH = 	
 JPG_LIB =	/usr/local/lib/libjpeg.a
 
 # additional system settings needed by LAMMPS package libraries
 # these settings are IGNORED if the corresponding LAMMPS package
 #   (e.g. gpu, meam) is NOT included in the LAMMPS build
 # SYSINC = settings to compile with
 # SYSLIB = libraries to link with
 # SYSPATH = paths to libraries
 
 gpu_SYSINC =
 meam_SYSINC =
 reax_SYSINC =
 user-atc_SYSINC =
+user-awpmd_SYSINC =
 
 gpu_SYSLIB =       -lcudart -lcuda
 meam_SYSLIB =      -lgfortran
 reax_SYSLIB =      -lgfortran
 user-atc_SYSLIB =  -lblas -llapack
+user-awpmd_SYSLIB =  -lblas -llapack
 
 gpu_SYSPATH =      -L/usr/local/cuda/lib64
 meam_SYSPATH =
 reax_SYSPATH =
 user-atc_SYSPATH = 	
+user-awpmd_SYSPATH = 	
 
 # ---------------------------------------------------------------------
 # build rules and dependencies
 # no need to edit this section
 
 include	Makefile.package
 
 EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC)
 EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH)
 EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB)
 
 # Link target
 
 $(EXE):	$(OBJ)
 	$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
 	$(SIZE) $(EXE)
 
 # Library target
 
 lib:	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
 
 # Compilation rules
 
 %.o:%.cpp
 	$(CC) $(CCFLAGS) $(EXTRA_INC) -c $<
 
 %.d:%.cpp
 	$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
 
 # Individual dependencies
 
 DEPENDS = $(OBJ:.o=.d)
 include $(DEPENDS)
diff --git a/src/Makefile b/src/Makefile
index a1b626f66..5b3246591 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,186 +1,186 @@
 # LAMMPS multiple-machine Makefile
 
 SHELL = /bin/sh
 #.IGNORE:
 
 # Definitions
 
 ROOT =	lmp
 EXE =	$(ROOT)_$@
 SRC =	$(wildcard *.cpp)
 INC =	$(wildcard *.h)
 OBJ = 	$(SRC:.cpp=.o)
 
 # Package variables
 
 PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
 	  kspace manybody meam molecule opt peri poems reax replica \
 	  shock srd xtc
 
-PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-cuda \
-	   user-eff user-ewaldn user-imd user-reaxc user-smd
+PACKUSER = user-ackland user-atc user-awpmd user-cd-eam user-cg-cmm \
+	   user-cuda user-eff user-ewaldn user-imd user-reaxc user-smd
 
 PACKALL = $(PACKAGE) $(PACKUSER)
 
 PACKAGEUC = $(shell echo $(PACKAGE) | tr a-z A-Z)
 PACKUSERUC = $(shell echo $(PACKUSER) | tr a-z A-Z)
 
 YESDIR = $(shell echo $(@:yes-%=%) | tr a-z A-Z)
 NODIR  = $(shell echo $(@:no-%=%) | tr a-z A-Z)
 
 # List of all targets
 
 help:
 	@echo ''
 	@echo 'make clean-all           delete all object files'
 	@echo 'make clean-machine       delete object files for one machine'
 	@echo 'make tar                 lmp_src.tar.gz of src dir and packages'
 	@echo 'make makelib             update Makefile.lib for library build'
 	@echo 'make makelist            update Makefile.list used by old makes'
 	@echo ''
 	@echo 'make package             list available packages'
 	@echo 'make package-status      status of all packages'
 	@echo 'make yes-package         install a single package in src dir'
 	@echo 'make no-package          remove a single package from src dir'
 	@echo 'make yes-all             install all packages in src dir'
 	@echo 'make no-all              remove all packages from src dir'
 	@echo 'make yes-standard        install all standard packages'
 	@echo 'make no-standard         remove all standard packages'
 	@echo 'make yes-user            install all user packages'
 	@echo 'make no-user             remove all user packages'
 	@echo ''
 	@echo 'make package-update      replace src files with package files'
 	@echo 'make package-overwrite   replace package files with src files'
 	@echo ''
 	@echo 'make machine             build LAMMPS where machine is one of:'
 	@echo ''
 	@files="`ls MAKE/Makefile.*`"; \
 	  for file in $$files; do head -1 $$file; done
 	@echo ''
 
 # Build the code
 
 .DEFAULT:
 	@test -f MAKE/Makefile.$@
 	@if [ ! -d Obj_$@ ]; then mkdir Obj_$@; fi
 	@$(SHELL) Make.sh style
 	@cp -p *.cpp *.h Obj_$@
 	@cp MAKE/Makefile.$@ Obj_$@/Makefile
 	@if [ ! -e Makefile.package ]; then make package-regenerate; fi
 	@cp Makefile.package Obj_$@
 	@cd Obj_$@; \
 	$(MAKE) $(MFLAGS) "OBJ = $(OBJ)" "INC = $(INC)" "EXE = ../$(EXE)" ../$(EXE)
 	@if [ -d Obj_$@ ]; then cd Obj_$@; rm -f $(SRC) $(INC) Makefile*; fi
 
 # Remove machine-specific object files
 
 clean:
 	@echo 'make clean-all           delete all object files'
 	@echo 'make clean-machine       delete object files for one machine'
 
 clean-all:
 	rm -rf Obj_*
 
 clean-%:
 	rm -rf Obj_$(@:clean-%=%)
 
 # Create a tarball of src dir and packages
 
 tar:
 	@cd STUBS; make clean
 	@cd ..; tar cvzf src/$(ROOT)_src.tar.gz \
 	  src/Make* src/Package.sh src/MAKE src/*.cpp src/*.h src/STUBS \
 	  $(patsubst %,src/%,$(PACKAGEUC)) $(patsubst %,src/%,$(PACKUSERUC)) \
           --exclude=*/.svn
 	@cd STUBS; make
 	@echo "Created $(ROOT)_src.tar.gz"
 
 # Update Makefile.lib and Makefile.list
 
 makelib:
 	@$(SHELL) Make.sh style
 	@$(SHELL) Make.sh Makefile.lib
 
 makelist:
 	@$(SHELL) Make.sh style
 	@$(SHELL) Make.sh Makefile.list
 
 # Package management
 
 package:
 	@echo 'Standard packages:' $(PACKAGE)
 	@echo ''
 	@echo 'User-contributed packages:' $(PACKUSER)
 	@echo ''
 	@echo 'make package             list available packages'
 	@echo 'make package-status      status of all packages'
 	@echo 'make yes-package         install a single package in src dir'
 	@echo 'make no-package          remove a single package from src dir'
 	@echo 'make yes-all             install all packages in src dir'
 	@echo 'make no-all              remove all packages from src dir'
 	@echo 'make yes-standard        install all standard packages'
 	@echo 'make no-standard         remove all standard packages'
 	@echo 'make yes-user            install all user packages'
 	@echo 'make no-user             remove all user packages'
 	@echo ''
 	@echo 'make package-update      replace src files with package files'
 	@echo 'make package-overwrite   replace package files with src files'
 
 yes-all:
 	@for p in $(PACKALL); do $(MAKE) yes-$$p; done
 
 no-all:
 	@for p in $(PACKALL); do $(MAKE) no-$$p; done
 
 yes-standard:
 	@for p in $(PACKAGE); do $(MAKE) yes-$$p; done
 
 no-standard:
 	@for p in $(PACKAGE); do $(MAKE) no-$$p; done
 
 yes-user:
 	@for p in $(PACKUSER); do $(MAKE) yes-$$p; done
 
 no-user:
 	@for p in $(PACKUSER); do $(MAKE) no-$$p; done
 
 yes-%:
 	@if [ ! -e $(YESDIR) ]; then \
 	  echo "Package $(@:yes-%=%) does not exist"; \
 	else \
 	  echo "Installing package $(@:yes-%=%)"; \
 	  cd $(YESDIR); $(SHELL) Install.sh 1; cd ..; $(SHELL) Depend.sh 1; \
 	fi;
 
 no-%:
 	@if [ ! -e $(NODIR) ]; then \
 	  echo "Package $(@:no-%=%) does not exist"; \
 	else \
 	  echo "Uninstalling package $(@:no-%=%)"; \
 	  cd $(NODIR); $(SHELL) Install.sh 0; cd ..; $(SHELL) Depend.sh 0; \
         fi;
 
 # status = list differences between src and package files
 # update = replace src files with newer package files
 # overwrite = overwrite package files with newer src files
 # regenerate = regenerate Makefile.package from Makefile.package.empty
 
 package-status:
 	@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p status; done
 	@echo ''
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p status; done
 
 package-update:
 	@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p update; done
 	@echo ''
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p update; done
 
 package-overwrite:
 	@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p overwrite; done
 	@echo ''
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p overwrite; done
 
 package-regenerate:
 	@cp Makefile.package.empty Makefile.package
 	@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p regenerate; done
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p regenerate; done
diff --git a/src/atom.cpp b/src/atom.cpp
index 103fd5bc5..b9c9c3e70 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1,1659 +1,1663 @@
 /* ----------------------------------------------------------------------
    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 "mpi.h"
 #include "math.h"
 #include "stdio.h"
 #include "stdlib.h"
 #include "string.h"
 #include "limits.h"
 #include "atom.h"
 #include "style_atom.h"
 #include "atom_vec.h"
 #include "atom_vec_ellipsoid.h"
 #include "comm.h"
 #include "neighbor.h"
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
 #include "output.h"
 #include "thermo.h"
 #include "update.h"
 #include "domain.h"
 #include "group.h"
 #include "accelerator_cuda.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 #define DELTA 1
 #define DELTA_MEMSTR 1024
 #define EPSILON 1.0e-6
 #define CUDA_CHUNK 3000
 
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
 
 /* ---------------------------------------------------------------------- */
 
 Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 {
   natoms = 0;
   nlocal = nghost = nmax = 0;
   ntypes = 0;
   nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
   nbonds = nangles = ndihedrals = nimpropers = 0;
   bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
   extra_bond_per_atom = 0;
 
   firstgroupname = NULL;
   sortfreq = 1000;
   nextsort = 0;
   userbinsize = 0.0;
   maxbin = maxnext = 0;
   binhead = NULL;
   next = permute = NULL;
 
   // initialize atom arrays
   // customize by adding new array
 
   tag = type = mask = image = NULL;
   x = v = f = NULL;
 
   molecule = NULL;
   q = NULL;
   mu = NULL;
   omega = angmom = torque = NULL;
   radius = rmass = NULL;
   vfrac = s0 = NULL;
   x0 = NULL;
   ellipsoid = NULL;
   spin = NULL;
   eradius = ervel = erforce = NULL;
+  cs = csforce = vforce = ervelforce = NULL;
+  etag = NULL;
 
   maxspecial = 1;
   nspecial = NULL;
   special = NULL;
 
   num_bond = NULL;
   bond_type = bond_atom = NULL;
 
   num_angle = NULL;
   angle_type = angle_atom1 = angle_atom2 = angle_atom3 = NULL;
 
   num_dihedral = NULL;
   dihedral_type = dihedral_atom1 = dihedral_atom2 = NULL;
   dihedral_atom3 = dihedral_atom4 = NULL;
 
   num_improper = NULL;
   improper_type = improper_atom1 = improper_atom2 = NULL;
   improper_atom3 = improper_atom4 = NULL;
 
   // initialize atom style and array existence flags
   // customize by adding new flag
 
   sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
+  wavepacket_flag = 0;
 
   molecule_flag = q_flag = mu_flag = 0;
   rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
   vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
+  cs_flag = csforce_flag = vforce_flag = ervelforce_flag= etag_flag = 0;
 
   // ntype-length arrays
 
   mass = NULL;
   mass_setflag = NULL;
 
   // callback lists & extra restart info
 
   nextra_grow = nextra_restart = 0;
   extra_grow = extra_restart = NULL;
   nextra_grow_max = nextra_restart_max = 0;
   nextra_store = 0;
   extra = NULL;
 
   // default mapping values and hash table primes
 
   tag_enable = 1;
   map_style = 0;
   map_tag_max = 0;
   map_nhash = 0;
   
   nprimes = 38;
   primes = new int[nprimes];
   int plist[] = {5041,10007,20011,30011,40009,50021,60013,70001,80021,
 		 90001,100003,110017,120011,130003,140009,150001,160001,
 		 170003,180001,190027,200003,210011,220009,230003,240007,
 		 250007,260003,270001,280001,290011,300007,310019,320009,
 		 330017,340007,350003,362881,3628801};
   for (int i = 0; i < nprimes; i++) primes[i] = plist[i];
 
   // default atom style = atomic
 
   atom_style = NULL;
   avec = NULL;
   create_avec("atomic",0,NULL);
 }
 
 /* ---------------------------------------------------------------------- */
 
 Atom::~Atom()
 {
   delete [] atom_style;
   delete avec;
 
   delete [] firstgroupname;
   memory->destroy(binhead);
   memory->destroy(next);
   memory->destroy(permute);
 
   // delete atom arrays
   // customize by adding new array
 
   memory->destroy(tag);
   memory->destroy(type);
   memory->destroy(mask);
   memory->destroy(image);
   memory->destroy(x);
   memory->destroy(v);
   memory->destroy(f);
 
   memory->destroy(q);
   memory->destroy(mu);
   memory->destroy(omega);
   memory->destroy(angmom);
   memory->destroy(torque);
   memory->destroy(radius);
   memory->destroy(rmass);
   memory->destroy(vfrac);
   memory->destroy(s0);
   memory->destroy(x0);
   memory->destroy(ellipsoid);
   memory->destroy(spin);
   memory->destroy(eradius);
   memory->destroy(ervel);
   memory->destroy(erforce);
 
   memory->destroy(molecule);
 
   memory->destroy(nspecial);
   memory->destroy(special);
 
   memory->destroy(num_bond);
   memory->destroy(bond_type);
   memory->destroy(bond_atom);
 
   memory->destroy(num_angle);
   memory->destroy(angle_type);
   memory->destroy(angle_atom1);
   memory->destroy(angle_atom2);
   memory->destroy(angle_atom3);
 
   memory->destroy(num_dihedral);
   memory->destroy(dihedral_type);
   memory->destroy(dihedral_atom1);
   memory->destroy(dihedral_atom2);
   memory->destroy(dihedral_atom3);
   memory->destroy(dihedral_atom4);
 
   memory->destroy(num_improper);
   memory->destroy(improper_type);
   memory->destroy(improper_atom1);
   memory->destroy(improper_atom2);
   memory->destroy(improper_atom3);
   memory->destroy(improper_atom4);
 
   // delete per-type arrays
 
   delete [] mass;
   delete [] mass_setflag;
 
   // delete extra arrays
 
   memory->destroy(extra_grow);
   memory->destroy(extra_restart);
   memory->destroy(extra);
 
   // delete mapping data structures
 
   map_delete();
   delete [] primes;
 }
 
 /* ----------------------------------------------------------------------
    copy modify settings from old Atom class to current Atom class
 ------------------------------------------------------------------------- */
 
 void Atom::settings(Atom *old)
 {
   map_style = old->map_style;
 }
 
 /* ----------------------------------------------------------------------
    create an AtomVec style
    called from input script, restart file, replicate
 ------------------------------------------------------------------------- */
 
 void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
 {
   delete [] atom_style;
   if (avec) delete avec;
 
   // unset atom style and array existence flags
   // may have been set by old avec
   // customize by adding new flag
 
   sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
 
   molecule_flag = q_flag = mu_flag = 0;
   rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
   vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
 
   int sflag;
   avec = new_avec(style,narg,arg,suffix,sflag);
 
   if (sflag) {
     char estyle[256];
     sprintf(estyle,"%s/%s",style,suffix);
     int n = strlen(estyle) + 1;
     atom_style = new char[n];
     strcpy(atom_style,estyle);
   } else {
     int n = strlen(style) + 1;
     atom_style = new char[n];
     strcpy(atom_style,style);
   }
 
   // if molecular system, default is to have array map
 
   molecular = avec->molecular;
   if (map_style == 0 && molecular) map_style = 1;
 }
 
 /* ----------------------------------------------------------------------
    generate an AtomVec class, first with suffix appended
 ------------------------------------------------------------------------- */
 
 AtomVec *Atom::new_avec(const char *style, int narg, char **arg,
 			char *suffix, int &sflag)
 {
   if (suffix && lmp->suffix_enable) {
     sflag = 1;
     char estyle[256];
     sprintf(estyle,"%s/%s",style,suffix);
 
     if (0) return NULL;
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
     else if (strcmp(estyle,#key) == 0) return new Class(lmp,narg,arg);
 #include "style_atom.h"
 #undef AtomStyle
 #undef ATOM_CLASS
 
   }
 
   sflag = 0;
 
   if (0) return NULL;
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
   else if (strcmp(style,#key) == 0) return new Class(lmp,narg,arg);
 #include "style_atom.h"
 #undef ATOM_CLASS
 
   else error->all("Invalid atom style");
 
   return NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Atom::init()
 {
   // delete extra array since it doesn't persist past first run
 
   if (nextra_store) {
     memory->destroy(extra);
     extra = NULL;
     nextra_store = 0;
   }
 
   // check arrays that are atom type in length
 
   check_mass();
 
   // setup of firstgroup
 
   if (firstgroupname) {
     firstgroup = group->find(firstgroupname);
     if (firstgroup < 0)
       error->all("Could not find atom_modify first group ID");
   } else firstgroup = -1;
 
   // init AtomVec
 
   avec->init();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Atom::setup()
 {
   // setup bins for sorting
   // cannot do this in init() because uses neighbor cutoff
 
   if (sortfreq > 0) setup_sort_bins();
 }
 
 /* ----------------------------------------------------------------------
    return ptr to AtomVec class if matches style or to matching hybrid sub-class
    return NULL if no match
 ------------------------------------------------------------------------- */
 
 AtomVec *Atom::style_match(const char *style)
 {
   if (strcmp(atom_style,style) == 0) return avec;
   else if (strcmp(atom_style,"hybrid") == 0) {
     AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec;
     for (int i = 0; i < avec_hybrid->nstyles; i++)
       if (strcmp(avec_hybrid->keywords[i],style) == 0)
 	return avec_hybrid->styles[i];
   }
   return NULL;
 }
 
 /* ----------------------------------------------------------------------
    modify parameters of the atom style
    some options can only be invoked before simulation box is defined
    first and sort options cannot be used together
 ------------------------------------------------------------------------- */
 
 void Atom::modify_params(int narg, char **arg)
 {
   if (narg == 0) error->all("Illegal atom_modify command");
 
   int iarg = 0;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"map") == 0) {
       if (iarg+2 > narg) error->all("Illegal atom_modify command");
       if (strcmp(arg[iarg+1],"array") == 0) map_style = 1;
       else if (strcmp(arg[iarg+1],"hash") == 0) map_style = 2;
       else error->all("Illegal atom_modify command");
       if (domain->box_exist) 
 	error->all("Atom_modify map command after simulation box is defined");
       iarg += 2;
     } else if (strcmp(arg[iarg],"first") == 0) {
       if (iarg+2 > narg) error->all("Illegal atom_modify command");
       if (strcmp(arg[iarg+1],"all") == 0) {
 	delete [] firstgroupname;
 	firstgroupname = NULL;
       } else {
 	int n = strlen(arg[iarg+1]) + 1;
 	firstgroupname = new char[n];
 	strcpy(firstgroupname,arg[iarg+1]);
 	sortfreq = 0;
       }
       iarg += 2;
     } else if (strcmp(arg[iarg],"sort") == 0) {
       if (iarg+3 > narg) error->all("Illegal atom_modify command");
       sortfreq = atoi(arg[iarg+1]);
       userbinsize = atof(arg[iarg+2]);
       if (sortfreq < 0 || userbinsize < 0.0)
 	error->all("Illegal atom_modify command");
       if (sortfreq >= 0 && firstgroupname) 
 	error->all("Atom_modify sort and first options "
 		   "cannot be used together");
       iarg += 3;
     } else error->all("Illegal atom_modify command");
   }
 }
 
 /* ----------------------------------------------------------------------
    allocate and initialize array or hash table for global -> local map
    set map_tag_max = largest atom ID (may be larger than natoms)
    for array option:
      array length = 1 to largest tag of any atom
      set entire array to -1 as initial values
    for hash option:
      map_nhash = length of hash table
      map_nbucket = # of hash buckets, prime larger than map_nhash
        so buckets will only be filled with 0 or 1 atoms on average
 ------------------------------------------------------------------------- */
 
 void Atom::map_init()
 {
   map_delete();
 
   if (tag_enable == 0)
     error->all("Cannot create an atom map unless atoms have IDs");
 
   int max = 0;
   for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
   MPI_Allreduce(&max,&map_tag_max,1,MPI_INT,MPI_MAX,world);
 
   if (map_style == 1) {
     memory->create(map_array,map_tag_max+1,"atom:map_array");
     for (int i = 0; i <= map_tag_max; i++) map_array[i] = -1;
 
   } else {
 
     // map_nhash = max of atoms/proc or total atoms, times 2, at least 1000
 
     int nper = static_cast<int> (natoms/comm->nprocs);
     map_nhash = MAX(nper,nmax);
     if (map_nhash > natoms) map_nhash = static_cast<int> (natoms);
     if (comm->nprocs > 1) map_nhash *= 2;
     map_nhash = MAX(map_nhash,1000);
 
     // map_nbucket = prime just larger than map_nhash
 
     int n = map_nhash/10000;
     n = MIN(n,nprimes-1);
     map_nbucket = primes[n];
     if (map_nbucket < map_nhash && n < nprimes-1) map_nbucket = primes[n+1];
 
     // set all buckets to empty
     // set hash to map_nhash in length
     // put all hash entries in free list and point them to each other
 
     map_bucket = new int[map_nbucket];
     for (int i = 0; i < map_nbucket; i++) map_bucket[i] = -1;
 
     map_hash = new HashElem[map_nhash];
     map_nused = 0;
     map_free = 0;
     for (int i = 0; i < map_nhash; i++) map_hash[i].next = i+1;
     map_hash[map_nhash-1].next = -1;
   }
 }
 
 /* ----------------------------------------------------------------------
    clear global -> local map for all of my own and ghost atoms
    for hash table option:
      global ID may not be in table if image atom was already cleared
 ------------------------------------------------------------------------- */
 
 void Atom::map_clear()
 {
   if (map_style == 1) {
     int nall = nlocal + nghost;
     for (int i = 0; i < nall; i++) map_array[tag[i]] = -1;
 
   } else {
     int previous,global,ibucket,index;
     int nall = nlocal + nghost;
     for (int i = 0; i < nall; i++) {
 
       // search for key
       // if don't find it, done
 
       previous = -1;
       global = tag[i];
       ibucket = global % map_nbucket;
       index = map_bucket[ibucket];
       while (index > -1) {
 	if (map_hash[index].global == global) break;
 	previous = index;
 	index = map_hash[index].next;
       }
       if (index == -1) continue;
 
       // delete the hash entry and add it to free list
       // special logic if entry is 1st in the bucket
     
       if (previous == -1) map_bucket[ibucket] = map_hash[index].next;
       else map_hash[previous].next = map_hash[index].next;
     
       map_hash[index].next = map_free;
       map_free = index;
       map_nused--;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    set global -> local map for all of my own and ghost atoms
    loop in reverse order so that nearby images take precedence over far ones
      and owned atoms take precedence over images
    this enables valid lookups of bond topology atoms 
    for hash table option:
      if hash table too small, re-init
      global ID may already be in table if image atom was set
 ------------------------------------------------------------------------- */
 
 void Atom::map_set()
 {
   if (map_style == 1) {
     int nall = nlocal + nghost;
     for (int i = nall-1; i >= 0 ; i--) map_array[tag[i]] = i;
 
   } else {
     int previous,global,ibucket,index;
     int nall = nlocal + nghost;
     if (nall > map_nhash) map_init();
 
     for (int i = nall-1; i >= 0 ; i--) {
 
       // search for key
       // if found it, just overwrite local value with index
 
       previous = -1;
       global = tag[i];
       ibucket = global % map_nbucket;
       index = map_bucket[ibucket];
       while (index > -1) {
 	if (map_hash[index].global == global) break;
 	previous = index;
 	index = map_hash[index].next;
       }
       if (index > -1) {
 	map_hash[index].local = i;
 	continue;
       }
 
       // take one entry from free list
       // add the new global/local pair as entry at end of bucket list
       // special logic if this entry is 1st in bucket
 
       index = map_free;
       map_free = map_hash[map_free].next;
       if (previous == -1) map_bucket[ibucket] = index;
       else map_hash[previous].next = index;
       map_hash[index].global = global;
       map_hash[index].local = i;
       map_hash[index].next = -1;
       map_nused++;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    set global to local map for one atom
    for hash table option:
      global ID may already be in table if atom was already set
 ------------------------------------------------------------------------- */
 
 void Atom::map_one(int global, int local)
 {
   if (map_style == 1) map_array[global] = local;
 
   else {
     // search for key
     // if found it, just overwrite local value with index
     
     int previous = -1;
     int ibucket = global % map_nbucket;
     int index = map_bucket[ibucket];
     while (index > -1) {
       if (map_hash[index].global == global) break;
       previous = index;
       index = map_hash[index].next;
     }
     if (index > -1) {
       map_hash[index].local = local;
       return;
     }
     
     // take one entry from free list
     // add the new global/local pair as entry at end of bucket list
     // special logic if this entry is 1st in bucket
     
     index = map_free;
     map_free = map_hash[map_free].next;
     if (previous == -1) map_bucket[ibucket] = index;
     else map_hash[previous].next = index;
     map_hash[index].global = global;
     map_hash[index].local = local;
     map_hash[index].next = -1;
     map_nused++;
   }
 }
 
 /* ----------------------------------------------------------------------
    free the array or hash table for global to local mapping
 ------------------------------------------------------------------------- */
 
 void Atom::map_delete()
 {
   if (map_style == 1) {
     if (map_tag_max) memory->destroy(map_array);
   } else {
     if (map_nhash) {
       delete [] map_bucket;
       delete [] map_hash;
     }
     map_nhash = 0;
   }
   map_tag_max = 0;
 }
 
 /* ----------------------------------------------------------------------
    lookup global ID in hash table, return local index
 ------------------------------------------------------------------------- */
 
 int Atom::map_find_hash(int global)
 {
   int local = -1;
   int index = map_bucket[global % map_nbucket];
   while (index > -1) {
     if (map_hash[index].global == global) {
       local = map_hash[index].local;
       break;
     }
     index = map_hash[index].next;
   }
   return local;
 }
 
 /* ----------------------------------------------------------------------
    add unique tags to any atoms with tag = 0
    new tags are grouped by proc and start after max current tag
    called after creating new atoms 
 ------------------------------------------------------------------------- */
 
 void Atom::tag_extend()
 {
   // maxtag_all = max tag for all atoms
 
   int maxtag = 0;
   for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
   int maxtag_all;
   MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world);
 
   // notag = # of atoms I own with no tag (tag = 0)
   // notag_sum = # of total atoms on procs <= me with no tag
 
   int notag = 0;
   for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++;
   int notag_sum;
   MPI_Scan(&notag,&notag_sum,1,MPI_INT,MPI_SUM,world);
 
   // itag = 1st new tag that my untagged atoms should use
 
   int itag = maxtag_all + notag_sum - notag + 1;
   for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++;
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs span range from 1 to Natoms
    return 0 if mintag != 1 or maxtag != Natoms
    return 1 if OK
    doesn't actually check if all tag values are used
 ------------------------------------------------------------------------- */
 
 int Atom::tag_consecutive()
 {
   int idmin = MAXTAGINT;
   int idmax = 0;
   
   for (int i = 0; i < nlocal; i++) {
     idmin = MIN(idmin,tag[i]);
     idmax = MAX(idmax,tag[i]);
   }
   int idminall,idmaxall;
   MPI_Allreduce(&idmin,&idminall,1,MPI_INT,MPI_MIN,world);
   MPI_Allreduce(&idmax,&idmaxall,1,MPI_INT,MPI_MAX,world);
 
   if (idminall != 1 || idmaxall != static_cast<int> (natoms)) return 0;
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    count and return words in a single line
    make copy of line before using strtok so as not to change line
    trim anything from '#' onward
 ------------------------------------------------------------------------- */
 
 int Atom::count_words(const char *line)
 {
   int n = strlen(line) + 1;
   char *copy;
   memory->create(copy,n,"atom:copy");
   strcpy(copy,line);
 
   char *ptr;
   if (ptr = strchr(copy,'#')) *ptr = '\0';
 
   if (strtok(copy," \t\n\r\f") == NULL) {
     memory->destroy(copy);
     return 0;
   }
   n = 1;
   while (strtok(NULL," \t\n\r\f")) n++;
 
   memory->destroy(copy);
   return n;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from Atom section of data file
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_atoms(int n, char *buf)
 {
   int m,imagedata,xptr,iptr;
   double xdata[3],lamda[3],sublo[3],subhi[3];
   double *coord;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3)
     error->all("Incorrect atom format in data file");
 
   char **values = new char*[nwords];
 
   // set bounds for my proc
   // if periodic and I am lo/hi proc, adjust bounds by EPSILON
   // insures all data atoms will be owned even with round-off
 
   int triclinic = domain->triclinic;
   if (triclinic == 0) {
     sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
     sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
     sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
   } else {
     sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
     sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
     sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
   }
 
   if (domain->xperiodic) {
     if (comm->myloc[0] == 0) sublo[0] -= EPSILON;
     if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += EPSILON;
   }
   if (domain->yperiodic) {
     if (comm->myloc[1] == 0) sublo[1] -= EPSILON;
     if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += EPSILON;
   }
   if (domain->zperiodic) {
     if (comm->myloc[2] == 0) sublo[2] -= EPSILON;
     if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += EPSILON;
   }
 
   // xptr = which word in line starts xyz coords
   // iptr = which word in line starts ix,iy,iz image flags
 
   xptr = avec->xcol_data - 1;
   int imageflag = 0;
   if (nwords > avec->size_data_atom) imageflag = 1;
   if (imageflag) iptr = nwords - 3;
 
   // loop over lines of atom data
   // tokenize the line into values
   // extract xyz coords and image flags
   // remap atom into simulation box
   // if atom is in my sub-domain, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     if (values[0] == NULL) error->all("Incorrect atom format in data file");
     for (m = 1; m < nwords; m++) {
       values[m] = strtok(NULL," \t\n\r\f");
       if (values[m] == NULL) error->all("Incorrect atom format in data file");
     }
 
     if (imageflag)
       imagedata = ((atoi(values[iptr+2]) + 512 & 1023) << 20) |
 	((atoi(values[iptr+1]) + 512 & 1023) << 10) |
 	(atoi(values[iptr]) + 512 & 1023);
     else imagedata = (512 << 20) | (512 << 10) | 512;
 
     xdata[0] = atof(values[xptr]);
     xdata[1] = atof(values[xptr+1]);
     xdata[2] = atof(values[xptr+2]);
     domain->remap(xdata,imagedata);
     if (triclinic) {
       domain->x2lamda(xdata,lamda);
       coord = lamda;
     } else coord = xdata;
 
     if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
 	coord[1] >= sublo[1] && coord[1] < subhi[1] &&
 	coord[2] >= sublo[2] && coord[2] < subhi[2])
       avec->data_atom(xdata,imagedata,values);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from Velocity section of data file
    check that atom IDs are > 0 and <= map_tag_max
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_vels(int n, char *buf)
 {
   int j,m,tagdata;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec->size_data_vel)
     error->all("Incorrect velocity format in data file");
 
   char **values = new char*[nwords];
 
   // loop over lines of atom velocities
   // tokenize the line into values
   // if I own atom tag, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     for (j = 1; j < nwords; j++)
       values[j] = strtok(NULL," \t\n\r\f");
 
     tagdata = atoi(values[0]);
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one("Invalid atom ID in Velocities section of data file");
     if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from atom-style specific section of data file
    check that atom IDs are > 0 and <= map_tag_max
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus)
 {
   int j,m,tagdata;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec_bonus->size_data_bonus)
     error->all("Incorrect bonus data format in data file");
 
   char **values = new char*[nwords];
 
   // loop over lines of bonus atom data
   // tokenize the line into values
   // if I own atom tag, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     for (j = 1; j < nwords; j++)
       values[j] = strtok(NULL," \t\n\r\f");
 
     tagdata = atoi(values[0]);
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one("Invalid atom ID in Bonus section of data file");
 
     // ok to call child's data_atom_bonus() method thru parent avec_bonus,
     // since data_bonus() was called with child ptr, and method is virtual
 
     if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_bonds(int n, char *buf)
 {
   int m,tmp,itype,atom1,atom2;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d %d %d",&tmp,&itype,&atom1,&atom2);
     if (atom1 <= 0 || atom1 > map_tag_max || 
 	atom2 <= 0 || atom2 > map_tag_max)
       error->one("Invalid atom ID in Bonds section of data file");
     if (itype <= 0 || itype > nbondtypes)
       error->one("Invalid bond type in Bonds section of data file");
     if ((m = map(atom1)) >= 0) {
       bond_type[m][num_bond[m]] = itype;
       bond_atom[m][num_bond[m]] = atom2;
       num_bond[m]++;
     }
     if (newton_bond == 0) {
       if ((m = map(atom2)) >= 0) {
 	bond_type[m][num_bond[m]] = itype;
 	bond_atom[m][num_bond[m]] = atom1;
 	num_bond[m]++;
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_angles(int n, char *buf)
 {
   int m,tmp,itype,atom1,atom2,atom3;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3);
     if (atom1 <= 0 || atom1 > map_tag_max || 
 	atom2 <= 0 || atom2 > map_tag_max || 
 	atom3 <= 0 || atom3 > map_tag_max)
       error->one("Invalid atom ID in Angles section of data file");
     if (itype <= 0 || itype > nangletypes)
       error->one("Invalid angle type in Angles section of data file");
     if ((m = map(atom2)) >= 0) {
       angle_type[m][num_angle[m]] = itype;
       angle_atom1[m][num_angle[m]] = atom1;
       angle_atom2[m][num_angle[m]] = atom2;
       angle_atom3[m][num_angle[m]] = atom3;
       num_angle[m]++;
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
 	angle_type[m][num_angle[m]] = itype;
 	angle_atom1[m][num_angle[m]] = atom1;
 	angle_atom2[m][num_angle[m]] = atom2;
 	angle_atom3[m][num_angle[m]] = atom3;
 	num_angle[m]++;
       }
       if ((m = map(atom3)) >= 0) {
 	angle_type[m][num_angle[m]] = itype;
 	angle_atom1[m][num_angle[m]] = atom1;
 	angle_atom2[m][num_angle[m]] = atom2;
 	angle_atom3[m][num_angle[m]] = atom3;
 	num_angle[m]++;
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_dihedrals(int n, char *buf)
 {
   int m,tmp,itype,atom1,atom2,atom3,atom4;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3,&atom4);
     if (atom1 <= 0 || atom1 > map_tag_max || 
 	atom2 <= 0 || atom2 > map_tag_max || 
 	atom3 <= 0 || atom3 > map_tag_max || 
 	atom4 <= 0 || atom4 > map_tag_max)
       error->one("Invalid atom ID in Dihedrals section of data file");
     if (itype <= 0 || itype > ndihedraltypes)
       error->one("Invalid dihedral type in Dihedrals section of data file");
     if ((m = map(atom2)) >= 0) {
       dihedral_type[m][num_dihedral[m]] = itype;
       dihedral_atom1[m][num_dihedral[m]] = atom1;
       dihedral_atom2[m][num_dihedral[m]] = atom2;
       dihedral_atom3[m][num_dihedral[m]] = atom3;
       dihedral_atom4[m][num_dihedral[m]] = atom4;
       num_dihedral[m]++;
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
 	dihedral_type[m][num_dihedral[m]] = itype;
 	dihedral_atom1[m][num_dihedral[m]] = atom1;
 	dihedral_atom2[m][num_dihedral[m]] = atom2;
 	dihedral_atom3[m][num_dihedral[m]] = atom3;
 	dihedral_atom4[m][num_dihedral[m]] = atom4;
 	num_dihedral[m]++;
       }
       if ((m = map(atom3)) >= 0) {
 	dihedral_type[m][num_dihedral[m]] = itype;
 	dihedral_atom1[m][num_dihedral[m]] = atom1;
 	dihedral_atom2[m][num_dihedral[m]] = atom2;
 	dihedral_atom3[m][num_dihedral[m]] = atom3;
 	dihedral_atom4[m][num_dihedral[m]] = atom4;
 	num_dihedral[m]++;
       }
       if ((m = map(atom4)) >= 0) {
 	dihedral_type[m][num_dihedral[m]] = itype;
 	dihedral_atom1[m][num_dihedral[m]] = atom1;
 	dihedral_atom2[m][num_dihedral[m]] = atom2;
 	dihedral_atom3[m][num_dihedral[m]] = atom3;
 	dihedral_atom4[m][num_dihedral[m]] = atom4;
 	num_dihedral[m]++;
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_impropers(int n, char *buf)
 {
   int m,tmp,itype,atom1,atom2,atom3,atom4;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3,&atom4);
     if (atom1 <= 0 || atom1 > map_tag_max || 
 	atom2 <= 0 || atom2 > map_tag_max || 
 	atom3 <= 0 || atom3 > map_tag_max || 
 	atom4 <= 0 || atom4 > map_tag_max)
       error->one("Invalid atom ID in Impropers section of data file");
     if (itype <= 0 || itype > nimpropertypes)
       error->one("Invalid improper type in Impropers section of data file");
     if ((m = map(atom2)) >= 0) {
       improper_type[m][num_improper[m]] = itype;
       improper_atom1[m][num_improper[m]] = atom1;
       improper_atom2[m][num_improper[m]] = atom2;
       improper_atom3[m][num_improper[m]] = atom3;
       improper_atom4[m][num_improper[m]] = atom4;
       num_improper[m]++;
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
 	improper_type[m][num_improper[m]] = itype;
 	improper_atom1[m][num_improper[m]] = atom1;
 	improper_atom2[m][num_improper[m]] = atom2;
 	improper_atom3[m][num_improper[m]] = atom3;
 	improper_atom4[m][num_improper[m]] = atom4;
 	num_improper[m]++;
       }
       if ((m = map(atom3)) >= 0) {
 	improper_type[m][num_improper[m]] = itype;
 	improper_atom1[m][num_improper[m]] = atom1;
 	improper_atom2[m][num_improper[m]] = atom2;
 	improper_atom3[m][num_improper[m]] = atom3;
 	improper_atom4[m][num_improper[m]] = atom4;
 	num_improper[m]++;
       }
       if ((m = map(atom4)) >= 0) {
 	improper_type[m][num_improper[m]] = itype;
 	improper_atom1[m][num_improper[m]] = atom1;
 	improper_atom2[m][num_improper[m]] = atom2;
 	improper_atom3[m][num_improper[m]] = atom3;
 	improper_atom4[m][num_improper[m]] = atom4;
 	num_improper[m]++;
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    allocate arrays of length ntypes
    only done after ntypes is set
 ------------------------------------------------------------------------- */
 
 void Atom::allocate_type_arrays()
 {
   if (avec->mass_type) {
     mass = new double[ntypes+1];
     mass_setflag = new int[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
   }
 }
 
 /* ----------------------------------------------------------------------
    set a mass and flag it as set
    called from reading of data file
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(const char *str)
 {
   if (mass == NULL) error->all("Cannot set mass for this atom style");
 
   int itype;
   double mass_one;
   int n = sscanf(str,"%d %lg",&itype,&mass_one);
   if (n != 2) error->all("Invalid mass line in data file");
 
   if (itype < 1 || itype > ntypes) error->all("Invalid type for mass set");
 
   mass[itype] = mass_one;
   mass_setflag[itype] = 1;
 
   if (mass[itype] <= 0.0) error->all("Invalid mass value");
 }
 
 /* ----------------------------------------------------------------------
    set a mass and flag it as set
    called from EAM pair routine
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(int itype, double value)
 {
   if (mass == NULL) error->all("Cannot set mass for this atom style");
   if (itype < 1 || itype > ntypes) error->all("Invalid type for mass set");
 
   mass[itype] = value;
   mass_setflag[itype] = 1;
 
   if (mass[itype] <= 0.0) error->all("Invalid mass value");
 }
 
 /* ----------------------------------------------------------------------
    set one or more masses and flag them as set
    called from reading of input script
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(int narg, char **arg)
 {
   if (mass == NULL) error->all("Cannot set mass for this atom style");
 
   int lo,hi;
   force->bounds(arg[0],ntypes,lo,hi);
   if (lo < 1 || hi > ntypes) error->all("Invalid type for mass set");
 
   for (int itype = lo; itype <= hi; itype++) {
     mass[itype] = atof(arg[1]);
     mass_setflag[itype] = 1;
 
     if (mass[itype] <= 0.0) error->all("Invalid mass value");
   }
 }
 
 /* ----------------------------------------------------------------------
    set all masses as read in from restart file
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(double *values)
 {
   for (int itype = 1; itype <= ntypes; itype++) {
     mass[itype] = values[itype];
     mass_setflag[itype] = 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    check that all masses have been set
 ------------------------------------------------------------------------- */
 
 void Atom::check_mass()
 {
   if (mass == NULL) return;
   for (int itype = 1; itype <= ntypes; itype++)
     if (mass_setflag[itype] == 0) error->all("All masses are not set");
 }
 
 /* ----------------------------------------------------------------------
    check that radii of all particles of itype are the same
    return 1 if true, else return 0
    also return the radius value for that type
 ------------------------------------------------------------------------- */
 
 int Atom::radius_consistency(int itype, double &rad)
 {
   double value = -1.0;
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
     if (type[i] != itype) continue;
     if (value < 0.0) value = radius[i];
     else if (value != radius[i]) flag = 1;
   }
 
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall) return 0;
 
   MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    check that shape of all particles of itype are the same
    return 1 if true, else return 0
    also return the 3 shape params for itype
 ------------------------------------------------------------------------- */
 
 int Atom::shape_consistency(int itype,
 			    double &shapex, double &shapey, double &shapez)
 {
   double zero[3] = {0.0, 0.0, 0.0};
   double one[3] = {-1.0, -1.0, -1.0};
   double *shape;
 
   AtomVecEllipsoid *avec_ellipsoid =
     (AtomVecEllipsoid *) style_match("ellipsoid");
   AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
 
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
     if (type[i] != itype) continue;
     if (ellipsoid[i] < 0) shape = zero;
     else shape = bonus[ellipsoid[i]].shape;
 
     if (one[0] < 0.0) {
       one[0] = shape[0];
       one[1] = shape[1];
       one[2] = shape[2];
     } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
       flag = 1;
   }
 
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall) return 0;
 
   double oneall[3];
   MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
   shapex = oneall[0];
   shapey = oneall[1];
   shapez = oneall[2];
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    reorder owned atoms so those in firstgroup appear first
    called by comm->exchange() if atom_modify first group is set
    only owned atoms exist at this point, no ghost atoms
 ------------------------------------------------------------------------- */
 
 void Atom::first_reorder()
 {
   // insure there is one extra atom location at end of arrays for swaps
 
   if (nlocal == nmax) avec->grow(0);
 
   // loop over owned atoms
   // nfirst = index of first atom not in firstgroup
   // when find firstgroup atom out of place, swap it with atom nfirst
 
   int bitmask = group->bitmask[firstgroup];
   nfirst = 0;
   while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & bitmask && i > nfirst) {
       avec->copy(i,nlocal,0);
       avec->copy(nfirst,i,0);
       avec->copy(nlocal,nfirst,0);
       while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform spatial sort of atoms within my sub-domain
    always called between comm->exchange() and comm->borders()
    don't have to worry about clearing/setting atom->map since done in comm
 ------------------------------------------------------------------------- */
 
 void Atom::sort()
 {
   int i,m,n,ix,iy,iz,ibin,empty;
 
   // set next timestep for sorting to take place
 
   nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
 
   // download data from GPU if necessary
 
   if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->downloadAll();
 
   // re-setup sort bins if needed
 
   if (domain->box_change) setup_sort_bins();
   if (nbins == 1) return;
 
   // reallocate per-atom vectors if needed
 
   if (nlocal > maxnext) {
     memory->destroy(next);
     memory->destroy(permute);
     maxnext = atom->nmax;
     memory->create(next,maxnext,"atom:next");
     memory->create(permute,maxnext,"atom:permute");
   }
 
   // insure there is one extra atom location at end of arrays for swaps
 
   if (nlocal == nmax) avec->grow(0);
 
   // bin atoms in reverse order so linked list will be in forward order
 
   for (i = 0; i < nbins; i++) binhead[i] = -1;
 
   for (i = nlocal-1; i >= 0; i--) {
     ix = static_cast<int> ((x[i][0]-bboxlo[0])*bininvx);
     iy = static_cast<int> ((x[i][1]-bboxlo[1])*bininvy);
     iz = static_cast<int> ((x[i][2]-bboxlo[2])*bininvz);
     ix = MAX(ix,0);
     iy = MAX(iy,0);
     iz = MAX(iz,0);
     ix = MIN(ix,nbinx-1);
     iy = MIN(iy,nbiny-1);
     iz = MIN(iz,nbinz-1);
     ibin = iz*nbiny*nbinx + iy*nbinx + ix;
     next[i] = binhead[ibin];
     binhead[ibin] = i;
   }
 
   // permute = desired permutation of atoms
   // permute[I] = J means Ith new atom will be Jth old atom
 
   n = 0;
   for (m = 0; m < nbins; m++) {
     i = binhead[m];
     while (i >= 0) {
       permute[n++] = i;
       i = next[i];
     }
   }
 
   // current = current permutation, just reuse next vector
   // current[I] = J means Ith current atom is Jth old atom
 
   int *current = next;
   for (i = 0; i < nlocal; i++) current[i] = i;
 
   // reorder local atom list, when done, current = permute
   // perform "in place" using copy() to extra atom location at end of list
   // inner while loop processes one cycle of the permutation
   // copy before inner-loop moves an atom to end of atom list
   // copy after inner-loop moves atom at end of list back into list
   // empty = location in atom list that is currently empty
 
   for (i = 0; i < nlocal; i++) {
     if (current[i] == permute[i]) continue;
     avec->copy(i,nlocal,0);
     empty = i;
     while (permute[empty] != i) {
       avec->copy(permute[empty],empty,0);
       empty = current[empty] = permute[empty];
     }      
     avec->copy(nlocal,empty,0);
     current[empty] = permute[empty];
   }
 
   // upload data back to GPU if necessary
 
   if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->uploadAll();
 
   // sanity check that current = permute
 
   //int flag = 0;
   //for (i = 0; i < nlocal; i++)
   //  if (current[i] != permute[i]) flag = 1;
   //int flagall;
   //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   //if (flagall) error->all("Atom sort did not operate correctly");
 }
 
 /* ----------------------------------------------------------------------
    setup bins for spatial sorting of atoms
 ------------------------------------------------------------------------- */
 
 void Atom::setup_sort_bins()
 {
   // binsize:
   // user setting if explicitly set
   // 1/2 of neighbor cutoff for non-CUDA
   // CUDA_CHUNK atoms/proc for CUDA
   // check if neighbor cutoff = 0.0
 
   double binsize;
   if (userbinsize > 0.0) binsize = userbinsize;
   else if (!lmp->cuda) binsize = 0.5 * neighbor->cutneighmax;
   else {
     if (domain->dimension == 3) {
       double vol = (domain->boxhi[0]-domain->boxlo[0]) * 
 	(domain->boxhi[1]-domain->boxlo[1]) * 
 	(domain->boxhi[2]-domain->boxlo[2]);
       binsize = pow(1.0*CUDA_CHUNK/natoms*vol,1.0/3.0);
     } else {
       double area = (domain->boxhi[0]-domain->boxlo[0]) * 
 	(domain->boxhi[1]-domain->boxlo[1]);
       binsize = pow(1.0*CUDA_CHUNK/natoms*area,1.0/2.0);
     }
   }
   if (binsize == 0.0) error->all("Atom sorting has bin size = 0.0");
 
   double bininv = 1.0/binsize;
 
   // nbin xyz = local bins
   // bbox lo/hi = bounding box of my sub-domain
  
   if (domain->triclinic)
     domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
   else {
     bboxlo[0] = domain->sublo[0];
     bboxlo[1] = domain->sublo[1];
     bboxlo[2] = domain->sublo[2];
     bboxhi[0] = domain->subhi[0];
     bboxhi[1] = domain->subhi[1];
     bboxhi[2] = domain->subhi[2];
   }
   
   nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) * bininv);
   nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) * bininv);
   nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) * bininv);
   if (domain->dimension == 2) nbinz = 1;
   
   if (nbinx == 0) nbinx = 1;
   if (nbiny == 0) nbiny = 1;
   if (nbinz == 0) nbinz = 1;
 
   bininvx = nbinx / (bboxhi[0]-bboxlo[0]);
   bininvy = nbiny / (bboxhi[1]-bboxlo[1]);
   bininvz = nbinz / (bboxhi[2]-bboxlo[2]);
 
   if (1.0*nbinx*nbiny*nbinz > INT_MAX) 
     error->one("Too many atom sorting bins");
 
   nbins = nbinx*nbiny*nbinz;
 
   // reallocate per-bin memory if needed
 
   if (nbins > maxbin) {
     memory->destroy(binhead);
     maxbin = nbins;
     memory->create(binhead,maxbin,"atom:binhead");
   }
 }
 
 /* ----------------------------------------------------------------------
    register a callback to a fix so it can manage atom-based arrays
    happens when fix is created
    flag = 0 for grow, 1 for restart
 ------------------------------------------------------------------------- */
 
 void Atom::add_callback(int flag)
 {
   int ifix;
 
   // find the fix
   // if find NULL ptr:
   //   it's this one, since it is being replaced and has just been deleted
   //   at this point in re-creation
   // if don't find NULL ptr:
   //   i is set to nfix = new one currently being added at end of list
 
   for (ifix = 0; ifix < modify->nfix; ifix++)
     if (modify->fix[ifix] == NULL) break;
 
   // add callback to lists, reallocating if necessary
 
   if (flag == 0) {
     if (nextra_grow == nextra_grow_max) {
       nextra_grow_max += DELTA;
       memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow");
     }
     extra_grow[nextra_grow] = ifix;
     nextra_grow++;
   } else if (flag == 1) {
     if (nextra_restart == nextra_restart_max) {
       nextra_restart_max += DELTA;
       memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart");
     }
     extra_restart[nextra_restart] = ifix;
     nextra_restart++;
   }
 }
 
 /* ----------------------------------------------------------------------
    unregister a callback to a fix
    happens when fix is deleted, called by its destructor
    flag = 0 for grow, 1 for restart
 ------------------------------------------------------------------------- */
 
 void Atom::delete_callback(const char *id, int flag)
 {
   int ifix;
   for (ifix = 0; ifix < modify->nfix; ifix++)
     if (strcmp(id,modify->fix[ifix]->id) == 0) break;
 
   // compact the list of callbacks
 
   if (flag == 0) {
     int match;
     for (match = 0; match < nextra_grow; match++)
       if (extra_grow[match] == ifix) break;
     for (int i = match; i < nextra_grow-1; i++)
       extra_grow[i] = extra_grow[i+1];
     nextra_grow--;
 
   } else if (flag == 1) {
     int match;
     for (match = 0; match < nextra_grow; match++)
       if (extra_restart[match] == ifix) break;
     for (int i = ifix; i < nextra_restart-1; i++)
       extra_restart[i] = extra_restart[i+1];
     nextra_restart--;
   }
 }
 
 /* ----------------------------------------------------------------------
    decrement ptrs in callback lists to fixes beyond the deleted ifix
    happens after fix is deleted
 ------------------------------------------------------------------------- */
 
 void Atom::update_callback(int ifix)
 {
   for (int i = 0; i < nextra_grow; i++)
     if (extra_grow[i] > ifix) extra_grow[i]--;
   for (int i = 0; i < nextra_restart; i++)
     if (extra_restart[i] > ifix) extra_restart[i]--;
 }
 
 /* ----------------------------------------------------------------------
    return a pointer to a named internal variable
    if don't recognize name, return NULL
    customize by adding names
 ------------------------------------------------------------------------- */
 
 void *Atom::extract(char *name)
 {
   if (strcmp(name,"id") == 0) return (void *) tag;
   if (strcmp(name,"type") == 0) return (void *) type;
   if (strcmp(name,"x") == 0) return (void *) x;
   if (strcmp(name,"v") == 0) return (void *) v;
   if (strcmp(name,"f") == 0) return (void *) f;
   if (strcmp(name,"mass") == 0) return (void *) mass;
   if (strcmp(name,"rmass") == 0) return (void *) rmass;
 
   return NULL;
 }
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
    call to avec tallies per-atom vectors
    add in global to local mapping storage
 ------------------------------------------------------------------------- */
 
 bigint Atom::memory_usage()
 {
   memlength = DELTA_MEMSTR;
   memory->create(memstr,memlength,"atom:memstr");
   memstr[0] = '\0';
   bigint bytes = avec->memory_usage();
   memory->destroy(memstr);
 
   if (map_style == 1)
     bytes += memory->usage(map_array,map_tag_max+1);
   else if (map_style == 2) {
     bytes += map_nbucket*sizeof(int);
     bytes += map_nhash*sizeof(HashElem);
   }
   if (maxnext) {
     bytes += memory->usage(next,maxnext);
     bytes += memory->usage(permute,maxnext);
   }
 
   return bytes;
 }
 
 /* ----------------------------------------------------------------------
    accumulate per-atom vec names in memstr, padded by spaces
    return 1 if padded str is not already in memlist, else 0
 ------------------------------------------------------------------------- */
 
 int Atom::memcheck(const char *str)
 {
   int n = strlen(str) + 3;
   char *padded = new char[n];
   strcpy(padded," ");
   strcat(padded,str);
   strcat(padded," ");
   
   if (strstr(memstr,padded)) {
     delete [] padded;
     return 0;
   }
 
   if (strlen(memstr) + n >= memlength) {
     memlength += DELTA_MEMSTR;
     memory->grow(memstr,memlength,"atom:memstr");
   }
 
   strcat(memstr,padded);
   delete [] padded;
   return 1;
 }
diff --git a/src/atom.h b/src/atom.h
index e4972b27e..9b56510ca 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -1,221 +1,225 @@
 /* ----------------------------------------------------------------------
    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.
 ------------------------------------------------------------------------- */
 
 #ifndef LMP_ATOM_H
 #define LMP_ATOM_H
 
 #include "pointers.h"
 
 namespace LAMMPS_NS {
 
 class Atom : protected Pointers {
  public:
   char *atom_style;
   class AtomVec *avec;
 
   // atom counts
 
   bigint natoms;                // total # of atoms in system, could be 0
   int nlocal,nghost;            // # of owned and ghost atoms on this proc
   int nmax;                     // max # of owned+ghost in arrays on this proc
   int tag_enable;               // 0/1 if atom ID tags are defined
   int molecular;                // 0 = atomic, 1 = molecular system
 
   bigint nbonds,nangles,ndihedrals,nimpropers;
   int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
   int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
   int extra_bond_per_atom;
 
   int firstgroup;               // store atoms in this group first, -1 if unset
   int nfirst;                   // # of atoms in first group on this proc
   char *firstgroupname;         // group-ID to store first, NULL if unset
 
   // per-atom arrays
   // customize by adding new array
 
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
 
   int *molecule;
   double *q,**mu;
   double **omega,**angmom,**torque;
   double *radius,*rmass,*vfrac,*s0;
   double **x0;
   int *ellipsoid;
   int *spin;
-  double *eradius,*ervel,*erforce;
+  double *eradius,*ervel,*erforce,*ervelforce;
+  double *cs,*csforce,*vforce;
+  int *etag;
 
   int **nspecial;               // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
   int **special;                // IDs of 1-2,1-3,1-4 neighs of each atom
   int maxspecial;               // special[nlocal][maxspecial]
 
   int *num_bond;
   int **bond_type;
   int **bond_atom;
 
   int *num_angle;
   int **angle_type;
   int **angle_atom1,**angle_atom2,**angle_atom3;
 
   int *num_dihedral;
   int **dihedral_type;
   int **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
 
   int *num_improper;
   int **improper_type;
   int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
 
   // atom style and per-atom array existence flags
   // customize by adding new flag
 
   int sphere_flag,ellipsoid_flag,peri_flag,electron_flag;
+  int wavepacket_flag;
 
   int molecule_flag,q_flag,mu_flag;
   int rmass_flag,radius_flag,omega_flag,torque_flag,angmom_flag;
   int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
+  int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag;
 
   // extra peratom info in restart file destined for fix & diag 
 
   double **extra;
 
   // per-type arrays
 
   double *mass;
   int *mass_setflag;
 
   // callback ptrs for atom arrays managed by fix classes
 
   int nextra_grow,nextra_restart;             // # of callbacks of each type
   int *extra_grow,*extra_restart;             // index of fix to callback to
   int nextra_grow_max,nextra_restart_max;     // size of callback lists
   int nextra_store;
 
   int map_style;                  // default or user-specified style of map
                                   // 0 = none, 1 = array, 2 = hash
 
   // spatial sorting of atoms
 
   int sortfreq;             // sort atoms every this many steps, 0 = off
   bigint nextsort;          // next timestep to sort on
 
   // functions
 
   Atom(class LAMMPS *);
   ~Atom();
 
   void settings(class Atom *);
   void create_avec(const char *, int, char **, char *suffix = NULL);
   class AtomVec *new_avec(const char *, int, char **, char *, int &);
   void init();
   void setup();
 
   class AtomVec *style_match(const char *);
   void modify_params(int, char **);
   void tag_extend();
   int tag_consecutive();
 
   int parse_data(const char *);
   int count_words(const char *);
 
   void data_atoms(int, char *);
   void data_vels(int, char *);
   void data_bonus(int, char *, class AtomVec *);
 
   void data_bonds(int, char *);
   void data_angles(int, char *);
   void data_dihedrals(int, char *);
   void data_impropers(int, char *);
 
   void allocate_type_arrays();
   void set_mass(const char *);
   void set_mass(int, double);
   void set_mass(int, char **);
   void set_mass(double *);
   void check_mass();
 
   int radius_consistency(int, double &);
   int shape_consistency(int, double &, double &, double &);
 
   void first_reorder();
   void sort();
 
   void add_callback(int);
   void delete_callback(const char *, int);
   void update_callback(int);
 
   void *extract(char *);
 
   inline int* get_map_array() {return map_array;};
   inline int get_map_size() {return map_tag_max+1;};
 
   bigint memory_usage();
   int memcheck(const char *);
 
   // functions for global to local ID mapping
   // map lookup function inlined for efficiency
   
   inline int map(int global) {
     if (map_style == 1) return map_array[global];
     else return map_find_hash(global);
   };
 
   void map_init();
   void map_clear();
   void map_set();
   void map_one(int, int);
   void map_delete();
   int map_find_hash(int);
 
  private:
 
   // global to local ID mapping
 
   int map_tag_max;
   int *map_array;
 
   struct HashElem {
     int global;                   // key to search on = global ID
     int local;                    // value associated with key = local index
     int next;                     // next entry in this bucket, -1 if last
   };
   int map_nhash;                  // # of entries hash table can hold
   int map_nused;                  // # of actual entries in hash table
   int map_free;                   // ptr to 1st unused entry in hash table
   int map_nbucket;                // # of hash buckets
   int *map_bucket;                // ptr to 1st entry in each bucket
   HashElem *map_hash;             // hash table
   int *primes;                    // table of prime #s for hashing
   int nprimes;                    // # of primes
 
   // spatial sorting of atoms
 
   int nbins;                      // # of sorting bins
   int nbinx,nbiny,nbinz;          // bins in each dimension
   int maxbin;                     // max # of bins
   int maxnext;                    // max size of next,permute
   int *binhead;                   // 1st atom in each bin
   int *next;                      // next atom in bin
   int *permute;                   // permutation vector
   double userbinsize;             // requested sort bin size
   double bininvx,bininvy,bininvz; // inverse actual bin sizes
   double bboxlo[3],bboxhi[3];     // bounding box of my sub-domain
 
   int memlength;                  // allocated size of memstr
   char *memstr;                   // string of array names already counted
 
   void setup_sort_bins();
 };
 
 }
 
 #endif
diff --git a/src/compute_slice.cpp b/src/compute_slice.cpp
index 15096b8c8..51e7747e3 100644
--- a/src/compute_slice.cpp
+++ b/src/compute_slice.cpp
@@ -1,312 +1,312 @@
 /* ----------------------------------------------------------------------
    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 "stdlib.h"
 #include "string.h"
 #include "compute_slice.h"
 #include "update.h"
 #include "modify.h"
 #include "fix.h"
 #include "group.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 enum{COMPUTE,FIX};
 
 #define INVOKED_VECTOR 2
 #define INVOKED_ARRAY 4
 
 /* ---------------------------------------------------------------------- */
 
 ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
   if (narg < 7) error->all("Illegal compute slice command");
 
   MPI_Comm_rank(world,&me);
 
   nstart = atoi(arg[3]);
   nstop = atoi(arg[4]);
   nskip = atoi(arg[5]);
 
   if (nstart < 1 || nstop < nstart || nskip < 1)
     error->all("Illegal compute slice command");
 
   // parse remaining values until one isn't recognized
 
   which = new int[narg-6];
   argindex = new int[narg-6];
   ids = new char*[narg-6];
   value2index = new int[narg-6];
   nvalues = 0;
 
   for (int iarg = 6; iarg < narg; iarg++) {
     if (strncmp(arg[iarg],"c_",2) == 0 || 
 	strncmp(arg[iarg],"f_",2) == 0) {
       if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
       else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
 
       int n = strlen(arg[iarg]);
       char *suffix = new char[n];
       strcpy(suffix,&arg[iarg][2]);
 
       char *ptr = strchr(suffix,'[');
       if (ptr) {
 	if (suffix[strlen(suffix)-1] != ']')
 	  error->all("Illegal compute slice command");
 	argindex[nvalues] = atoi(ptr+1);
 	*ptr = '\0';
       } else argindex[nvalues] = 0;
 
       n = strlen(suffix) + 1;
       ids[nvalues] = new char[n];
       strcpy(ids[nvalues],suffix);
       nvalues++;
       delete [] suffix;
 
     } else error->all("Illegal compute slice command");
   }
 
   // setup and error check
 
   for (int i = 0; i < nvalues; i++) {
     if (which[i] == COMPUTE) {
       int icompute = modify->find_compute(ids[i]);
       if (icompute < 0)
 	error->all("Compute ID for compute slice does not exist");
       if (modify->compute[icompute]->vector_flag) {
 	if (argindex[i])
 	  error->all("Compute slice compute does not calculate a global array");
 	if (nstop > modify->compute[icompute]->size_vector)
 	  error->all("Compute slice compute vector is accessed out-of-range");
       } else if (modify->compute[icompute]->array_flag) {
 	if (argindex[i] == 0)
 	  error->all("Compute slice compute does not calculate a global vector");
 	if (argindex[i] > modify->compute[icompute]->size_array_cols)
 	  error->all("Compute slice compute array is accessed out-of-range");
 	if (nstop > modify->compute[icompute]->size_array_rows)
 	  error->all("Compute slice compute array is accessed out-of-range");
       } else error->all("Compute slice compute does not calculate "
 			"global vector or array");
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
 	error->all("Fix ID for compute slice does not exist");
       if (modify->fix[ifix]->vector_flag) {
 	if (argindex[i])
 	  error->all("Compute slice fix does not calculate a global array");
 	if (nstop > modify->fix[ifix]->size_vector)
 	  error->all("Compute slice fix vector is accessed out-of-range");
       } else if (modify->fix[ifix]->array_flag) {
 	if (argindex[i] == 0) 
 	  error->all("Compute slice fix does not calculate a global vector");
 	if (argindex[i] > modify->fix[ifix]->size_array_cols)
 	  error->all("Compute slice fix array is accessed out-of-range");
 	if (nstop > modify->fix[ifix]->size_array_rows)
 	  error->all("Compute slice fix array is accessed out-of-range");
       } else error->all("Compute slice fix does not calculate "
 			"global vector or array");
     }
   }
 
   // this compute produces either a vector or array
   // for vector, set intensive/extensive to mirror input values
   // for array, set intensive if all input values are intensive, else extensive
 
   vector = NULL;
   array = NULL;
   extlist = NULL;
 
   if (nvalues == 1) {
     vector_flag = 1;
     size_vector = (nstop-nstart) / nskip;
     memory->create(vector,size_vector,"slice:vector");
 
     if (which[0] == COMPUTE) {
       int icompute = modify->find_compute(ids[0]);
       if (argindex[0] == 0) {
 	extvector = modify->compute[icompute]->extvector;
 	if (modify->compute[icompute]->extvector == -1) {
 	  extlist = new int[size_vector];
 	  int j = 0;
 	  for (int i = nstart; i < nstop; i += nskip)
 	    extlist[j++] = modify->compute[icompute]->extlist[i-1];
 	}
       } else extvector = modify->compute[icompute]->extarray;
-    } else if (which[0] = FIX) {
+    } else if (which[0] == FIX) {
       int ifix = modify->find_fix(ids[0]);
       if (argindex[0] == 0) {
 	extvector = modify->fix[ifix]->extvector;
 	if (modify->fix[ifix]->extvector == -1) {
 	  extlist = new int[size_vector];
 	  int j = 0;
 	  for (int i = nstart; i < nstop; i += nskip)
 	    extlist[j++] = modify->fix[ifix]->extlist[i-1];
 	}
       } else extvector = modify->fix[ifix]->extarray;
     }
 
   } else {
     array_flag = 1;
     size_array_rows = (nstop-nstart) / nskip;
     size_array_cols = nvalues;
     memory->create(array,size_array_rows,size_array_cols,"slice:array");
 
     extarray = 0;
     for (int i = 0; i < nvalues; i++) {
       if (which[i] == COMPUTE) {
 	int icompute = modify->find_compute(ids[i]);
 	if (argindex[i] == 0) {
 	  if (modify->compute[icompute]->extvector == 1) extarray = 1;
 	  if (modify->compute[icompute]->extvector == -1) {
 	    for (int j = 0; j < modify->compute[icompute]->size_vector; j++)
 	      if (modify->compute[icompute]->extlist[j]) extarray = 1;
 	  }
 	} else {
 	  if (modify->compute[icompute]->extarray) extarray = 1;
 	}
       } else if (which[i] == FIX) {
 	int ifix = modify->find_fix(ids[i]);
 	if (argindex[i] == 0) {
 	  if (modify->fix[ifix]->extvector == 1) extarray = 1;
 	  if (modify->fix[ifix]->extvector == -1) {
 	    for (int j = 0; j < modify->fix[ifix]->size_vector; j++)
 	      if (modify->fix[ifix]->extlist[j]) extarray = 1;
 	  }
 	} else {
 	  if (modify->fix[ifix]->extarray) extarray = 1;
 	}
       }
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeSlice::~ComputeSlice()
 {
   delete [] which;
   delete [] argindex;
   for (int m = 0; m < nvalues; m++) delete [] ids[m];
   delete [] ids;
   delete [] value2index;
 
   memory->destroy(vector);
   memory->destroy(array);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeSlice::init()
 {
   // set indices and check validity of all computes,fixes
 
   for (int m = 0; m < nvalues; m++) {
     if (which[m] == COMPUTE) {
       int icompute = modify->find_compute(ids[m]);
       if (icompute < 0)
 	error->all("Compute ID for compute slice does not exist");
       value2index[m] = icompute;
     } else if (which[m] == FIX) {
       int ifix = modify->find_fix(ids[m]);
       if (ifix < 0) 
 	error->all("Fix ID for compute slice does not exist");
       value2index[m] = ifix;
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeSlice::compute_vector()
 {
   invoked_vector = update->ntimestep;
 
   extract_one(0,vector,1);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeSlice::compute_array()
 {
   invoked_array = update->ntimestep;
 
   for (int m = 0; m < nvalues; m++)
     extract_one(0,&array[m][0],nvalues);
 }
 
 /* ----------------------------------------------------------------------
    calculate sliced value for one input M and return it in vec
    vec may be array so that returned values are with stride
 ------------------------------------------------------------------------- */
 
 void ComputeSlice::extract_one(int m, double *vec, int stride)
 {
   int i,j;
 
   // invoke the appropriate compute if needed
 
   if (which[m] == COMPUTE) {
     Compute *compute = modify->compute[value2index[m]];
 
     if (argindex[m] == 0) {
       if (!(compute->invoked_flag & INVOKED_VECTOR)) {
 	compute->compute_vector();
 	compute->invoked_flag |= INVOKED_VECTOR;
       }
       double *cvector = compute->vector;
       j = 0;
       for (i = nstart; i < nstop; i += nskip) {
 	vec[j] = cvector[i-1];
 	j += stride;
       }
       
     } else {
       if (!(compute->invoked_flag & INVOKED_ARRAY)) {
 	compute->compute_array();
 	compute->invoked_flag |= INVOKED_ARRAY;
       }
       double **carray = compute->array;
       int icol = argindex[m]-1;
       j = 0;
       for (i = nstart; i < nstop; i += nskip) {
 	vec[j] = carray[i-1][icol];
 	j += stride;
       }
     }
 
   // access fix fields, check if fix frequency is a match
     
   } else if (which[m] == FIX) {
     if (update->ntimestep % modify->fix[value2index[m]]->global_freq)
       error->all("Fix used in compute slice not computed at compatible time");
     Fix *fix = modify->fix[value2index[m]];
 
     if (argindex[m] == 0) {
       j = 0;
       for (i = nstart; i < nstop; i += nskip) {
 	vec[j] = fix->compute_vector(i-1);
 	j += stride;
       }
     } else {
       int icol = argindex[m]-1;
       j = 0;
       for (i = nstart; i < nstop; i += nskip) {
 	vec[j] = fix->compute_array(i-1,icol);
 	j += stride;
       }
     }
   }
 }
diff --git a/src/force.h b/src/force.h
index 9f20b4c66..df214c974 100644
--- a/src/force.h
+++ b/src/force.h
@@ -1,95 +1,99 @@
 /* ----------------------------------------------------------------------
    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.
 ------------------------------------------------------------------------- */
 
 #ifndef LMP_FORCE_H
 #define LMP_FORCE_H
 
 #include "pointers.h"
 
 namespace LAMMPS_NS {
 
 class Force : protected Pointers {
  public:
   double boltz;                      // Boltzmann constant (eng/degree-K)
   double mvv2e;                      // conversion of mv^2 to energy
   double ftm2v;                      // conversion of ft/m to velocity
   double mv2d;                       // conversion of mass/volume to density
   double nktv2p;                     // conversion of NkT/V to pressure
   double qqr2e;                      // conversion of q^2/r to energy
   double qe2f;                       // conversion of qE to force
   double vxmu2f;                     // conversion of vx dynamic-visc to force
   double xxt2kmu;                    // conversion of xx/t to kinematic-visc
   double dielectric;                 // dielectric constant
   double qqrd2e;                     // q^2/r to energy w/ dielectric constant
+  double e_mass;                     // electron mass
+  double hhmrr2e;                    // conversion of (hbar)^2/(mr^2) to energy
+  double mvh2r;                      // conversion of mv/hbar to distance 
+                                     // hbar = h/(2*pi)
 
   int newton,newton_pair,newton_bond;   // Newton's 3rd law settings
 
   class Pair *pair;
   char *pair_style;
 
   class Bond *bond;
   char *bond_style;
 
   class Angle *angle;
   char *angle_style;
 
   class Dihedral *dihedral;
   char *dihedral_style;
 
   class Improper *improper;
   char *improper_style;
 
   class KSpace *kspace;
   char *kspace_style;
                              // index [0] is not used in these arrays
   double special_lj[4];      // 1-2, 1-3, 1-4 prefactors for LJ
   double special_coul[4];    // 1-2, 1-3, 1-4 prefactors for Coulombics
   int special_angle;         // 0 if defined angles are ignored
                              // 1 if only weight 1,3 atoms if in an angle
   int special_dihedral;      // 0 if defined dihedrals are ignored
                              // 1 if only weight 1,4 atoms if in a dihedral
   int special_extra;         // extra space for added bonds
 
   Force(class LAMMPS *);
   ~Force();
   void init();
 
   void create_pair(const char *, char *suffix = NULL);
   class Pair *new_pair(const char *, char *, int &);
   class Pair *pair_match(const char *, int);
 
   void create_bond(const char *);
   class Bond *new_bond(const char *);
   class Bond *bond_match(const char *); 
 
   void create_angle(const char *);
   class Angle *new_angle(const char *);
 
   void create_dihedral(const char *);
   class Dihedral *new_dihedral(const char *);
 
   void create_improper(const char *);
   class Improper *new_improper(const char *);
 
   void create_kspace(int, char **);
 
   void set_special(int, char **);
   void bounds(char *, int, int &, int &);
   double numeric(char *);
   int inumeric(char *);
   bigint memory_usage();
 };
 
 }
 
 #endif
diff --git a/src/update.cpp b/src/update.cpp
index cc6f49b16..5718368d1 100644
--- a/src/update.cpp
+++ b/src/update.cpp
@@ -1,355 +1,379 @@
 /* ----------------------------------------------------------------------
    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 "lmptype.h"
 #include "string.h"
 #include "stdlib.h"
 #include "update.h"
 #include "style_integrate.h"
 #include "style_minimize.h"
 #include "neighbor.h"
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
 #include "domain.h"
 #include "region.h"
 #include "compute.h"
 #include "output.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 Update::Update(LAMMPS *lmp) : Pointers(lmp)
 {
   int n;
   char *str;
 
   ntimestep = 0;
   first_update = 0;
 
   whichflag = 0;
   firststep = laststep = 0;
   beginstep = endstep = 0;
   setupflag = 0;
   multireplica = 0;
 
   restrict_output = 0;
 
   eflag_global = vflag_global = -1;
 
   unit_style = NULL;
   set_units("lj");
 
   integrate_style = NULL;
   integrate = NULL;
   minimize_style = NULL;
   minimize = NULL;
 
   if (lmp->cuda) {
     str = (char *) "verlet/cuda";
     create_integrate(1,&str,NULL);
   } else {
     str = (char *) "verlet";
     create_integrate(1,&str,NULL);
   }
 
   str = (char *) "cg";
   create_minimize(1,&str);
 }
 
 /* ---------------------------------------------------------------------- */
 
 Update::~Update()
 {
   delete [] unit_style;
 
   delete [] integrate_style;
   delete integrate;
 
   delete [] minimize_style;
   delete minimize;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Update::init()
 {
   // if USER-CUDA mode is enabled:
   // integrate/minimize style must be CUDA variant
 
   if (whichflag == 1 && lmp->cuda)
     if (strstr(integrate_style,"cuda") == NULL)
       error->all("USER-CUDA mode requires CUDA variant of run style");
   if (whichflag == 2 && lmp->cuda)
     if (strstr(minimize_style,"cuda") == NULL)
       error->all("USER-CUDA mode requires CUDA variant of min style");
 
   // init the appropriate integrate and/or minimize class
   // if neither (e.g. from write_restart) then just return
 
   if (whichflag == 0) return;
   if (whichflag == 1) integrate->init();
   else if (whichflag == 2) minimize->init();
 
   // only set first_update if a run or minimize is being performed
 
   first_update = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Update::set_units(const char *style)
 {
   // physical constants from:
   // http://physics.nist.gov/cuu/Constants/Table/allascii.txt
   // using thermochemical calorie = 4.184 J
   
   if (strcmp(style,"lj") == 0) {
     force->boltz = 1.0;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
     force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 1.0;
     force->qe2f = 1.0;
     force->vxmu2f = 1.0;
     force->xxt2kmu = 1.0;
+    force->e_mass = 0.0;    // not yet set
+    force->hhmrr2e = 0.0;
+    force->mvh2r = 0.0;
+
     dt = 0.005;
     neighbor->skin = 0.3;
     
   } else if (strcmp(style,"real") == 0) {
     force->boltz = 0.0019872067;
     force->mvv2e = 48.88821291 * 48.88821291;
     force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
     force->mv2d = 1.0 / 0.602214179;
     force->nktv2p = 68568.415;
     force->qqr2e = 332.06371;
     force->qe2f = 23.060549; 
     force->vxmu2f = 1.4393264316e4;
     force->xxt2kmu = 0.1;
+    force->e_mass = 1.0/1836.1527556560675;
+    force->hhmrr2e = 0.0957018663603261;
+    force->mvh2r = 1.5339009481951;
+
     dt = 1.0;
     neighbor->skin = 2.0;
 
   } else if (strcmp(style,"metal") == 0) {
     force->boltz = 8.617343e-5;
     force->mvv2e = 1.0364269e-4;
     force->ftm2v = 1.0 / 1.0364269e-4;
     force->mv2d = 1.0 / 0.602214179;
     force->nktv2p = 1.6021765e6;
     force->qqr2e = 14.399645;
     force->qe2f = 1.0;
     force->vxmu2f = 0.6241509647;
     force->xxt2kmu = 1.0e-4;
+    force->e_mass = 0.0;    // not yet set
+    force->hhmrr2e = 0.0;
+    force->mvh2r = 0.0;
+
     dt = 0.001;
     neighbor->skin = 2.0;
 
   } else if (strcmp(style,"si") == 0) {
     force->boltz = 1.3806504e-23;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
     force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 8.9876e9;
     force->qe2f = 1.0;
     force->vxmu2f = 1.0;
     force->xxt2kmu = 1.0;
+    force->e_mass = 0.0;    // not yet set
+    force->hhmrr2e = 0.0;
+    force->mvh2r = 0.0;
+
     dt = 1.0e-8;
     neighbor->skin = 0.001;
 
   } else if (strcmp(style,"cgs") == 0) {
     force->boltz = 1.3806504e-16;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
     force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 1.0;
     force->qe2f = 1.0;
     force->vxmu2f = 1.0;
     force->xxt2kmu = 1.0;
+    force->e_mass = 0.0;    // not yet set
+    force->hhmrr2e = 0.0;
+    force->mvh2r = 0.0;
+
     dt = 1.0e-8;
     neighbor->skin = 0.1;
 
   } else if (strcmp(style,"electron") == 0) {
     force->boltz = 3.16681534e-6;  
     force->mvv2e = 1.06657236;     
     force->ftm2v = 0.937582899;    
     force->mv2d = 1.0;             
     force->nktv2p = 2.94210108e13; 
     force->qqr2e = 1.0;            
     force->qe2f = 1.94469051e-10;  
     force->vxmu2f = 3.39893149e1;  
     force->xxt2kmu = 3.13796367e-2;
+    force->e_mass = 0.0;    // not yet set
+    force->hhmrr2e = 0.0;
+    force->mvh2r = 0.0;
+
     dt = 0.001;
     neighbor->skin = 2.0;
     
   } else error->all("Illegal units command");
 
   delete [] unit_style;
   int n = strlen(style) + 1;
   unit_style = new char[n];
   strcpy(unit_style,style);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Update::create_integrate(int narg, char **arg, char *suffix)
 {
   if (narg < 1) error->all("Illegal run_style command");
 
   delete [] integrate_style;
   delete integrate;
 
   int sflag;
   new_integrate(arg[0],narg-1,&arg[1],suffix,sflag);
 
   if (sflag) {
     char estyle[256];
     sprintf(estyle,"%s/%s",arg[0],suffix);
     int n = strlen(estyle) + 1;
     integrate_style = new char[n];
     strcpy(integrate_style,estyle);
   } else {
     int n = strlen(arg[0]) + 1;
     integrate_style = new char[n];
     strcpy(integrate_style,arg[0]);
   }
 }
 
 /* ----------------------------------------------------------------------
    create the Integrate style, first with suffix appended
 ------------------------------------------------------------------------- */
 
 void Update::new_integrate(char *style, int narg, char **arg,
 			   char *suffix, int &sflag)
 {
   int success = 0;
 
   if (suffix && lmp->suffix_enable) {
     sflag = 1;
     char estyle[256];
     sprintf(estyle,"%s/%s",style,suffix);
     success = 1;
 
     if (0) return;
 
 #define INTEGRATE_CLASS
 #define IntegrateStyle(key,Class) \
     else if (strcmp(estyle,#key) == 0) integrate = new Class(lmp,narg,arg);
 #include "style_integrate.h"
 #undef IntegrateStyle
 #undef INTEGRATE_CLASS
 
     else success = 0;
   }
 
   if (!success) {
     sflag = 0;
 
     if (0) return;
 
 #define INTEGRATE_CLASS
 #define IntegrateStyle(key,Class) \
     else if (strcmp(style,#key) == 0) integrate = new Class(lmp,narg,arg);
 #include "style_integrate.h"
 #undef IntegrateStyle
 #undef INTEGRATE_CLASS
 
     else error->all("Illegal integrate style");
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Update::create_minimize(int narg, char **arg)
 {
   if (narg != 1) error->all("Illegal min_style command");
 
   delete [] minimize_style;
   delete minimize;
 
   if (0) return;      // dummy line to enable else-if macro expansion
 
 #define MINIMIZE_CLASS
 #define MinimizeStyle(key,Class) \
   else if (strcmp(arg[0],#key) == 0) minimize = new Class(lmp);
 #include "style_minimize.h"
 #undef MINIMIZE_CLASS
 
   else error->all("Illegal min_style command");
 
   int n = strlen(arg[0]) + 1;
   minimize_style = new char[n];
   strcpy(minimize_style,arg[0]);
 }
 
 /* ----------------------------------------------------------------------
    reset timestep from input script
    do not allow dump files or a restart to be defined
    do not allow any timestep-dependent fixes to be defined
    do not allow any dynamic regions to be defined
    reset eflag/vflag global so nothing will think eng/virial are current
    reset invoked flags of computes,
      so nothing will think they are current between runs
    clear timestep list of computes that store future invocation times
 ------------------------------------------------------------------------- */
 
 void Update::reset_timestep(int narg, char **arg)
 {
   if (narg != 1) error->all("Illegal reset_timestep command");
 
   for (int i = 0; i < output->ndump; i++)
     if (output->last_dump[i] >= 0)
       error->all("Cannot reset timestep with dump file already written to");
   if (output->restart && output->last_restart >= 0)
     error->all("Cannot reset timestep with restart file already written");
 
   for (int i = 0; i < modify->nfix; i++)
     if (modify->fix[i]->time_depend)
       error->all("Cannot reset timestep with a time-dependent fix defined");
 
   for (int i = 0; i < domain->nregion; i++)
     if (domain->regions[i]->dynamic_check())
       error->all("Cannot reset timestep with a dynamic region defined");
 
   eflag_global = vflag_global = -1;
 
   for (int i = 0; i < modify->ncompute; i++) {
     modify->compute[i]->invoked_scalar = -1;
     modify->compute[i]->invoked_vector = -1;
     modify->compute[i]->invoked_array = -1;
     modify->compute[i]->invoked_peratom = -1;
     modify->compute[i]->invoked_local = -1;
   }
 
   for (int i = 0; i < modify->ncompute; i++)
     if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
 
   ntimestep = ATOBIGINT(arg[0]);
   if (ntimestep < 0) error->all("Timestep must be >= 0");
   if (ntimestep > MAXBIGINT) error->all("Too big a timestep");
 }
 
 /* ----------------------------------------------------------------------
    memory usage of update and integrate/minimize
 ------------------------------------------------------------------------- */
 
 bigint Update::memory_usage()
 {
   bigint bytes = 0;
   if (whichflag == 1) bytes += integrate->memory_usage();
   else if (whichflag == 2) bytes += minimize->memory_usage();
   return bytes;
 }