diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index 86174a69ff30cb4843f1b981fc288ee01c7e25c1..e852e0abd452b707fc0eb25f0783cfb14bf85507 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -1856,9 +1856,21 @@ internal LAMMPS operations.  Note that LAMMPS classes are defined
 within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
 application.
 
-Library.cpp contains these functions for creating and destroying an
-instance of LAMMPS and sending it commands to execute.  See the
-documentation in the src/library.cpp file for details:
+The examples/COUPLE and python/examples directories have example C++
+and C and Python codes which show how a driver code can link to LAMMPS
+as a library, run LAMMPS on a subset of processors, grab data from
+LAMMPS, change it, and put it back into LAMMPS.
+
+The file src/library.cpp contains the following functions for creating
+and destroying an instance of LAMMPS and sending it commands to
+execute.  See the documentation in the src/library.cpp file for
+details.
+
+NOTE: You can write code for additional functions as needed to define
+how your code talks to LAMMPS and add them to src/library.cpp and
+src/library.h, as well as to the "Python
+interface"_Section_python.html.  The added functions can access or
+change any internal LAMMPS data you wish.
 
 void lammps_open(int, char **, MPI_Comm, void **)
 void lammps_open_no_mpi(int, char **, void **)
@@ -1932,11 +1944,12 @@ the library allocated previously via other function calls.  See
 comments in src/library.cpp file for which other functions need this
 clean-up.
 
-Library.cpp also contains these functions for extracting information
-from LAMMPS and setting value within LAMMPS.  Again, see the
-documentation in the src/library.cpp file for details, including
+The file src/library.cpp also contains these functions for extracting
+information from LAMMPS and setting value within LAMMPS.  Again, see
+the documentation in the src/library.cpp file for details, including
 which quantities can be queried by name:
 
+int lammps_extract_setting(void *, char *)
 void *lammps_extract_global(void *, char *)
 void lammps_extract_box(void *, double *, double *,
                         double *, double *, double *, int *, int *)
@@ -1945,55 +1958,77 @@ void *lammps_extract_compute(void *, char *, int, int)
 void *lammps_extract_fix(void *, char *, int, int, int, int)
 void *lammps_extract_variable(void *, char *, char *) :pre
 
-void lammps_reset_box(void *, double *, double *, double, double, double)
-int lammps_set_variable(void *, char *, char *) :pre
+The extract_setting() function returns info on the size
+of data types (e.g. 32-bit or 64-bit atom IDs) used
+by the LAMMPS executable (a compile-time choice).
+
+The other extract functions return a pointer to various global or
+per-atom quantities stored in LAMMPS or to values calculated by a
+compute, fix, or variable.  The pointer returned by the
+extract_global() function can be used as a permanent reference to a
+value which may change.  For the extract_atom() method, see the
+extract() method in the src/atom.cpp file for a list of valid per-atom
+properties.  New names could easily be added if the property you want
+is not listed.  For the other extract functions, the underlying
+storage may be reallocated as LAMMPS runs, so you need to re-call the
+function to assure a current pointer or returned value(s).
 
 double lammps_get_thermo(void *, char *)
-int lammps_get_natoms(void *)
-void lammps_gather_atoms(void *, double *)
-void lammps_scatter_atoms(void *, double *) :pre
-void lammps_create_atoms(void *, int, tagint *, int *, double *, double *,
-                         imageint *, int) :pre
+int lammps_get_natoms(void *) :pre
 
-The extract functions return a pointer to various global or per-atom
-quantities stored in LAMMPS or to values calculated by a compute, fix,
-or variable.  The pointer returned by the extract_global() function
-can be used as a permanent reference to a value which may change.  For
-the extract_atom() method, see the extract() method in the
-src/atom.cpp file for a list of valid per-atom properties.  New names
-could easily be added if the property you want is not listed.  For the
-other extract functions, the underlying storage may be reallocated as
-LAMMPS runs, so you need to re-call the function to assure a current
-pointer or returned value(s).
+int lammps_set_variable(void *, char *, char *)
+void lammps_reset_box(void *, double *, double *, double, double, double) :pre
 
-The lammps_reset_box() function resets the size and shape of the
-simulation box, e.g. as part of restoring a previously extracted and
-saved state of a simulation.
+The lammps_get_thermo() function returns the current value of a thermo
+keyword as a double precision value.
+
+The lammps_get_natoms() function returns the total number of atoms in
+the system and can be used by the caller to allocate memory for the
+lammps_gather_atoms() and lammps_scatter_atoms() functions.
 
 The lammps_set_variable() function can set an existing string-style
 variable to a new string value, so that subsequent LAMMPS commands can
 access the variable.
 
-The lammps_get_thermo() function returns the current value of a thermo
-keyword as a double precision value.
+The lammps_reset_box() function resets the size and shape of the
+simulation box, e.g. as part of restoring a previously extracted and
+saved state of a simulation.
 
-The lammps_get_natoms() function returns the total number of atoms in
-the system and can be used by the caller to allocate space for the
-lammps_gather_atoms() and lammps_scatter_atoms() functions.  The
-gather function collects peratom info of the requested type (atom
-coords, types, forces, etc) from all processors, orders them by atom
-ID, and returns a full list to each calling processor.  The scatter
-function does the inverse.  It distributes the same peratom values,
-passed by the caller, to each atom owned by individual processors.
-Both methods are thus a means to extract or assign (overwrite) any
-peratom quantities within LAMMPS.  See the extract() method in the
-src/atom.cpp file for a list of valid per-atom properties.  New names
-could easily be added if the property you want is not listed.
-A special treatment is applied for accessing image flags via the
-"image" property. Image flags are stored in a packed format with all
-three image flags stored in a single integer. When signaling to access
-the image flags as 3 individual values per atom instead of 1, the data
-is transparently packed or unpacked by the library interface.
+void lammps_gather_atoms(void *, char *, int, int, void *)
+void lammps_gather_atoms_concat(void *, char *, int, int, void *)
+void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *)
+void lammps_scatter_atoms(void *, char *, int, int, void *)
+void lammps_scatter_atoms_subset(void *, char *, int, int, int, int *, void *) :pre
+
+void lammps_create_atoms(void *, int, tagint *, int *, double *, double *,
+                         imageint *, int) :pre
+
+The gather functions collect peratom info of the requested type (atom
+coords, atom types, forces, etc) from all processors, and returns the
+same vector of values to each callling processor.  The scatter
+functions do the inverse.  They distribute a vector of peratom values,
+passed by all calling processors, to invididual atoms, which may be
+owned by different processos.
+
+The lammps_gather_atoms() function does this for all N atoms in the
+system, ordered by atom ID, from 1 to N.  The
+lammps_gather_atoms_concat() function does it for all N atoms, but
+simply concatenates the subset of atoms owned by each processor.  The
+resulting vector is not ordered by atom ID.  Atom IDs can be requetsed
+by the same function if the caller needs to know the ordering.  The
+lammps_gather_subset() function allows the caller to request values
+for only a subset of atoms (identified by ID).
+For all 3 gather function, per-atom image flags can be retrieved in 2 ways.
+If the count is specified as 1, they are returned 
+in a packed format with all three image flags stored in a single integer.
+If the count is specified as 3, the values are unpacked into xyz flags
+by the library before returning them.
+
+The lammps_scatter_atoms() function takes a list of values for all N
+atoms in the system, ordered by atom ID, from 1 to N, and assigns
+those values to each atom in the system.  The
+lammps_scatter_atoms_subset() function takes a subset of IDs as an
+argument and only scatters those values to the owning atoms.
 
 The lammps_create_atoms() function takes a list of N atoms as input
 with atom types and coords (required), an optionally atom IDs and
@@ -2005,17 +2040,6 @@ of a simulation.  Additional properties for the new atoms can then be
 assigned via the lammps_scatter_atoms() or lammps_extract_atom()
 functions.
 
-The examples/COUPLE and python directories have example C++ and C and
-Python codes which show how a driver code can link to LAMMPS as a
-library, run LAMMPS on a subset of processors, grab data from LAMMPS,
-change it, and put it back into LAMMPS.
-
-NOTE: You can write code for additional functions as needed to define
-how your code talks to LAMMPS and add them to src/library.cpp and
-src/library.h, as well as to the "Python
-interface"_Section_python.html.  The added functions can access or
-change any LAMMPS data you wish.
-
 :line
 
 6.20 Calculating thermal conductivity :link(howto_20),h4
diff --git a/doc/src/Section_python.txt b/doc/src/Section_python.txt
index 09a84fa3957c49fb6e02d6df9c69b725371d0506..5427cd67f28d36bb0344012e6f7cfe593f7aca88 100644
--- a/doc/src/Section_python.txt
+++ b/doc/src/Section_python.txt
@@ -551,11 +551,14 @@ Python script, as follows:
 from lammps import lammps :pre
 
 These are the methods defined by the lammps module.  If you look at
-the files src/library.cpp and src/library.h you will see that they
+the files src/library.cpp and src/library.h you will see they
 correspond one-to-one with calls you can make to the LAMMPS library
 from a C++ or C or Fortran program, and which are described in
 "Section 6.19"_Section_howto.html#howto_19 of the manual.
 
+The python/examples directory has Python scripts which show how Python
+can run LAMMPS, grab data, change it, and put it back into LAMMPS.
+
 lmp = lammps()           # create a LAMMPS object using the default liblammps.so library
                          # 4 optional args are allowed: name, cmdargs, ptr, comm
 lmp = lammps(ptr=lmpptr) # use lmpptr as previously created LAMMPS object
@@ -565,18 +568,22 @@ lmp = lammps(name="g++",cmdargs=list)    # add LAMMPS command-line args, e.g. li
 
 lmp.close()              # destroy a LAMMPS object :pre
 
-version = lmp.version()  # return the numerical version id, e.g. LAMMPS 2 Sep 2015 -> 20150902
+version = lmp.version()  # return the numerical version id, e.g. LAMMPS 2 Sep 2015 -> 20150902 :pre
 
 lmp.file(file)           # run an entire input script, file = "in.lj"
-lmp.command(cmd)         # invoke a single LAMMPS command, cmd = "run 100" :pre
+lmp.command(cmd)         # invoke a single LAMMPS command, cmd = "run 100"
 lmp.commands_list(cmdlist)     # invoke commands in cmdlist = ["run 10", "run 20"]
-lmp.commands_string(multicmd)  # invoke commands in multicmd = "run 10\nrun 20"
+lmp.commands_string(multicmd)  # invoke commands in multicmd = "run 10\nrun 20" :pre
+
+size = lmp.extract_setting(name)     # return data type info :pre
 
 xlo = lmp.extract_global(name,type)  # extract a global quantity
                                      # name = "boxxlo", "nlocal", etc
                                      # type = 0 = int
                                      #        1 = double :pre
 
+boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()  # extract box info :pre
+
 coords = lmp.extract_atom(name,type)      # extract a per-atom quantity
                                           # name = "x", "type", etc
                                           # type = 0 = vector of ints
@@ -601,16 +608,23 @@ var = lmp.extract_variable(name,group,flag)  # extract value(s) from a variable
                                              # flag = 0 = equal-style variable
                                              #        1 = atom-style variable :pre
 
-flag = lmp.set_variable(name,value)       # set existing named string-style variable to value, flag = 0 if successful
 value = lmp.get_thermo(name)              # return current value of a thermo keyword
+natoms = lmp.get_natoms()                 # total # of atoms as int :pre
+
+flag = lmp.set_variable(name,value)       # set existing named string-style variable to value, flag = 0 if successful
+lmp.reset_box(boxlo,boxhi,xy,yz,xz)       # reset the simulation box size :pre
 
-natoms = lmp.get_natoms()                 # total # of atoms as int
 data = lmp.gather_atoms(name,type,count)  # return per-atom property of all atoms gathered into data, ordered by atom ID
                                           # name = "x", "charge", "type", etc
-                                          # count = # of per-atom values, 1 or 3, etc
+data = lmp.gather_atoms_concat(name,type,count)  # ditto, but concatenated atom values from each proc (unordered)
+data = lmp.gather_atoms_subset(name,type,count,ndata,ids)  # ditto, but for subset of Ndata atoms with IDs :pre
+
 lmp.scatter_atoms(name,type,count,data)   # scatter per-atom property to all atoms from data, ordered by atom ID
                                           # name = "x", "charge", "type", etc
                                           # count = # of per-atom values, 1 or 3, etc :pre
+lmp.scatter_atoms_subset(name,type,count,ndata,ids,data)  # ditto, but for subset of Ndata atoms with IDs :pre
+
+lmp.create_atoms(n,ids,types,x,v,image,shrinkexceed)   # create N atoms with IDs, types, x, v, and image flags :pre
 
 :line
 
@@ -655,9 +669,10 @@ The file(), command(), commands_list(), commands_string() methods
 allow an input script, a single command, or multiple commands to be
 invoked.
 
-The extract_global(), extract_atom(), extract_compute(),
-extract_fix(), and extract_variable() methods return values or
-pointers to data structures internal to LAMMPS.
+The extract_setting(), extract_global(), extract_box(),
+extract_atom(), extract_compute(), extract_fix(), and
+extract_variable() methods return values or pointers to data
+structures internal to LAMMPS.
 
 For extract_global() see the src/library.cpp file for the list of
 valid names.  New names could easily be added.  A double or integer is
@@ -697,60 +712,46 @@ doubles is returned, one value per atom, which you can use via normal
 Python subscripting. The values will be zero for atoms not in the
 specified group.
 
+The get_thermo() method returns returns the current value of a thermo
+keyword as a float.
+
 The get_natoms() method returns the total number of atoms in the
 simulation, as an int.
 
-The gather_atoms() method allows any per-atom property (coordinates,
-velocities, etc) to be extracted from LAMMPS.  It returns a ctypes
-vector of ints or doubles as specified by type, of length
-count*natoms, for the named property for all atoms in the simulation.
-The data is ordered by count and then by atom ID.  See the extract()
-method in the src/atom.cpp file for a list of valid names.  Again, new
-names could easily be added if the property you want is missing.  The
-vector can be used via normal Python subscripting.  If atom IDs are
-not consecutively ordered within LAMMPS, a None is returned as
-indication of an error. A special treatment is applied for image flags
-stored in the "image" property. All three image flags are stored in
-a packed format in a single integer, so count would be 1 to retrieve
-that integer, however also a count value of 3 can be used and then
-the image flags will be unpacked into 3 individual integers, ordered
-in a similar fashion as coordinates.
-
-Note that the data structure gather_atoms("x") returns is different
-from the data structure returned by extract_atom("x") in four ways.
-(1) Gather_atoms() returns a vector which you index as x\[i\];
-extract_atom() returns an array which you index as x\[i\]\[j\].  (2)
-Gather_atoms() orders the atoms by atom ID while extract_atom() does
-not.  (3) Gathert_atoms() returns a list of all atoms in the
-simulation; extract_atoms() returns just the atoms local to each
-processor.  (4) Finally, the gather_atoms() data structure is a copy
-of the atom coords stored internally in LAMMPS, whereas extract_atom()
-returns an array that effectively points directly to the internal
-data.  This means you can change values inside LAMMPS from Python by
-assigning a new values to the extract_atom() array.  To do this with
-the gather_atoms() vector, you need to change values in the vector,
-then invoke the scatter_atoms() method.
-
-The scatter_atoms() method allows any per-atom property (coordinates,
-velocities, etc) to be inserted into LAMMPS, overwriting the current
-property.  It takes a vector of ints or doubles as specified by type,
-of length count*natoms, for the named property for all atoms in the
-simulation.  The data should be ordered by count and then by atom ID.
-See the extract() method in the src/atom.cpp file for a list of valid
-names.  Again, new names could easily be added if the property you
-want is missing.  It uses the vector of data to overwrite the
-corresponding properties for each atom inside LAMMPS.  This requires
-LAMMPS to have its "map" option enabled; see the
-"atom_modify"_atom_modify.html command for details.  If it is not, or
-if atom IDs are not consecutively ordered, no coordinates are reset.
-Similar as for gather_atoms() a special treatment is applied for image
-flags, which can be provided in packed (count = 1) or unpacked (count = 3)
-format and in the latter case, they will be packed before applied to
-atoms.
-
-The array of coordinates passed to scatter_atoms() must be a ctypes
-vector of ints or doubles, allocated and initialized something like
-this:
+The set_variable() methosd sets an existing string-style variable to a
+new string value, so that subsequent LAMMPS commands can access the
+variable.
+
+The reset_box() emthods resets the size and shape of the simulation
+box, e.g. as part of restoring a previously extracted and saved state
+of a simulation.
+
+The gather methods collect peratom info of the requested type (atom
+coords, atom types, forces, etc) from all processors, and returns the
+same vector of values to each callling processor.  The scatter
+functions do the inverse.  They distribute a vector of peratom values,
+passed by all calling processors, to invididual atoms, which may be
+owned by different processos.
+
+Note that the data returned by the gather methods,
+e.g. gather_atoms("x"), is different from the data structure returned
+by extract_atom("x") in four ways.  (1) Gather_atoms() returns a
+vector which you index as x\[i\]; extract_atom() returns an array
+which you index as x\[i\]\[j\].  (2) Gather_atoms() orders the atoms
+by atom ID while extract_atom() does not.  (3) Gather_atoms() returns
+a list of all atoms in the simulation; extract_atoms() returns just
+the atoms local to each processor.  (4) Finally, the gather_atoms()
+data structure is a copy of the atom coords stored internally in
+LAMMPS, whereas extract_atom() returns an array that effectively
+points directly to the internal data.  This means you can change
+values inside LAMMPS from Python by assigning a new values to the
+extract_atom() array.  To do this with the gather_atoms() vector, you
+need to change values in the vector, then invoke the scatter_atoms()
+method.
+
+For the scatter methods, the array of coordinates passed to must be a
+ctypes vector of ints or doubles, allocated and initialized something
+like this:
 
 from ctypes import *
 natoms = lmp.get_natoms()
@@ -765,7 +766,7 @@ x\[n3-1\] = z coord of atom with ID natoms
 lmp.scatter_atoms("x",1,3,x) :pre
 
 Alternatively, you can just change values in the vector returned by
-gather_atoms("x",1,3), since it is a ctypes vector of doubles.
+the gather methods, since they are also ctypes vectors.
 
 :line
 
diff --git a/python/examples/simple.py b/python/examples/simple.py
index ff777c8a0904089ec9e900474775d2759f8ea045..632dde9ed3bcc1195bc19ce7587cedcfde9019ab 100755
--- a/python/examples/simple.py
+++ b/python/examples/simple.py
@@ -80,12 +80,59 @@ lmp.commands_list(cmds)
 # initial thermo should be same as step 20
 
 natoms = lmp.get_natoms()
-type = natoms*[1]
+atype = natoms*[1]
 
 lmp.command("delete_atoms group all");
-lmp.create_atoms(natoms,None,type,x,v);
+lmp.create_atoms(natoms,None,atype,x,v);
 lmp.command("run 10");
 
+############
+# test of new gather/scatter and box extract/reset methods
+# can try this in parallel and with/without atom_modify sort enabled
+
+lmp.command("write_dump all custom tmp.simple id type x y z fx fy fz");
+
+x = lmp.gather_atoms("x",1,3)
+f = lmp.gather_atoms("f",1,3)
+
+if me == 0: print("Gather XF:",x[3],x[9],f[3],f[9])
+
+ids = lmp.gather_atoms_concat("id",0,1)
+x = lmp.gather_atoms_concat("x",1,3)
+f = lmp.gather_atoms_concat("f",1,3)
+
+if me == 0: print("Gather concat XF:",ids[0],ids[1],x[0],x[3],f[0],f[3])
+
+ids = (2*ctypes.c_int)()
+ids[0] = 2
+ids[1] = 4
+x = lmp.gather_atoms_subset("x",1,3,2,ids)
+f = lmp.gather_atoms_subset("f",1,3,2,ids)
+
+if me == 0: print("Gather subset XF:",x[0],x[3],f[0],f[3])
+
+x[0] = -1.0
+x[1] = 0.0
+x[2] = 0.0
+x[3] = -2.0
+x[4] = 0.0
+x[5] = 0.0
+ids[0] = 100
+ids[1] = 200
+lmp.scatter_atoms_subset("x",1,3,2,ids,x)
+
+x = lmp.gather_atoms("x",1,3)
+if me == 0: print("Gather post scatter subset:",
+                  x[3],x[9],x[297],x[298],x[299],x[597],x[598],x[599])
+
+boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()
+if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change)
+
+lmp.reset_box([0,0,0],[10,10,8],0,0,0)
+
+boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()
+if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change)
+
 # uncomment if running in parallel via Pypar
 #print("Proc %d out of %d procs has" % (me,nprocs), lmp)
 
diff --git a/python/lammps.py b/python/lammps.py
index fdb898488a4d23faab66301b65c40eae9abc4017..e798ef60718d7cd9763641f799b0dd8330da1f53 100644
--- a/python/lammps.py
+++ b/python/lammps.py
@@ -94,6 +94,39 @@ class lammps(object):
         if not name: self.lib = CDLL("liblammps.so",RTLD_GLOBAL)
         else: self.lib = CDLL("liblammps_%s.so" % name,RTLD_GLOBAL)
 
+    # define ctypes API for each library method
+    # NOTE: should add one of these for each lib function
+
+    self.lib.lammps_extract_box.argtypes = \
+      [c_void_p,POINTER(c_double),POINTER(c_double),
+       POINTER(c_double),POINTER(c_double),POINTER(c_double),
+       POINTER(c_int),POINTER(c_int)]
+    self.lib.lammps_extract_box.restype = None
+
+    self.lib.lammps_reset_box.argtypes = \
+      [c_void_p,POINTER(c_double),POINTER(c_double),c_double,c_double,c_double]
+    self.lib.lammps_reset_box.restype = None
+
+    self.lib.lammps_gather_atoms.argtypes = \
+      [c_void_p,c_char_p,c_int,c_int,c_void_p]
+    self.lib.lammps_gather_atoms.restype = None
+    
+    self.lib.lammps_gather_atoms_concat.argtypes = \
+      [c_void_p,c_char_p,c_int,c_int,c_void_p]
+    self.lib.lammps_gather_atoms_concat.restype = None
+    
+    self.lib.lammps_gather_atoms_subset.argtypes = \
+      [c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]
+    self.lib.lammps_gather_atoms_subset.restype = None
+    
+    self.lib.lammps_scatter_atoms.argtypes = \
+      [c_void_p,c_char_p,c_int,c_int,c_void_p]
+    self.lib.lammps_scatter_atoms.restype = None
+
+    self.lib.lammps_scatter_atoms_subset.argtypes = \
+      [c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]
+    self.lib.lammps_scatter_atoms_subset.restype = None
+
     # if no ptr provided, create an instance of LAMMPS
     #   don't know how to pass an MPI communicator from PyPar
     #   but we can pass an MPI communicator from mpi4py v2.0.0 and later
@@ -177,6 +210,8 @@ class lammps(object):
     self.c_tagint = get_ctypes_int(self.extract_setting("tagint"))
     self.c_imageint = get_ctypes_int(self.extract_setting("imageint"))
 
+  # shut-down LAMMPS instance
+  
   def __del__(self):
     if self.lmp and self.opened:
       self.lib.lammps_close(self.lmp)
@@ -219,10 +254,16 @@ class lammps(object):
   # send a string of commands
 
   def commands_string(self,multicmd):
-    if type(multicmd) is str:
-        multicmd = multicmd.encode()
+    if type(multicmd) is str: multicmd = multicmd.encode()
     self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd))
     
+  # extract lammps type byte sizes
+
+  def extract_setting(self, name):
+    if name: name = name.encode()
+    self.lib.lammps_extract_atom.restype = c_int
+    return int(self.lib.lammps_extract_setting(self.lmp,name))
+
   # extract global info
     
   def extract_global(self,name,type):
@@ -235,8 +276,35 @@ class lammps(object):
     ptr = self.lib.lammps_extract_global(self.lmp,name)
     return ptr[0]
 
+  # extract global info
+    
+  def extract_box(self):
+    boxlo = (3*c_double)()
+    boxhi = (3*c_double)()
+    xy = c_double()
+    yz = c_double()
+    xz = c_double()
+    periodicity = (3*c_int)()
+    box_change = c_int()
+    
+    self.lib.lammps_extract_box(self.lmp,boxlo,boxhi,
+                                byref(xy),byref(yz),byref(xz),
+                                periodicity,byref(box_change))
+    
+    boxlo = boxlo[:3]
+    boxhi = boxhi[:3]
+    xy = xy.value
+    yz = yz.value
+    xz = xz.value
+    periodicity = periodicity[:3]
+    box_change = box_change.value
+    
+    return boxlo,boxhi,xy,yz,xz,periodicity,box_change
+    
   # extract per-atom info
-  
+  # NOTE: need to insure are converting to/from correct Python type
+  #   e.g. for Python list or NumPy or ctypes
+
   def extract_atom(self,name,type):
     if name: name = name.encode()
     if type == 0:
@@ -250,14 +318,7 @@ class lammps(object):
     else: return None
     ptr = self.lib.lammps_extract_atom(self.lmp,name)
     return ptr
-
-  # extract lammps type byte sizes
-
-  def extract_setting(self, name):
-    if name: name = name.encode()
-    self.lib.lammps_extract_atom.restype = c_int
-    return int(self.lib.lammps_extract_setting(self.lmp,name))
-
+    
   @property
   def numpy(self):
     if not self._numpy:
@@ -378,15 +439,6 @@ class lammps(object):
       return result
     return None
 
-  # set variable value
-  # value is converted to string
-  # returns 0 for success, -1 if failed
-
-  def set_variable(self,name,value):
-    if name: name = name.encode()
-    if value: value = str(value).encode()
-    return self.lib.lammps_set_variable(self.lmp,name,value)
-
   # return current value of thermo keyword
 
   def get_thermo(self,name):
@@ -399,14 +451,30 @@ class lammps(object):
   def get_natoms(self):
     return self.lib.lammps_get_natoms(self.lmp)
 
-  # return vector of atom properties gathered across procs, ordered by atom ID
+  # set variable value
+  # value is converted to string
+  # returns 0 for success, -1 if failed
+
+  def set_variable(self,name,value):
+    if name: name = name.encode()
+    if value: value = str(value).encode()
+    return self.lib.lammps_set_variable(self.lmp,name,value)
+
+  # reset simulation box size
+
+  def reset_box(self,boxlo,boxhi,xy,yz,xz):
+    cboxlo = (3*c_double)(*boxlo)
+    cboxhi = (3*c_double)(*boxhi)
+    self.lib.lammps_reset_box(self.lmp,cboxlo,cboxhi,xy,yz,xz)
+    
+  # return vector of atom properties gathered across procs
+  # 3 variants to match src/library.cpp
   # name = atom property recognized by LAMMPS in atom->extract()
   # type = 0 for integer values, 1 for double values
   # count = number of per-atom valus, 1 for type or charge, 3 for x or f
   # returned data is a 1d vector - doc how it is ordered?
-  # NOTE: how could we insure are converting to correct Python type
-  #   e.g. for Python list or NumPy, etc
-  #   ditto for extract_atom() above
+  # NOTE: need to insure are converting to/from correct Python type
+  #   e.g. for Python list or NumPy or ctypes
   
   def gather_atoms(self,name,type,count):
     if name: name = name.encode()
@@ -419,19 +487,47 @@ class lammps(object):
       self.lib.lammps_gather_atoms(self.lmp,name,type,count,data)
     else: return None
     return data
+    
+  def gather_atoms_concat(self,name,type,count):
+    if name: name = name.encode()
+    natoms = self.lib.lammps_get_natoms(self.lmp)
+    if type == 0:
+      data = ((count*natoms)*c_int)()
+      self.lib.lammps_gather_atoms_concat(self.lmp,name,type,count,data)
+    elif type == 1:
+      data = ((count*natoms)*c_double)()
+      self.lib.lammps_gather_atoms_concat(self.lmp,name,type,count,data)
+    else: return None
+    return data
+
+  def gather_atoms_subset(self,name,type,count,ndata,ids):
+    if name: name = name.encode()
+    if type == 0:
+      data = ((count*ndata)*c_int)()
+      self.lib.lammps_gather_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
+    elif type == 1:
+      data = ((count*ndata)*c_double)()
+      self.lib.lammps_gather_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
+    else: return None
+    return data
 
-  # scatter vector of atom properties across procs, ordered by atom ID
+  # scatter vector of atom properties across procs
+  # 2 variants to match src/library.cpp
   # name = atom property recognized by LAMMPS in atom->extract()
   # type = 0 for integer values, 1 for double values
   # count = number of per-atom valus, 1 for type or charge, 3 for x or f
   # assume data is of correct type and length, as created by gather_atoms()
-  # NOTE: how could we insure are passing correct type to LAMMPS
-  # e.g. for Python list or NumPy, etc
-    
+  # NOTE: need to insure are converting to/from correct Python type
+  #   e.g. for Python list or NumPy or ctypes
+  
   def scatter_atoms(self,name,type,count,data):
     if name: name = name.encode()
     self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data)
 
+  def scatter_atoms_subset(self,name,type,count,ndata,ids,data):
+    if name: name = name.encode()
+    self.lib.lammps_scatter_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
+
   # create N atoms on all procs
   # N = global number of atoms
   # id = ID of each atom (optional, can be None)
@@ -457,8 +553,8 @@ class lammps(object):
 
     type_lmp = (c_int * n)()
     type_lmp[:] = type
-    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v,image_lmp,shrinkexceed)
-
+    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v,image_lmp,
+                                 shrinkexceed)
 
   @property
   def uses_exceptions(self):
diff --git a/src/library.cpp b/src/library.cpp
index 56b4bc596d86e7033d05d4a25e13511d82d14115..07c4cae9c16327b5efd618391d3096f1551ac4e3 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -687,33 +687,40 @@ void *lammps_extract_variable(void *ptr, char *name, char *group)
   return NULL;
 }
 
-
 /* ----------------------------------------------------------------------
-   reset simulation box parameters
-   see domain.h for definition of these arguments
-   assumes domain->set_initial_box() has been invoked previously
+   return the current value of a thermo keyword as a double
+   unlike lammps_extract_global() this does not give access to the
+     storage of the data in question
+   instead it triggers the Thermo class to compute the current value
+     and returns it
 ------------------------------------------------------------------------- */
 
-void lammps_reset_box(void *ptr, double *boxlo, double *boxhi,
-                      double xy, double yz, double xz)
+double lammps_get_thermo(void *ptr, char *name)
 {
   LAMMPS *lmp = (LAMMPS *) ptr;
-  Domain *domain = lmp->domain;
+  double dval = 0.0;
 
-  domain->boxlo[0] = boxlo[0];
-  domain->boxlo[1] = boxlo[1];
-  domain->boxlo[2] = boxlo[2];
-  domain->boxhi[0] = boxhi[0];
-  domain->boxhi[1] = boxhi[1];
-  domain->boxhi[2] = boxhi[2];
+  BEGIN_CAPTURE
+  {
+    lmp->output->thermo->evaluate_keyword(name,&dval);
+  }
+  END_CAPTURE
 
-  domain->xy = xy;
-  domain->yz = yz;
-  domain->xz = xz;
+  return dval;
+}
 
-  domain->set_global_box();
-  lmp->comm->set_proc_grid();
-  domain->set_local_box();
+/* ----------------------------------------------------------------------
+   return the total number of atoms in the system
+   useful before call to lammps_get_atoms() so can pre-allocate vector
+------------------------------------------------------------------------- */
+
+int lammps_get_natoms(void *ptr)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  if (lmp->atom->natoms > MAXSMALLINT) return 0;
+  int natoms = static_cast<int> (lmp->atom->natoms);
+  return natoms;
 }
 
 /* ----------------------------------------------------------------------
@@ -737,50 +744,52 @@ int lammps_set_variable(void *ptr, char *name, char *str)
 }
 
 /* ----------------------------------------------------------------------
-   return the current value of a thermo keyword as a double
-   unlike lammps_extract_global() this does not give access to the
-     storage of the data in question
-   instead it triggers the Thermo class to compute the current value
-     and returns it
+   reset simulation box parameters
+   see domain.h for definition of these arguments
+   assumes domain->set_initial_box() has been invoked previously
 ------------------------------------------------------------------------- */
 
-double lammps_get_thermo(void *ptr, char *name)
+void lammps_reset_box(void *ptr, double *boxlo, double *boxhi,
+                      double xy, double yz, double xz)
 {
   LAMMPS *lmp = (LAMMPS *) ptr;
-  double dval = 0.0;
-
-  BEGIN_CAPTURE
-  {
-    lmp->output->thermo->evaluate_keyword(name,&dval);
-  }
-  END_CAPTURE
-
-  return dval;
-}
+  Domain *domain = lmp->domain;
 
-/* ----------------------------------------------------------------------
-   return the total number of atoms in the system
-   useful before call to lammps_get_atoms() so can pre-allocate vector
-------------------------------------------------------------------------- */
+  domain->boxlo[0] = boxlo[0];
+  domain->boxlo[1] = boxlo[1];
+  domain->boxlo[2] = boxlo[2];
+  domain->boxhi[0] = boxhi[0];
+  domain->boxhi[1] = boxhi[1];
+  domain->boxhi[2] = boxhi[2];
 
-int lammps_get_natoms(void *ptr)
-{
-  LAMMPS *lmp = (LAMMPS *) ptr;
+  domain->xy = xy;
+  domain->yz = yz;
+  domain->xz = xz;
 
-  if (lmp->atom->natoms > MAXSMALLINT) return 0;
-  int natoms = static_cast<int> (lmp->atom->natoms);
-  return natoms;
+  domain->set_global_box();
+  lmp->comm->set_proc_grid();
+  domain->set_local_box();
 }
 
 /* ----------------------------------------------------------------------
-   gather the named atom-based entity across all processors
-   atom IDs must be consecutive from 1 to N
+   gather the named atom-based entity for all atoms
+     return it in user-allocated data
+   data will be ordered by atom ID
+     requirement for consecutive atom IDs (1 to N)
+   see gather_atoms_concat() to return data for all atoms, unordered
+   see gather_atoms_subset() to return data for only a subset of atoms
    name = desired quantity, e.g. x or charge
    type = 0 for integer values, 1 for double values
    count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+     use count = 3 with "image" if want single image flag unpacked into xyz
    return atom-based values in 1d data, ordered by count, then by atom ID
      e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
      data must be pre-allocated by caller to correct length
+     correct length = count*Natoms, as queried by get_natoms()
+   method:
+     alloc and zero count*Natom length vector
+     loop over Nlocal to fill vector with my values
+     Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
 void lammps_gather_atoms(void *ptr, char *name,
@@ -790,7 +799,10 @@ void lammps_gather_atoms(void *ptr, char *name,
 
   BEGIN_CAPTURE
   {
+    int i,j,offset;
+
     // error if tags are not defined or not consecutive
+    // NOTE: test that name = image or ids is not a 64-bit int in code?
 
     int flag = 0;
     if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0)
@@ -804,11 +816,10 @@ void lammps_gather_atoms(void *ptr, char *name,
 
     int natoms = static_cast<int> (lmp->atom->natoms);
 
-    int i,j,offset;
     void *vptr = lmp->atom->extract(name);
-    if(vptr == NULL) {
-        lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name");
-        return;
+    if (vptr == NULL) {
+      lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name");
+      return;
     }
 
     // copy = Natom length vector of per-atom values
@@ -830,10 +841,11 @@ void lammps_gather_atoms(void *ptr, char *name,
       tagint *tag = lmp->atom->tag;
       int nlocal = lmp->atom->nlocal;
 
-      if (count == 1)
+      if (count == 1) {
         for (i = 0; i < nlocal; i++)
           copy[tag[i]-1] = vector[i];
-      else if (imgunpack) {
+
+      } else if (imgunpack) {
         for (i = 0; i < nlocal; i++) {
           offset = count*(tag[i]-1);
           const int image = vector[i];
@@ -841,12 +853,14 @@ void lammps_gather_atoms(void *ptr, char *name,
           copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX;
           copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX;
         }
-      } else
+
+      } else {
         for (i = 0; i < nlocal; i++) {
           offset = count*(tag[i]-1);
           for (j = 0; j < count; j++)
             copy[offset++] = array[i][j];
         }
+      }
 
       MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
       lmp->memory->destroy(copy);
@@ -867,6 +881,7 @@ void lammps_gather_atoms(void *ptr, char *name,
       if (count == 1) {
         for (i = 0; i < nlocal; i++)
           copy[tag[i]-1] = vector[i];
+
       } else {
         for (i = 0; i < nlocal; i++) {
           offset = count*(tag[i]-1);
@@ -883,13 +898,305 @@ void lammps_gather_atoms(void *ptr, char *name,
 }
 
 /* ----------------------------------------------------------------------
-   scatter the named atom-based entity across all processors
-   atom IDs must be consecutive from 1 to N
+   gather the named atom-based entity for all atoms
+     return it in user-allocated data
+   data will be a concatenation of chunks of each proc's atoms,
+     in whatever order the atoms are on each proc
+     no requirement for consecutive atom IDs (1 to N)
+     can do a gather_atoms_concat for "id" if need to know atom IDs
+   see gather_atoms() to return data ordered by consecutive atom IDs
+   see gather_atoms_subset() to return data for only a subset of atoms
    name = desired quantity, e.g. x or charge
    type = 0 for integer values, 1 for double values
    count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+     use count = 3 with "image" if want single image flag unpacked into xyz
+   return atom-based values in 1d data, ordered by count, then by atom
+     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     data must be pre-allocated by caller to correct length
+     correct length = count*Natoms, as queried by get_natoms()
+   method:
+     Allgather Nlocal atoms from each proc into data
+------------------------------------------------------------------------- */
+
+void lammps_gather_atoms_concat(void *ptr, char *name,
+                                int type, int count, void *data)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  BEGIN_CAPTURE
+  {
+    int i,j,offset;
+
+    // error if tags are not defined
+    // NOTE: test that name = image or ids is not a 64-bit int in code?
+
+    int flag = 0;
+    if (lmp->atom->tag_enable == 0) flag = 1;
+    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
+    if (flag) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
+      return;
+    }
+
+    int natoms = static_cast<int> (lmp->atom->natoms);
+
+    void *vptr = lmp->atom->extract(name);
+    if (vptr == NULL) {
+      lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name");
+      return;
+    }
+
+    // perform MPI_Allgatherv on each proc's chunk of Nlocal atoms
+    
+    int nprocs = lmp->comm->nprocs;
+
+    int *recvcounts,*displs;
+    lmp->memory->create(recvcounts,nprocs,"lib/gather:recvcounts");
+    lmp->memory->create(displs,nprocs,"lib/gather:displs");
+
+    if (type == 0) {
+      int *vector = NULL;
+      int **array = NULL;
+      const int imgunpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgunpack) vector = (int *) vptr;
+      else array = (int **) vptr;
+
+      int *copy;
+      lmp->memory->create(copy,count*natoms,"lib/gather:copy");
+      for (i = 0; i < count*natoms; i++) copy[i] = 0;
+
+      tagint *tag = lmp->atom->tag;
+      int nlocal = lmp->atom->nlocal;
+
+      if (count == 1) {
+        MPI_Allgather(&nlocal,1,MPI_INT,recvcounts,1,MPI_INT,lmp->world);
+        displs[0] = 0;
+        for (i = 1; i < nprocs; i++)
+          displs[i] = displs[i-1] + recvcounts[i-1];
+        MPI_Allgatherv(vector,nlocal,MPI_INT,data,recvcounts,displs,
+                       MPI_INT,lmp->world);
+
+      } else if (imgunpack) {
+        int *copy;
+        lmp->memory->create(copy,count*nlocal,"lib/gather:copy");
+        offset = 0;
+        for (i = 0; i < nlocal; i++) {
+          const int image = vector[i];
+          copy[offset++] = (image & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX;
+        }
+        int n = count*nlocal;
+        MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,lmp->world);
+        displs[0] = 0;
+        for (i = 1; i < nprocs; i++)
+          displs[i] = displs[i-1] + recvcounts[i-1];
+        MPI_Allgatherv(copy,count*nlocal,MPI_INT,
+                       data,recvcounts,displs,MPI_INT,lmp->world);
+        lmp->memory->destroy(copy);
+
+      } else {
+        int n = count*nlocal;
+        MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,lmp->world);
+        displs[0] = 0;
+        for (i = 1; i < nprocs; i++)
+          displs[i] = displs[i-1] + recvcounts[i-1];
+        MPI_Allgatherv(&array[0][0],count*nlocal,MPI_INT,
+                       data,recvcounts,displs,MPI_INT,lmp->world);
+      }
+
+    } else {
+      double *vector = NULL;
+      double **array = NULL;
+      if (count == 1) vector = (double *) vptr;
+      else array = (double **) vptr;
+
+      int nlocal = lmp->atom->nlocal;
+
+      if (count == 1) {
+        MPI_Allgather(&nlocal,1,MPI_INT,recvcounts,1,MPI_INT,lmp->world);
+        displs[0] = 0;
+        for (i = 1; i < nprocs; i++)
+          displs[i] = displs[i-1] + recvcounts[i-1];
+        MPI_Allgatherv(vector,nlocal,MPI_DOUBLE,data,recvcounts,displs,
+                       MPI_DOUBLE,lmp->world);
+
+      } else {
+        int n = count*nlocal;
+        MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,lmp->world);
+        displs[0] = 0;
+        for (i = 1; i < nprocs; i++)
+          displs[i] = displs[i-1] + recvcounts[i-1];
+        MPI_Allgatherv(&array[0][0],count*nlocal,MPI_DOUBLE,
+                       data,recvcounts,displs,MPI_DOUBLE,lmp->world);
+      }
+    }
+
+    lmp->memory->destroy(recvcounts);
+    lmp->memory->destroy(displs);
+  }
+  END_CAPTURE
+}
+
+/* ----------------------------------------------------------------------
+   gather the named atom-based entity for a subset of atoms
+     return it in user-allocated data
+   data will be ordered by requested atom IDs
+     no requirement for consecutive atom IDs (1 to N)
+   see gather_atoms() to return data for all atoms, ordered by consecutive IDs
+   see gather_atoms_concat() to return data for all atoms, unordered
+   name = desired quantity, e.g. x or charge
+   type = 0 for integer values, 1 for double values
+   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+     use count = 3 with "image" if want single image flag unpacked into xyz
+   ndata = # of atoms to return data for (could be all atoms)
+   ids = list of Ndata atom IDs to return data for
+   return atom-based values in 1d data, ordered by count, then by atom
+     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     data must be pre-allocated by caller to correct length
+     correct length = count*Ndata
+   method:
+     alloc and zero count*Ndata length vector
+     loop over Ndata to fill vector with my values
+     Allreduce to sum vector into data across all procs
+------------------------------------------------------------------------- */
+
+void lammps_gather_atoms_subset(void *ptr, char *name,
+                                int type, int count,
+                                int ndata, int *ids, void *data)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  BEGIN_CAPTURE
+  { 
+    int i,j,m,offset;
+    tagint id;
+
+    // error if tags are not defined
+    // NOTE: test that name = image or ids is not a 64-bit int in code?
+
+    int flag = 0;
+    if (lmp->atom->tag_enable == 0) flag = 1;
+    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
+    if (flag) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms_subset");
+      return;
+    }
+
+    void *vptr = lmp->atom->extract(name);
+    if (vptr == NULL) {
+      lmp->error->warning(FLERR,"lammps_gather_atoms_subset: "
+                          "unknown property name");
+      return;
+    }
+
+    // copy = Ndata length vector of per-atom values
+    // use atom ID to insert each atom's values into copy
+    // MPI_Allreduce with MPI_SUM to merge into data
+
+    if (type == 0) {
+      int *vector = NULL;
+      int **array = NULL;
+      const int imgunpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgunpack) vector = (int *) vptr;
+      else array = (int **) vptr;
+
+      int *copy;
+      lmp->memory->create(copy,count*ndata,"lib/gather:copy");
+      for (i = 0; i < count*ndata; i++) copy[i] = 0;
+
+      tagint *tag = lmp->atom->tag;
+      int nlocal = lmp->atom->nlocal;
+
+      if (count == 1) {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0 && m < nlocal)
+            copy[i] = vector[m];
+        }
+
+      } else if (imgunpack) {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0 && m < nlocal) {
+            offset = count*i;
+            const int image = vector[m];
+            copy[offset++] = (image & IMGMASK) - IMGMAX;
+            copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX;
+            copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX;
+          }
+        }
+
+      } else {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0 && m < nlocal) {
+            offset = count*i;
+            for (j = 0; j < count; j++)
+              copy[offset++] = array[m][j];
+          }
+        }
+      }
+
+      MPI_Allreduce(copy,data,count*ndata,MPI_INT,MPI_SUM,lmp->world);
+      lmp->memory->destroy(copy);
+
+    } else {
+      double *vector = NULL;
+      double **array = NULL;
+      if (count == 1) vector = (double *) vptr;
+      else array = (double **) vptr;
+
+      double *copy;
+      lmp->memory->create(copy,count*ndata,"lib/gather:copy");
+      for (i = 0; i < count*ndata; i++) copy[i] = 0.0;
+
+      tagint *tag = lmp->atom->tag;
+      int nlocal = lmp->atom->nlocal;
+
+      if (count == 1) {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0 && m < nlocal)
+            copy[i] = vector[m];
+        }
+
+      } else {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0 && m < nlocal) {
+            offset = count*i;
+            for (j = 0; j < count; j++)
+              copy[offset++] = array[m][j];
+          }
+        }
+      }
+
+      MPI_Allreduce(copy,data,count*ndata,MPI_DOUBLE,MPI_SUM,lmp->world);
+      lmp->memory->destroy(copy);
+    }
+  }
+  END_CAPTURE
+}
+
+/* ----------------------------------------------------------------------
+   scatter the named atom-based entity in data to all atoms
+   data is ordered by atom ID
+     requirement for consecutive atom IDs (1 to N)
+   see scatter_atoms_subset() to scatter data for some (or all) atoms, unordered
+   name = desired quantity, e.g. x or charge
+   type = 0 for integer values, 1 for double values
+   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+     use count = 3 with "image" for xyz to be packed into single image flag
    data = atom-based values in 1d data, ordered by count, then by atom ID
      e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     data must be correct length = count*Natoms, as queried by get_natoms()
+   method:
+     loop over Natoms, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
 
 void lammps_scatter_atoms(void *ptr, char *name,
@@ -899,7 +1206,10 @@ void lammps_scatter_atoms(void *ptr, char *name,
 
   BEGIN_CAPTURE
   {
+    int i,j,m,offset;
+
     // error if tags are not defined or not consecutive or no atom map
+    // NOTE: test that name = image or ids is not a 64-bit int in code?
 
     int flag = 0;
     if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0)
@@ -914,7 +1224,6 @@ void lammps_scatter_atoms(void *ptr, char *name,
 
     int natoms = static_cast<int> (lmp->atom->natoms);
 
-    int i,j,m,offset;
     void *vptr = lmp->atom->extract(name);
     if(vptr == NULL) {
         lmp->error->warning(FLERR,
@@ -939,6 +1248,7 @@ void lammps_scatter_atoms(void *ptr, char *name,
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0)
             vector[m] = dptr[i];
+
       } else if (imgpack) {
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0) {
@@ -948,6 +1258,7 @@ void lammps_scatter_atoms(void *ptr, char *name,
             image += (dptr[offset++] + IMGMAX) << IMG2BITS;
             vector[m] = image;
           }
+
       } else {
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0) {
@@ -956,6 +1267,7 @@ void lammps_scatter_atoms(void *ptr, char *name,
               array[m][j] = dptr[offset++];
           }
       }
+
     } else {
       double *vector = NULL;
       double **array = NULL;
@@ -967,6 +1279,7 @@ void lammps_scatter_atoms(void *ptr, char *name,
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0)
             vector[m] = dptr[i];
+
       } else {
         for (i = 0; i < natoms; i++) {
           if ((m = lmp->atom->map(i+1)) >= 0) {
@@ -981,6 +1294,125 @@ void lammps_scatter_atoms(void *ptr, char *name,
   END_CAPTURE
 }
 
+/* ----------------------------------------------------------------------
+   scatter the named atom-based entity in data to a subset of atoms
+   data is ordered by provided atom IDs
+     no requirement for consecutive atom IDs (1 to N)
+   see scatter_atoms() to scatter data for all atoms, ordered by consecutive IDs
+   name = desired quantity, e.g. x or charge
+   type = 0 for integer values, 1 for double values
+   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+     use count = 3 with "image" for xyz to be packed into single image flag
+   ndata = # of atoms in ids and data (could be all atoms)
+   ids = list of Ndata atom IDs to scatter data to
+   data = atom-based values in 1d data, ordered by count, then by atom ID
+     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     data must be correct length = count*Ndata
+   method:
+     loop over Ndata, if I own atom ID, set its values from data
+------------------------------------------------------------------------- */
+
+void lammps_scatter_atoms_subset(void *ptr, char *name,
+                                 int type, int count, 
+                                 int ndata, int *ids, void *data)
+{ 
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  BEGIN_CAPTURE
+  {
+    int i,j,m,offset;
+    tagint id;
+
+    // error if tags are not defined or no atom map
+    // NOTE: test that name = image or ids is not a 64-bit int in code?
+
+    int flag = 0;
+    if (lmp->atom->tag_enable == 0) flag = 1;
+    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
+    if (lmp->atom->map_style == 0) flag = 1;
+    if (flag) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms_subset");
+      return;
+    }
+
+    void *vptr = lmp->atom->extract(name);
+    if(vptr == NULL) {
+        lmp->error->warning(FLERR,
+                            "lammps_scatter_atoms_subset: unknown property name");
+        return;
+    }
+
+    // copy = Natom length vector of per-atom values
+    // use atom ID to insert each atom's values into copy
+    // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
+
+    if (type == 0) {
+      int *vector = NULL;
+      int **array = NULL;
+      const int imgpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgpack) vector = (int *) vptr;
+      else array = (int **) vptr;
+      int *dptr = (int *) data;
+
+      if (count == 1) {
+        for (i = 0; i < ndata; i++)
+          if ((m = lmp->atom->map(i+1)) >= 0)
+            vector[m] = dptr[i];
+
+      } else if (imgpack) {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0) {
+            offset = count*i;
+            int image = dptr[offset++] + IMGMAX;
+            image += (dptr[offset++] + IMGMAX) << IMGBITS;
+            image += (dptr[offset++] + IMGMAX) << IMG2BITS;
+            vector[m] = image;
+          }
+        }
+
+      } else {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0) {
+            offset = count*i;
+            for (j = 0; j < count; j++)
+              array[m][j] = dptr[offset++];
+          }
+        }
+      }
+
+    } else {
+      double *vector = NULL;
+      double **array = NULL;
+      if (count == 1) vector = (double *) vptr;
+      else array = (double **) vptr;
+      double *dptr = (double *) data;
+
+      if (count == 1) {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0)
+            vector[m] = dptr[i];
+        }
+
+      } else {
+        for (i = 0; i < ndata; i++) {
+          id = ids[i];
+          if ((m = lmp->atom->map(id)) >= 0) {
+            offset = count*i;
+            for (j = 0; j < count; j++)
+              array[m][j] = dptr[offset++];
+          }
+        }
+      }
+    }
+  }
+  END_CAPTURE
+}
+
 /* ----------------------------------------------------------------------
    create N atoms and assign them to procs based on coords
    id = atom IDs (optional, NULL will generate 1 to N)
diff --git a/src/library.h b/src/library.h
index 1cb3f9cbf38233c6437fe93bf9f3114879c1d731..93d9f0ef382036c9e3d0c5e4c36a697306ae26b6 100644
--- a/src/library.h
+++ b/src/library.h
@@ -43,13 +43,17 @@ void *lammps_extract_compute(void *, char *, int, int);
 void *lammps_extract_fix(void *, char *, int, int, int, int);
 void *lammps_extract_variable(void *, char *, char *);
 
-void lammps_reset_box(void *, double *, double *, double, double, double);
-int lammps_set_variable(void *, char *, char *);
 double lammps_get_thermo(void *, char *);
-
 int lammps_get_natoms(void *);
+
+int lammps_set_variable(void *, char *, char *);
+void lammps_reset_box(void *, double *, double *, double, double, double);
+
 void lammps_gather_atoms(void *, char *, int, int, void *);
+void lammps_gather_atoms_concat(void *, char *, int, int, void *);
+void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
 void lammps_scatter_atoms(void *, char *, int, int, void *);
+void lammps_scatter_atoms_subset(void *, char *, int, int, int, int *, void *);
 
 // lammps_create_atoms() takes tagint and imageint as args
 // ifdef insures they are compatible with rest of LAMMPS