diff --git a/python/examples/matplotlib_plot.py b/python/examples/matplotlib_plot.py
new file mode 100755
index 0000000000000000000000000000000000000000..d4e304aa19f193e628369540b7ccc59cb5f86fba
--- /dev/null
+++ b/python/examples/matplotlib_plot.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python -i
+# preceding line should have path for Python on your machine
+
+# matplotlib_plot.py
+# Purpose: plot Temp of running LAMMPS simulation via matplotlib
+# Syntax:  plot.py in.lammps Nfreq Nsteps compute-ID
+#          in.lammps = LAMMPS input script
+#          Nfreq = plot data point every this many steps
+#          Nsteps = run for this many steps
+#          compute-ID = ID of compute that calculates temperature
+#                       (or any other scalar quantity)
+
+from __future__ import print_function
+import sys
+sys.path.append("./pizza")
+import matplotlib
+matplotlib.use('tkagg')
+import matplotlib.pyplot as plt
+
+# parse command line
+
+argv = sys.argv
+if len(argv) != 5:
+  print("Syntax: plot.py in.lammps Nfreq Nsteps compute-ID")
+  sys.exit()
+
+infile = sys.argv[1]
+nfreq = int(sys.argv[2])
+nsteps = int(sys.argv[3])
+compute = sys.argv[4]
+
+me = 0
+# uncomment if running in parallel via Pypar
+#import pypar
+#me = pypar.rank()
+#nprocs = pypar.size()
+
+from lammps import lammps
+lmp = lammps()
+
+# run infile all at once
+# assumed to have no run command in it
+
+lmp.file(infile)
+lmp.command("thermo %d" % nfreq)
+
+# initial 0-step run to generate initial 1-point plot
+
+lmp.command("run 0 pre yes post no")
+value = lmp.extract_compute(compute,0,0)
+ntimestep = 0
+xaxis = [ntimestep]
+yaxis = [value]
+
+# create matplotlib plot
+# just proc 0 handles plotting
+
+if me == 0:
+  fig = plt.figure()
+  line, = plt.plot(xaxis, yaxis)
+  plt.xlim([0, nsteps])
+  plt.title(compute)
+  plt.xlabel("Timestep")
+  plt.ylabel("Temperature")
+  plt.show(block=False)
+
+# run nfreq steps at a time w/out pre/post, query compute, refresh plot
+import time
+
+while ntimestep < nsteps:
+  lmp.command("run %d pre no post no" % nfreq)
+  ntimestep += nfreq
+  value = lmp.extract_compute(compute,0,0)
+  xaxis.append(ntimestep)
+  yaxis.append(value)
+  if me == 0:
+    line.set_xdata(xaxis)
+    line.set_ydata(yaxis)
+    ax = plt.gca()
+    ax.relim()
+    ax.autoscale_view(True, True, True)
+    fig.canvas.draw()
+
+lmp.command("run 0 pre no post yes")
+
+# uncomment if running in parallel via Pypar
+#print("Proc %d out of %d procs has" % (me,nprocs), lmp)
+#pypar.finalize()
+
+if sys.version_info[0] == 3:
+    input("Press Enter to exit...")
+else:
+    raw_input("Press Enter to exit...")
diff --git a/python/examples/pizza/dump.py b/python/examples/pizza/dump.py
index ec18cba62eefee12b189f1e7f4d2f64527dc5857..14645ac1f3ad13bb33b8b515d58b969238677f27 100644
--- a/python/examples/pizza/dump.py
+++ b/python/examples/pizza/dump.py
@@ -146,6 +146,8 @@ d.extra(data)				   extract bond/tri/line list from data
 # History
 #   8/05, Steve Plimpton (SNL): original version
 #   12/09, David Hart (SNL): allow use of NumPy or Numeric
+#   03/17, Richard Berger (Temple U): improve Python 3 compatibility,
+#                                     simplify read_snapshot by using reshape
 
 # ToDo list
 #   try to optimize this line in read_snap: words += f.readline().split()
@@ -224,7 +226,7 @@ class dump:
     self.flist = []
     for word in words: self.flist += glob.glob(word)
     if len(self.flist) == 0 and len(list) == 1:
-      raise StandardError("no dump file specified")
+      raise Exception("no dump file specified")
     
     if len(list) == 1:
       self.increment = 0
@@ -299,7 +301,7 @@ class dump:
 
   def next(self):
 
-    if not self.increment: raise StandardError("cannot read incrementally")
+    if not self.increment: raise Exception("cannot read incrementally")
 
     # read next snapshot in current file using eof as pointer
     # if fail, try next file
@@ -344,13 +346,13 @@ class dump:
     try:
       snap = Snap()
       item = f.readline()
-      snap.time = int(f.readline().split()[0])    # just grab 1st field
+      snap.time = int(f.readline().decode().split()[0])    # just grab 1st field
       item = f.readline()
-      snap.natoms = int(f.readline())
+      snap.natoms = int(f.readline().decode())
 
       snap.aselect = np.zeros(snap.natoms)
 
-      item = f.readline()
+      item = f.readline().decode()
       words = f.readline().split()
       snap.xlo,snap.xhi = float(words[0]),float(words[1])
       words = f.readline().split()
@@ -358,7 +360,7 @@ class dump:
       words = f.readline().split()
       snap.zlo,snap.zhi = float(words[0]),float(words[1])
 
-      item = f.readline()
+      item = f.readline().decode()
       if len(self.names) == 0:
         words = item.split()[2:]
         if len(words):
@@ -372,24 +374,22 @@ class dump:
             else: self.names[words[i]] = i
 
       if snap.natoms:
-        words = f.readline().split()
+        words = f.readline().decode().split()
         ncol = len(words)
         for i in range(1,snap.natoms):
-          words += f.readline().split()
+          words += f.readline().decode().split()
         floats = map(float,words)
-        if oldnumeric: atoms = np.zeros((snap.natoms,ncol),np.Float)
-        else: atoms = np.zeros((snap.natoms,ncol),np.float)
-        start = 0
-        stop = ncol
-        for i in range(snap.natoms):
-          atoms[i] = floats[start:stop]
-          start = stop
-          stop += ncol
-      else: atoms = None
-      snap.atoms = atoms
+        if oldnumeric:
+          atom_data = np.array(list(floats),np.Float)
+        else:
+          atom_data = np.array(list(floats),np.float)
+
+        snap.atoms = atom_data.reshape((snap.natoms, ncol))
+      else:
+        snap.atoms = None
       return snap
     except:
-      return 0
+      return None
 
   # --------------------------------------------------------------------
   # decide if snapshot i is scaled/unscaled from coords of first and last atom
@@ -417,7 +417,7 @@ class dump:
   
   def map(self,*pairs):
     if len(pairs) % 2 != 0:
-      raise StandardError("dump map() requires pairs of mappings")
+      raise Exception("dump map() requires pairs of mappings")
     for i in range(0,len(pairs),2):
       j = i + 1
       self.names[pairs[j]] = pairs[i]-1
@@ -734,7 +734,7 @@ class dump:
     for snap in self.snaps:
       if not snap.tselect: continue
       if snap.nselect != len(vec):
-        raise StandardError("vec length does not match # of selected atoms")
+        raise Exception("vec length does not match # of selected atoms")
       atoms = snap.atoms
       m = 0
       for i in range(snap.natoms):
@@ -800,7 +800,7 @@ class dump:
 
   def atom(self,n,*list):
     if len(list) == 0:
-      raise StandardError("no columns specified")
+      raise Exception("no columns specified")
     columns = []
     values = []
     for name in list:
@@ -816,7 +816,7 @@ class dump:
       for i in range(snap.natoms):
         if atoms[i][id] == n: break
       if atoms[i][id] != n:
-        raise StandardError("could not find atom ID in snapshot")
+        raise Exception("could not find atom ID in snapshot")
       for j in range(ncol):
         values[j][m] = atoms[i][columns[j]]
       m += 1
@@ -831,7 +831,7 @@ class dump:
     snap = self.snaps[self.findtime(n)]
     
     if len(list) == 0:
-      raise StandardError("no columns specified")
+      raise Exception("no columns specified")
     columns = []
     values = []
     for name in list:
@@ -957,9 +957,9 @@ class dump:
   # --------------------------------------------------------------------
 
   def findtime(self,n):
-    for i in range(self.nsnaps):
-      if self.snaps[i].time == n: return i
-    raise StandardError("no step %d exists" % n)
+    for i, snap in enumerate(self.snaps):
+      if snap.time == n: return i
+    raise Exception("no step %d exists" % n)
 
   # --------------------------------------------------------------------
   # return maximum box size across all selected snapshots
@@ -1008,7 +1008,7 @@ class dump:
         nbonds = int(f.readline())
         item = f.readline()
         if not re.search("BONDS",item):
-          raise StandardError("could not read bonds from dump file")
+          raise Exception("could not read bonds from dump file")
 
         words = f.readline().split()
         ncol = len(words)
@@ -1031,7 +1031,7 @@ class dump:
           self.bondflag = 1
           self.bondlist = bondlist
       except:
-        raise StandardError("could not read from bond dump file")
+        raise Exception("could not read from bond dump file")
       
     # request bonds from data object
     
@@ -1047,7 +1047,7 @@ class dump:
           self.bondflag = 1
           self.bondlist = bondlist
       except:
-        raise StandardError("could not extract bonds from data object")
+        raise Exception("could not extract bonds from data object")
 
     # request tris/lines from cdata object
     
@@ -1061,7 +1061,7 @@ class dump:
           self.lineflag = 1
           self.linelist = lines
       except:
-        raise StandardError("could not extract tris/lines from cdata object")
+        raise Exception("could not extract tris/lines from cdata object")
 
     # request tris from mdump object
     
@@ -1070,10 +1070,10 @@ class dump:
         self.triflag = 2
         self.triobj = arg
       except:
-        raise StandardError("could not extract tris from mdump object")
+        raise Exception("could not extract tris from mdump object")
 
     else:
-      raise StandardError("unrecognized argument to dump.extra()")
+      raise Exception("unrecognized argument to dump.extra()")
       
   # --------------------------------------------------------------------
 
diff --git a/python/examples/pizza/gl.py b/python/examples/pizza/gl.py
index 765653c9836d7be4cde8e3fdc9c5a76f685b1e69..908db67052e1db71ef0f4fbfcfea2c8f2cc9a955 100644
--- a/python/examples/pizza/gl.py
+++ b/python/examples/pizza/gl.py
@@ -638,7 +638,7 @@ class gl:
                           fraction*(self.scale_stop - self.scale_start)
           self.viewupright()
 
-	if n == nstart or self.panflag: self.center = compute_center(box)
+        if n == nstart or self.panflag: self.center = compute_center(box)
 
         if bonds: self.bonds_augment(bonds)
 
diff --git a/python/examples/pizza/pdbfile.py b/python/examples/pizza/pdbfile.py
index 2ca85b64da81ec4afcf183cfeabc91490cdcc3bc..51c88604249491dae6f5772be1be949f71e801fe 100644
--- a/python/examples/pizza/pdbfile.py
+++ b/python/examples/pizza/pdbfile.py
@@ -54,6 +54,7 @@ index,time,flag = p.iterator(1)
 
 # History
 #   8/05, Steve Plimpton (SNL): original version
+#   3/17, Richard Berger (Temple U): improve Python 3 compatibility
 
 # ToDo list
 #   for generic PDB file (no template) from a LJ unit system,
@@ -68,6 +69,12 @@ index,time,flag = p.iterator(1)
 # Imports and external programs
 
 import sys, types, glob, urllib
+PY3 = sys.version_info[0] == 3
+
+if PY3:
+    string_types = str,
+else:
+    string_types = basestring
 
 # Class definition
 
@@ -77,7 +84,7 @@ class pdbfile:
 
   def __init__(self,*args):
     if len(args) == 1:
-      if type(args[0]) is types.StringType:
+      if type(args[0]) is string_types:
         filestr = args[0]
         self.data = None
       else:
diff --git a/python/examples/viz_atomeye.py b/python/examples/viz_atomeye.py
index d9cf4dcef274f3675d8891d6f64608893d0165a4..d7f83f274dfd57f985322d8f3c05db7d9506f779 100755
--- a/python/examples/viz_atomeye.py
+++ b/python/examples/viz_atomeye.py
@@ -42,7 +42,7 @@ lmp = lammps()
 
 lmp.file(infile)
 lmp.command("thermo %d" % nfreq)
-lmp.command("dump python all cfg %d tmp.cfg.* id type xs ys zs" % nfreq)
+lmp.command("dump python all cfg %d tmp.cfg.* mass type xs ys zs id" % nfreq)
 
 # initial 0-step run to generate dump file and image
 
diff --git a/python/examples/vizplotgui_atomeye.py b/python/examples/vizplotgui_atomeye.py
index 02f119b7260c0411316f13270c65d53aa642b306..8f5f4e3dfff638264cf8dd17ce5f07edad185528 100755
--- a/python/examples/vizplotgui_atomeye.py
+++ b/python/examples/vizplotgui_atomeye.py
@@ -73,7 +73,7 @@ lmp = lammps()
 
 lmp.file(infile)
 lmp.command("thermo %d" % nfreq)
-lmp.command("dump python all cfg %d tmp.cfg.* id type xs ys zs" % nfreq)
+lmp.command("dump python all cfg %d tmp.cfg.* mass type xs ys zs id" % nfreq)
 
 # initial 0-step run to generate initial 1-point plot, dump file, and image
 
diff --git a/python/lammps.py b/python/lammps.py
index 7344a16be30fb8ab47e3ffc1f6cd939dae95b4fe..a36abb87e8fd408f953f58f0f5b14b2cb9f38e55 100644
--- a/python/lammps.py
+++ b/python/lammps.py
@@ -190,12 +190,15 @@ class lammps(object):
   # send a list of commands
 
   def commands_list(self,cmdlist):
-    args = (c_char_p * len(cmdlist))(*cmdlist)
+    cmds = [x.encode() for x in cmdlist if type(x) is str]
+    args = (c_char_p * len(cmdlist))(*cmds)
     self.lib.lammps_commands_list(self.lmp,len(cmdlist),args)
     
   # send a string of commands
 
   def commands_string(self,multicmd):
+    if type(multicmd) is str:
+        multicmd = multicmd.encode()
     self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd))
     
   # extract global info
@@ -359,19 +362,27 @@ class lammps(object):
   #   e.g. for Python list or NumPy, etc
   #   ditto for gather_atoms() above
 
-  def create_atoms(self,n,id,type,x,v):
+  def create_atoms(self,n,id,type,x,v,image=None,shrinkexceed=False):
     if id:
       id_lmp = (c_int * n)()
       id_lmp[:] = id
-    else: id_lmp = id
+    else:
+      id_lmp = id
+
+    if image:
+      image_lmp = (c_int * n)()
+      image_lmp[:] = image
+    else:
+      image_lmp = image
+
     type_lmp = (c_int * n)()
     type_lmp[:] = type
-    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v)
+    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v,image_lmp,shrinkexceed)
+
 
-  # document this?
-    
   @property
   def uses_exceptions(self):
+    """ Return whether the LAMMPS shared library was compiled with C++ exceptions handling enabled """
     try:
       if self.lib.lammps_has_error:
         return True
diff --git a/src/library.cpp b/src/library.cpp
index 5980b135eb1bd8ee00f417b662e63339a9f9b0c2..da491f71528a319a69a2e28f822ac72e6ad14399 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -971,11 +971,9 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
       xdata[0] = x[3*i];
       xdata[1] = x[3*i+1];
       xdata[2] = x[3*i+2];
-      if (image) {
-        if (!domain->ownatom(id[i],xdata,&image[i],shrinkexceed)) continue;
-      } else {
-        if (!domain->ownatom(id[i],xdata,NULL,shrinkexceed)) continue;
-      }
+      imageint * img = image ? &image[i] : NULL;
+      tagint     tag = id    ? id[i]     : -1;
+      if (!domain->ownatom(tag, xdata, img, shrinkexceed)) continue;
   
       atom->avec->create_atom(type[i],xdata);
       if (id) atom->tag[nlocal] = id[i];