diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst
index b6f65a9525..54adb7d721 100644
--- a/doc/src/pair_mliap.rst
+++ b/doc/src/pair_mliap.rst
@@ -10,7 +10,7 @@ Syntax
 
    pair_style mliap ... keyword values ...
 
-* two keyword/value pairs must be appended
+* one or two keyword/value pairs must be appended
 * keyword = *model* or *descriptor* or *unified*
 
   .. parsed-literal::
@@ -122,10 +122,10 @@ The detail of *nn* module implementation can be found at :ref:`(Yanxon) <Yanxon2
    will be loaded using pytorch; otherwise, it will be loaded as a pickle file.To load
    a model from memory (i.e. an existing python object), specify the filename as
    "LATER", and then call `lammps.mliap.load_model(model)` from python
-   before using the pair style. When using lammps via the library mode, you
+   before using the pair style. When using LAMMPS via the library mode, you
    will need to call `lammps.mliappy.activate_mliappy(lmp)` on the active
-   lammps object before the pair style is defined. This call locates and
-   loads the mliap-specific python module that is built into lammps.
+   LAMMPS object before the pair style is defined. This call locates and
+   loads the mliap-specific python module that is built into LAMMPS.
 
 The *descriptor* keyword is followed by a descriptor style, and additional arguments.
 Currently two descriptor styles are available: *sna* and *so3*.
diff --git a/examples/mliap/README b/examples/mliap/README
index c16fa218e0..20792d2c2c 100644
--- a/examples/mliap/README
+++ b/examples/mliap/README
@@ -109,4 +109,70 @@ in.mliap.nn.cu
 -------------------------
 Run a neural network potential for Cu, a combination of SNAP descriptors and the "nn" model style
 
+in.mliap.unified.lj.Ar
+-----------------------
+This reproduces the output of examples/melt/in.melt,
+but using the ``unified`` keyword and 
+the Python coupling to PyTorch.
 
+This example can be run in two different ways:
+
+1: Running a LAMMPS executable: in.mliap.unified.lj.Ar
+
+First run ``python ./pickle_mliap_unified_lj_Ar.py``. It creates
+an MLIAP unified model that replicates the ``lj/cut`` pair
+style defined in examples/melt/in/melt
+and saves it in the file "mliap_unified_lj_Ar.pkl".
+
+You can then run the example as follows
+
+`lmp -in in.mliap.unified.lj.Ar -echo both`
+
+The resultant log.lammps output should be identical to that generated
+by in.melt.
+
+If this fails, see the instructions for building the MLIAP package
+with Python support enabled. Also, confirm that the
+LAMMPS Python embedded Python interpreter is
+working by running ../examples/in.python.
+
+2: Running a Python script: mliap_unified_lj_Ar.py
+
+Before testing this, ensure that the previous method
+(running a LAMMPS executable) works.
+
+You can run the example in serial: 
+
+`python mliap_unified_lj_Ar.py`
+
+or in parallel:
+
+`mpirun -np 4 python mliap_unified_lj_Ar.py`
+
+The resultant log.lammps output should be identical to that generated
+by in.melt and in.mliap.unified.lj.Ar.
+
+in.mliap.unified.hippynn.Al
+---------------------------
+
+Demonstrate MLIAP interface to HIPNN Al potential. Requires the HIPNN code.
+
+in.mliap.unified.hippynn.Al.ghostneigh
+--------------------------------------
+
+Demonstrate MLIAP interface to HIPNN Al potential with ghost neighbors. Requires the HIPNN code.
+
+in.mliap.unified.hippynn.InP
+----------------------------
+
+Demonstrate MLIAP interface to HIPNN InP potential. Requires the HIPNN code.
+
+in.mliap.so3.Ni_Mo
+------------------
+
+Example of linear model with SO3 descriptors for Ni/Mo
+
+in.mliap.so3.nn.Si
+------------------
+
+Example of NN model with SO3 descriptors for Si
diff --git a/examples/mliap/in.mliap.unified.hippynn.Al b/examples/mliap/in.mliap.unified.hippynn.Al
index 2be0af158e..4efcba8b37 100644
--- a/examples/mliap/in.mliap.unified.hippynn.Al
+++ b/examples/mliap/in.mliap.unified.hippynn.Al
@@ -1,4 +1,5 @@
 # Demonstrate MLIAP interface to HIPNN Al potential
+# Example taken from https://github.com/lanl/hippynn
 
 # Initialize simulation
 
diff --git a/examples/mliap/in.mliap.unified.hippynn.Al.ghostneigh b/examples/mliap/in.mliap.unified.hippynn.Al.ghostneigh
index 923568f25f..7dfc982cbc 100644
--- a/examples/mliap/in.mliap.unified.hippynn.Al.ghostneigh
+++ b/examples/mliap/in.mliap.unified.hippynn.Al.ghostneigh
@@ -1,4 +1,5 @@
 # Demonstrate MLIAP interface to HIPNN Al potential
+# Example taken from https://github.com/lanl/hippynn
 
 # Initialize simulation
 
diff --git a/examples/mliap/in.mliap.unified.hippynn.InP b/examples/mliap/in.mliap.unified.hippynn.InP
index 2300f1b15b..0cba7f2416 100644
--- a/examples/mliap/in.mliap.unified.hippynn.InP
+++ b/examples/mliap/in.mliap.unified.hippynn.InP
@@ -1,4 +1,5 @@
 # Demonstrate MLIAP interface to HIPNN InP potential
+# Example taken from https://github.com/lanl/hippynn
 
 # Initialize simulation
 
diff --git a/src/ML-IAP/mliap_unified.cpp b/src/ML-IAP/mliap_unified.cpp
index c31cdd2c81..98befeb7c4 100644
--- a/src/ML-IAP/mliap_unified.cpp
+++ b/src/ML-IAP/mliap_unified.cpp
@@ -14,6 +14,7 @@
 /* ----------------------------------------------------------------------
    Contributing author: Steven Ray Anaya (LANL)
 ------------------------------------------------------------------------- */
+
 #ifdef MLIAP_PYTHON
 
 #include <Python.h>
@@ -38,9 +39,14 @@ MLIAPDummyDescriptor::MLIAPDummyDescriptor(LAMMPS *lmp) :
 MLIAPDummyDescriptor::~MLIAPDummyDescriptor() {
   memory->destroy(radelem);
   memory->destroy(cutsq);
+  // manually decrement borrowed reference from Python
   Py_DECREF(unified_interface);
 }
 
+/* ----------------------------------------------------------------------
+   invoke compute_descriptors from Cython interface
+   ---------------------------------------------------------------------- */
+
 void MLIAPDummyDescriptor::compute_descriptors(class MLIAPData *data) {
   PyGILState_STATE gstate = PyGILState_Ensure();
   compute_descriptors_python(unified_interface, data);
@@ -53,6 +59,10 @@ void MLIAPDummyDescriptor::compute_descriptors(class MLIAPData *data) {
   PyGILState_Release(gstate);
 }
 
+/* ----------------------------------------------------------------------
+   invoke compute_forces from Cython interface
+   ---------------------------------------------------------------------- */
+
 void MLIAPDummyDescriptor::compute_forces(class MLIAPData *data) {
   PyGILState_STATE gstate = PyGILState_Ensure();
   compute_forces_python(unified_interface, data);
@@ -65,10 +75,12 @@ void MLIAPDummyDescriptor::compute_forces(class MLIAPData *data) {
   PyGILState_Release(gstate);
 }
 
+// not implemented
 void MLIAPDummyDescriptor::compute_force_gradients(class MLIAPData *) {
   error->all(FLERR, "compute_force_gradients not implemented");
 }
 
+// not implemented
 void MLIAPDummyDescriptor::compute_descriptor_gradients(class MLIAPData *) {
   error->all(FLERR, "compute_descriptor_gradients not implemented");
 }
@@ -113,7 +125,7 @@ MLIAPDummyModel::MLIAPDummyModel(LAMMPS *lmp, char *coefffilename) :
 }
 
 MLIAPDummyModel::~MLIAPDummyModel() {
-  // MLIAPPY_unload_model(this);
+  // manually decrement borrowed reference from Python
   Py_DECREF(unified_interface);
 }
 
@@ -126,6 +138,10 @@ int MLIAPDummyModel::get_gamma_nnz(class MLIAPData *) {
   return 0;
 }
 
+/* ----------------------------------------------------------------------
+   invoke compute_gradients from Cython interface
+   ---------------------------------------------------------------------- */
+
 void MLIAPDummyModel::compute_gradients(class MLIAPData *data) {
   PyGILState_STATE gstate = PyGILState_Ensure();
   compute_gradients_python(unified_interface, data);
@@ -138,24 +154,33 @@ void MLIAPDummyModel::compute_gradients(class MLIAPData *data) {
   PyGILState_Release(gstate);
 }
 
+// not implemented
 void MLIAPDummyModel::compute_gradgrads(class MLIAPData *) {
   error->all(FLERR, "compute_gradgrads not implemented");
 }
 
+// not implemented
 void MLIAPDummyModel::compute_force_gradients(class MLIAPData *) {
   error->all(FLERR, "compute_force_gradients not implemented");
 }
 
+/* ----------------------------------------------------------------------
+   memory usage unclear due to Cython/Python implementation
+   ---------------------------------------------------------------------- */
+
 double MLIAPDummyModel::memory_usage() {
   // TODO: implement memory usage in Cython(?)
   return 0;
 }
 
+// not implemented
 void MLIAPDummyModel::read_coeffs(char *) {
   error->all(FLERR, "read_coeffs not implemented");
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   build the unified interface object, connect to dummy model and descriptor
+   ---------------------------------------------------------------------- */
 
 MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *data, LAMMPS *lmp,
     char *coefffilename) {
@@ -190,6 +215,7 @@ MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *dat
     lmp->error->all(FLERR, "Running mliappy unified module failure.");
   }
 
+  // Borrowed references must be manually incremented
   model->unified_interface = unified_interface;
   Py_INCREF(unified_interface);
   descriptor->unified_interface = unified_interface;
@@ -203,6 +229,10 @@ MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *dat
 
 /* ---------------------------------------------------------------------- */
 
+/* ----------------------------------------------------------------------
+   set energy for ij atom pairs
+   ---------------------------------------------------------------------- */
+
 void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij) {
     double e_total = 0;
     for (int ii = 0; ii < data->nlistatoms; ii++)
@@ -218,6 +248,10 @@ void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij) {
     data->energy = e_total;
 }
 
+/* ----------------------------------------------------------------------
+   set forces for ij atom pairs
+   ---------------------------------------------------------------------- */
+
 void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij) {
     double **f = data->f;
     for (int ii = 0; ii < data->npairs; ii++) {
diff --git a/src/ML-IAP/mliap_unified.h b/src/ML-IAP/mliap_unified.h
index cf39165353..dae4824d06 100644
--- a/src/ML-IAP/mliap_unified.h
+++ b/src/ML-IAP/mliap_unified.h
@@ -33,7 +33,7 @@ class MLIAPDummyDescriptor : public MLIAPDescriptor {
   virtual void init();
   void set_elements(char **, int);
 
-  PyObject *unified_interface;
+  PyObject *unified_interface;  // MLIAPUnifiedInterface
   double rcutfac;
 };
 
diff --git a/src/ML-IAP/mliap_unified_couple.pyx b/src/ML-IAP/mliap_unified_couple.pyx
index 6ec2de42f9..3fde99a25e 100644
--- a/src/ML-IAP/mliap_unified_couple.pyx
+++ b/src/ML-IAP/mliap_unified_couple.pyx
@@ -106,6 +106,8 @@ def write_only_property(fset):
     return property(fget=None, fset=fset)
 
 
+# Cython implementation of MLIAPData
+# Automatically converts between C arrays and numpy when needed
 cdef class MLIAPDataPy:
     cdef MLIAPData * data
 
@@ -296,6 +298,7 @@ cdef class MLIAPDataPy:
         return self.data.vflag
 
 
+# Interface between C and Python compute functions
 cdef class MLIAPUnifiedInterface:
     cdef MLIAPDummyModel * model
     cdef MLIAPDummyDescriptor * descriptor
@@ -334,6 +337,7 @@ cdef public void compute_forces_python(unified_int, MLIAPData *data) except * wi
     unified_int.compute_forces(pydata)
 
 
+# Create a MLIAPUnifiedInterface and connect it to the dummy model, descriptor
 cdef public object mliap_unified_connect(char *fname, MLIAPDummyModel * model,
                                          MLIAPDummyDescriptor * descriptor) with gil:
     str_fname = fname.decode('utf-8')
@@ -383,6 +387,7 @@ cdef public object mliap_unified_connect(char *fname, MLIAPDummyModel * model,
     return unified_int
 
 
+# For pre-loading a Python model
 def load_from_python(unified):
     global LOADED_MODEL
     LOADED_MODEL = unified