diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index 7b8ed2d9a6..33cbfa9586 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -1864,8 +1864,8 @@ void lammps_close(void *)
 int lammps_version(void *)
 void lammps_file(void *, char *)
 char *lammps_command(void *, char *)
-char *lammps_commands_list(void *, int, char **)
-char *lammps_commands_concatenated(void *, char *)
+void lammps_commands_list(void *, int, char **)
+void lammps_commands_string(void *, char *)
 void lammps_free(void *) :pre
 
 The lammps_open() function is used to initialize LAMMPS, passing in a
@@ -1901,51 +1901,66 @@ changes to the LAMMPS command syntax between versions. The returned
 LAMMPS version code is an integer (e.g. 2 Sep 2015 results in
 20150902) that grows with every new LAMMPS version.
 
-The lammps_file(), lammps_command(), lammps_commands_list(),
-lammps_commands_concatenated() functions are used to pass one or more
-commands to LAMMPS to execute as if they were coming from an input
-script.
+The lammps_file(), lammps_command(), lammps_commands_list(), and
+lammps_commands_string() functions are used to pass one or more
+commands to LAMMPS to execute, the same as if they were coming from an
+input script.
+
+Via these functions, the calling code can read or generate a series of
+LAMMPS commands one or multiple at a time and pass it thru the library
+interface to setup a problem and then run it in stages.  The caller
+can interleave the command function calls with operations it performs,
+calls to extract information from or set information within LAMMPS, or
+calls to another code's library.
 
 The lammps_file() function passes the filename of an input script.
 The lammps_command() function passes a single command as a string.
-The lammps_commands_multiple() function passes multiple commands in a
-char** list.  The lammps_commands_concatentaed() function passes
-multiple commands concatenated into one long string, separated by
-newline characters.  A single command can be spread across multiple
-lines, if the last printable character of all but the last line is
-"&", the same as if the lines appeared in an input script.
+The lammps_commands_list() function passes multiple commands in a
+char** list.  In both lammps_command() and lammps_commands_list(),
+individual commands may or may not have a trailing newline.  The
+lammps_commands_string() function passes multiple commands
+concatenated into one long string, separated by newline characters.
+In both lammps_commands_list() and lammps_commands_string(), a single
+command can be spread across multiple lines, if the last printable
+character of all but the last line is "&", the same as if the lines
+appeared in an input script.
 
-Via these library functions, the calling code can read or generate a
-series of LAMMPS commands one or multiple at a time and pass it thru
-the library interface to setup a problem and then run it, interleaving
-the command function calls with operations performed within the
-faller, calls to extract information from LAMMPS, or calls to another
-code's library.
-
-The lammps_free() is a clean-up function to free memory that the library
-allocated previously.
+The lammps_free() function is a clean-up function to free memory that
+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:
+documentation in the src/library.cpp file for details, including
+which quantities can be queried by name:
 
 void *lammps_extract_global(void *, char *)
 void *lammps_extract_atom(void *, char *)
 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_extract_variable(void *, char *, char *) :pre
+
 int lammps_set_variable(void *, char *, char *)
-double lammps_get_thermo(void *, char *)
+double lammps_get_thermo(void *, char *) :pre
+
 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 *) :pre
 
-These functions can extract various global or per-atom quantities from
-LAMMPS as well as values calculated by a compute, fix, or variable.
-The lammps_set_variable() function can set an existing string-style variable
-to a new value, so that subsequent LAMMPS commands can access the
-variable.  The lammps_get_thermo() function returns the current
-value of a thermo keyword.
+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 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).
+
+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.
 
 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
@@ -1956,15 +1971,22 @@ returns a full list to each calling processor.  The scatter function
 does the inverse.  It distributes the same kinds of values, 
 passed by the caller, to each atom owned by individual processors.
 
+The lammps_create_atoms() function takes a list of N atoms as input
+with atom types and coords (required), an optionally atom IDs and
+velocities.  It uses the coords of each atom to assign it as a new
+atom to the processor that owns it.  Additional properties for the new
+atoms can 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 your own additional functions as needed to define
+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 routines you add can access or
+interface"_Section_python.html.  The added functions can access or
 change any LAMMPS data you wish.
 
 :line
diff --git a/examples/COUPLE/simple/simple.c b/examples/COUPLE/simple/simple.c
index c50f289b7a..fd78ae5dcc 100644
--- a/examples/COUPLE/simple/simple.c
+++ b/examples/COUPLE/simple/simple.c
@@ -71,8 +71,8 @@ int main(int narg, char **arg)
      (could just send it to proc 0 of comm_lammps and let it Bcast)
      all LAMMPS procs call lammps_command() on the line */
 
-  void *ptr;
-  if (lammps == 1) lammps_open(0,NULL,comm_lammps,&ptr);
+  void *lmp = NULL;
+  if (lammps == 1) lammps_open(0,NULL,comm_lammps,&lmp);
 
   int n;
   char line[1024];
@@ -85,7 +85,7 @@ int main(int narg, char **arg)
     MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
     if (n == 0) break;
     MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
-    if (lammps == 1) lammps_command(ptr,line);
+    if (lammps == 1) lammps_command(lmp,line);
   }
 
   /* run 10 more steps
@@ -94,23 +94,62 @@ int main(int narg, char **arg)
      put coords back into LAMMPS
      run a single step with changed coords */
 
-  if (lammps == 1) {
-    lammps_command(ptr,"run 10");
+  double *x = NULL;
+  double *v = NULL;
 
-    int natoms = lammps_get_natoms(ptr);
-    double *x = (double *) malloc(3*natoms*sizeof(double));
-    lammps_gather_atoms(ptr,"x",1,3,x);
+  if (lammps == 1) {
+    lammps_command(lmp,"run 10");
+
+    int natoms = lammps_get_natoms(lmp);
+    x = (double *) malloc(3*natoms*sizeof(double));
+    lammps_gather_atoms(lmp,"x",1,3,x);
+    v = (double *) malloc(3*natoms*sizeof(double));
+    lammps_gather_atoms(lmp,"v",1,3,v);
     double epsilon = 0.1;
     x[0] += epsilon;
-    lammps_scatter_atoms(ptr,"x",1,3,x);
-    free(x);
+    lammps_scatter_atoms(lmp,"x",1,3,x);
 
-    lammps_command(ptr,"run 1");
+    lammps_command(lmp,"run 1");
   }
 
-  if (lammps == 1) lammps_close(ptr);
+  /* use commands_list() and commands_string() to invoke more commands */
+
+  char *strtwo = "run 10\nrun 20";
+  if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+  char *cmds[2];
+  cmds[0] = "run 10";
+  cmds[1] = "run 20";
+  if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+  /* delete all atoms
+     create_atoms() to create new ones with old coords, vels
+     initial thermo should be same as step 20 */
+
+  int *type = NULL;
+
+  if (lammps == 1) {
+    int natoms = lammps_get_natoms(lmp);
+    type = (int *) malloc(natoms*sizeof(double));
+    int i;
+    for (i = 0; i < natoms; i++) type[i] = 1;
+
+    lammps_command(lmp,"delete_atoms group all");
+    lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+    lammps_command(lmp,"run 10");
+  }
+
+  if (x) free(x);
+  if (v) free(v);
+  if (type) free(type);
+
+  // close down LAMMPS
+
+  lammps_close(lmp);
 
   /* close down MPI */
 
+  if (lammps == 1) MPI_Comm_free(&comm_lammps);
+  MPI_Barrier(MPI_COMM_WORLD);
   MPI_Finalize();
 }
diff --git a/examples/COUPLE/simple/simple.cpp b/examples/COUPLE/simple/simple.cpp
index a440bcb1ae..e10a67d370 100644
--- a/examples/COUPLE/simple/simple.cpp
+++ b/examples/COUPLE/simple/simple.cpp
@@ -77,7 +77,7 @@ int main(int narg, char **arg)
   // (could just send it to proc 0 of comm_lammps and let it Bcast)
   // all LAMMPS procs call input->one() on the line
   
-  LAMMPS *lmp;
+  LAMMPS *lmp = NULL;
   if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);
 
   int n;
@@ -91,7 +91,7 @@ int main(int narg, char **arg)
     MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
     if (n == 0) break;
     MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
-    if (lammps == 1) lmp->input->one(line);
+    if (lammps == 1) lammps_command(lmp,line);
   }
 
   // run 10 more steps
@@ -100,23 +100,64 @@ int main(int narg, char **arg)
   // put coords back into LAMMPS
   // run a single step with changed coords
 
+  double *x = NULL;
+  double *v = NULL;
+
   if (lammps == 1) {
     lmp->input->one("run 10");
 
     int natoms = static_cast<int> (lmp->atom->natoms);
-    double *x = new double[3*natoms];
+    x = new double[3*natoms];
+    v = new double[3*natoms];
     lammps_gather_atoms(lmp,"x",1,3,x);
+    lammps_gather_atoms(lmp,"v",1,3,v);
     double epsilon = 0.1;
     x[0] += epsilon;
     lammps_scatter_atoms(lmp,"x",1,3,x);
-    delete [] x;
 
+    // these 2 lines are the same
+
+    // lammps_command(lmp,"run 1");
     lmp->input->one("run 1");
   }
 
-  if (lammps == 1) delete lmp;
+  // use commands_list() and commands_string() to invoke more commands
+
+  char *strtwo = "run 10\nrun 20";
+  if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+  char *cmds[2];
+  cmds[0] = "run 10";
+  cmds[1] = "run 20";
+  if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+  // delete all atoms
+  // create_atoms() to create new ones with old coords, vels
+  // initial thermo should be same as step 20
+
+  int *type = NULL;
+
+  if (lammps == 1) {
+    int natoms = static_cast<int> (lmp->atom->natoms);
+    type = new int[natoms];
+    for (int i = 0; i < natoms; i++) type[i] = 1;
+
+    lmp->input->one("delete_atoms group all");
+    lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+    lmp->input->one("run 10");
+  }
+
+  delete [] x;
+  delete [] v;
+  delete [] type;
+
+  // close down LAMMPS
+
+  delete lmp;
 
   // close down MPI
 
+  if (lammps == 1) MPI_Comm_free(&comm_lammps);
+  MPI_Barrier(MPI_COMM_WORLD);
   MPI_Finalize();
 }
diff --git a/python/examples/simple.py b/python/examples/simple.py
index 9315a555fc..1832ba6274 100755
--- a/python/examples/simple.py
+++ b/python/examples/simple.py
@@ -51,17 +51,39 @@ for line in lines: lmp.command(line)
 
 lmp.command("run 10")
 x = lmp.gather_atoms("x",1,3)
+v = lmp.gather_atoms("v",1,3)
 epsilon = 0.1
 x[0] += epsilon
 lmp.scatter_atoms("x",1,3,x)
 lmp.command("run 1");
 
+# extract force on single atom
+
 f = lmp.extract_atom("f",3)
+print("F",f)
 print("Force on 1 atom via extract_atom: ",f[0][0])
 
 fx = lmp.extract_variable("fx","all",1)
+print("FX",fx)
 print("Force on 1 atom via extract_variable:",fx[0])
 
+# use commands_list() to invoke more commands
+
+#strtwo = ["run 10","run 20"]
+#lmp.commands_list(strtwo)
+#lmp.commands_list(strtwo)
+
+# delete all atoms
+# create_atoms() to create new ones with old coords, vels
+# initial thermo should be same as step 20
+
+#natoms = lmp.get_natoms()
+#type = natoms*[1]
+
+#lmp.command("delete_atoms group all");
+#lmp.create_atoms(natoms,None,type,x,v);
+#lmp.command("run 10");
+
 # uncomment if running in parallel via Pypar
 #print("Proc %d out of %d procs has" % (me,nprocs), lmp)
 #pypar.finalize()
diff --git a/python/lammps.py b/python/lammps.py
index 8c57928064..f67938d036 100644
--- a/python/lammps.py
+++ b/python/lammps.py
@@ -172,6 +172,8 @@ class lammps(object):
     if file: file = file.encode()
     self.lib.lammps_file(self.lmp,file)
 
+  # send a single command
+    
   def command(self,cmd):
     if cmd: cmd = cmd.encode()
     self.lib.lammps_command(self.lmp,cmd)
@@ -185,6 +187,13 @@ class lammps(object):
         raise MPIAbortException(error_msg)
       raise Exception(error_msg)
 
+  # send a list of commands, one at a time
+
+  def commands_list(self,cmdlist):
+    for cmd in cmdlist: self.command(cmd)
+
+  # extract global info
+    
   def extract_global(self,name,type):
     if name: name = name.encode()
     if type == 0:
@@ -195,6 +204,8 @@ class lammps(object):
     ptr = self.lib.lammps_extract_global(self.lmp,name)
     return ptr[0]
 
+  # extract per-atom info
+  
   def extract_atom(self,name,type):
     if name: name = name.encode()
     if type == 0:
@@ -209,6 +220,8 @@ class lammps(object):
     ptr = self.lib.lammps_extract_atom(self.lmp,name)
     return ptr
 
+  # extract compute info
+  
   def extract_compute(self,id,style,type):
     if id: id = id.encode()
     if type == 0:
@@ -226,6 +239,7 @@ class lammps(object):
       return ptr
     return None
 
+  # extract fix info
   # in case of global datum, free memory for 1 double via lammps_free()
   # double was allocated by library interface function
 
@@ -249,6 +263,7 @@ class lammps(object):
     else:
       return None
 
+  # extract variable info
   # free memory for 1 double or 1 vector of doubles via lammps_free()
   # for vector, must copy nlocal returned values to local c_double vector
   # memory was allocated by library interface function
@@ -296,7 +311,14 @@ class lammps(object):
     return self.lib.lammps_get_natoms(self.lmp)
 
   # return vector of atom properties gathered across procs, ordered by atom ID
-
+  # 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 extact_atom() above
+  
   def gather_atoms(self,name,type,count):
     if name: name = name.encode()
     natoms = self.lib.lammps_get_natoms(self.lmp)
@@ -310,12 +332,38 @@ class lammps(object):
     return data
 
   # scatter vector of atom properties across procs, ordered by atom ID
-  # assume vector is of correct type and length, as created by gather_atoms()
-
+  # 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
+    
   def scatter_atoms(self,name,type,count,data):
     if name: name = name.encode()
     self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data)
 
+  # create N atoms on all procs
+  # N = global number of atoms
+  # id = ID of each atom (optional, can be None)
+  # type = type of each atom (1 to Ntypes) (required)
+  # x = coords of each atom as (N,3) array (required)
+  # v = velocity of each atom as (N,3) array (optional, can be None)
+  # NOTE: how could we insure are passing correct type to LAMMPS
+  #   e.g. for Python list or NumPy, etc
+  #   ditto for gather_atoms() above
+
+  def create_atoms(self,n,id,type,x,v):
+    if id:
+      id_lmp = (c_int * n)()
+      id_lmp[:] = id
+    else: id_lmp = id
+    type_lmp = (c_int * n)()
+    type_lmp[:] = type
+    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v)
+
+  # document this?
+    
   @property
   def uses_exceptions(self):
     try: