diff --git a/.gitignore b/.gitignore
index 29648886a7..ae708ff184 100644
--- a/.gitignore
+++ b/.gitignore
@@ -44,6 +44,7 @@ Thumbs.db
 /build*
 /CMakeCache.txt
 /CMakeFiles/
+/Testing
 /Makefile
 /Testing
 /cmake_install.cmake
diff --git a/doc/src/Library_utility.rst b/doc/src/Library_utility.rst
index b2f3666f88..32fac6bcc8 100644
--- a/doc/src/Library_utility.rst
+++ b/doc/src/Library_utility.rst
@@ -8,7 +8,11 @@ functions.  They do not directly call the LAMMPS library.
 - :cpp:func:`lammps_decode_image_flags`
 - :cpp:func:`lammps_set_fix_external_callback`
 - :cpp:func:`lammps_fix_external_set_energy_global`
+- :cpp:func:`lammps_fix_external_set_energy_peratom`
 - :cpp:func:`lammps_fix_external_set_virial_global`
+- :cpp:func:`lammps_fix_external_set_virial_peratom`
+- :cpp:func:`lammps_fix_external_set_vector_length`
+- :cpp:func:`lammps_fix_external_set_vector`
 - :cpp:func:`lammps_free`
 - :cpp:func:`lammps_is_running`
 - :cpp:func:`lammps_force_timeout`
@@ -33,7 +37,7 @@ where such memory buffers were allocated that require the use of
 
 -----------------------
 
-.. doxygenfunction:: lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*)
+.. doxygenfunction:: lammps_set_fix_external_callback(void *, const char *, FixExternalFnPtr, void*)
    :project: progguide
 
 -----------------------
@@ -43,11 +47,31 @@ where such memory buffers were allocated that require the use of
 
 -----------------------
 
+.. doxygenfunction:: lammps_fix_external_set_energy_peratom
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_fix_external_set_virial_global
    :project: progguide
 
 -----------------------
 
+.. doxygenfunction:: lammps_fix_external_set_virial_peratom
+   :project: progguide
+
+-----------------------
+
+.. doxygenfunction:: lammps_fix_external_set_vector_length
+   :project: progguide
+
+-----------------------
+
+.. doxygenfunction:: lammps_fix_external_set_vector
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_free
    :project: progguide
 
diff --git a/python/lammps/core.py b/python/lammps/core.py
index 2f101f4eab..d981243503 100644
--- a/python/lammps/core.py
+++ b/python/lammps/core.py
@@ -295,9 +295,16 @@ class lammps(object):
 
     self.lib.lammps_extract_variable.argtypes = [c_void_p, c_char_p, c_char_p]
 
-    # TODO: NOT IMPLEMENTED IN PYTHON WRAPPER
-    self.lib.lammps_fix_external_set_energy_global = [c_void_p, c_char_p, c_double]
-    self.lib.lammps_fix_external_set_virial_global = [c_void_p, c_char_p, POINTER(c_double)]
+    self.lib.lammps_fix_external_get_force.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_fix_external_get_force.restype = POINTER(POINTER(c_double))
+
+    self.lib.lammps_fix_external_set_energy_global.argtypes = [c_void_p, c_char_p, c_double]
+    self.lib.lammps_fix_external_set_virial_global.argtypes = [c_void_p, c_char_p, POINTER(c_double)]
+    self.lib.lammps_fix_external_set_energy_peratom.argtypes = [c_void_p, c_char_p, POINTER(c_double)]
+    self.lib.lammps_fix_external_set_virial_peratom.argtypes = [c_void_p, c_char_p, POINTER(POINTER(c_double))]
+
+    self.lib.lammps_fix_external_set_vector_length.argtypes = [c_void_p, c_char_p, c_int]
+    self.lib.lammps_fix_external_set_vector.argtypes = [c_void_p, c_char_p, c_int, c_double]
 
     # detect if Python is using a version of mpi4py that can pass communicators
     # only needed if LAMMPS has been compiled with MPI support.
@@ -1725,7 +1732,33 @@ class lammps(object):
 
   # -------------------------------------------------------------------------
 
-  def set_fix_external_callback(self, fix_name, callback, caller=None):
+  def set_fix_external_callback(self, fix_id, callback, caller=None):
+    """Set the callback function for a fix external instance with a given fix ID.
+
+    Optionally also set a reference to the calling object.
+
+    This is a wrapper around the :cpp:func:`lammps_set_fix_external_callback` function
+    of the C-library interface.  However this is set up to call a Python function with
+    the following arguments.
+
+    .. code-block: python
+
+       def func(object, ntimestep, nlocal, tag, x, f):
+
+    - object is the value of the "caller" argument
+    - ntimestep is the current timestep
+    - nlocal is the number of local atoms on the current MPI process
+    - tag is a 1d NumPy array of integers representing the atom IDs of the local atoms
+    - x is a 2d NumPy array of doubles of the coordinates of the local atoms
+    - f is a 2d NumPy array of doubles of the forces on the local atoms that will be added
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param callback: Python function that will be called from fix external
+    :type: function
+    :param caller: reference to some object passed to the callback function
+    :type: object, optional
+    """
     import numpy as np
 
     def callback_wrapper(caller, ntimestep, nlocal, tag_ptr, x_ptr, fext_ptr):
@@ -1737,10 +1770,148 @@ class lammps(object):
     cFunc   = self.FIX_EXTERNAL_CALLBACK_FUNC(callback_wrapper)
     cCaller = caller
 
-    self.callback[fix_name] = { 'function': cFunc, 'caller': caller }
+    self.callback[fix_id] = { 'function': cFunc, 'caller': caller }
     with ExceptionCheck(self):
-      self.lib.lammps_set_fix_external_callback(self.lmp, fix_name.encode(), cFunc, cCaller)
+      self.lib.lammps_set_fix_external_callback(self.lmp, fix_id.encode(), cFunc, cCaller)
 
+  # -------------------------------------------------------------------------
+
+  def fix_external_get_force(self, fix_id):
+    """Get access to the array with per-atom forces of a fix external instance with a given fix ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_get_force` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :return: requested data
+    :rtype: ctypes.POINTER(ctypes.POINTER(ctypes.double))
+    """
+
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_get_force(self.lmp, fix_id.encode())
+    return None
+
+  # -------------------------------------------------------------------------
+
+  def fix_external_set_energy_global(self, fix_id, eng):
+    """Set the global energy contribution for a fix external instance with the given ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_energy_global` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param eng:     potential energy value to be added by fix external
+    :type: float
+    """
+
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_energy_global(self.lmp, fix_id.encode(), eng)
+
+  # -------------------------------------------------------------------------
+
+  def fix_external_set_virial_global(self, fix_id, virial):
+    """Set the global virial contribution for a fix external instance with the given ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_virial_global` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param eng:     list of 6 floating point numbers with the virial to be added by fix external
+    :type: float
+    """
+
+    cvirial = (6*c_double)(*virial)
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_virial_global(self.lmp, fix_id.encode(), cvirial)
+
+  # -------------------------------------------------------------------------
+
+  def fix_external_set_energy_peratom(self, fix_id, eatom):
+    """Set the per-atom energy contribution for a fix external instance with the given ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_energy_peratom` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param eatom:   list of potential energy values for local atoms to be added by fix external
+    :type: float
+    """
+
+    nlocal = self.extract_setting('nlocal')
+    if len(eatom) < nlocal:
+      raise Exception('per-atom energy list length must be at least nlocal')
+    ceatom = (nlocal*c_double)(*eatom)
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_energy_peratom(self.lmp, fix_id.encode(), ceatom)
+
+  # -------------------------------------------------------------------------
+
+  def fix_external_set_virial_peratom(self, fix_id, vatom):
+    """Set the per-atom virial contribution for a fix external instance with the given ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_virial_peratom` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param vatom:     list of natoms lists with 6 floating point numbers to be added by fix external
+    :type: float
+    """
+
+    # copy virial data to C compatible buffer
+    nlocal = self.extract_setting('nlocal')
+    if len(vatom) < nlocal:
+      raise Exception('per-atom virial first dimension must be at least nlocal')
+    if len(vatom[0]) != 6:
+      raise Exception('per-atom virial second dimension must be 6')
+    vbuf = (c_double * 6)
+    vptr = POINTER(c_double)
+    c_virial = (vptr * nlocal)()
+    for i in range(nlocal):
+        c_virial[i] = vbuf()
+        for j in range(6):
+          c_virial[i][j] = vatom[i][j]
+
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_virial_peratom(self.lmp, fix_id.encode(), c_virial)
+
+  # -------------------------------------------------------------------------
+  def fix_external_set_vector_length(self, fix_id, length):
+    """Set the vector length for a global vector stored with fix external for analysis
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_vector_length` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param length:  length of the global vector
+    :type: int
+    """
+
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_vector_length(self.lmp, fix_id.encode(), length)
+
+  # -------------------------------------------------------------------------
+  def fix_external_set_vector(self, fix_id, idx, val):
+    """Store a global vector value for a fix external instance with the given ID.
+
+    This is a wrapper around the :cpp:func:`lammps_fix_external_set_vector` function
+    of the C-library interface.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param idx:     1-based index of the value in the global vector
+    :type: int
+    :param val:     value to be stored in the global vector
+    :type: float
+    """
+
+    with ExceptionCheck(self):
+      return self.lib.lammps_fix_external_set_vector(self.lmp, fix_id.encode(), idx, val)
 
   # -------------------------------------------------------------------------
 
diff --git a/python/lammps/numpy_wrapper.py b/python/lammps/numpy_wrapper.py
index 20fdf4cc68..6f4503a9c8 100644
--- a/python/lammps/numpy_wrapper.py
+++ b/python/lammps/numpy_wrapper.py
@@ -17,7 +17,7 @@
 ################################################################################
 
 import warnings
-from ctypes import POINTER, c_double, c_int, c_int32, c_int64, cast
+from ctypes import POINTER, c_void_p, c_char_p, c_double, c_int, c_int32, c_int64, cast
 
 
 from .constants import *                # lgtm [py/polluting-import]
@@ -246,7 +246,85 @@ class numpy_wrapper:
       return np.ctypeslib.as_array(value)
     return value
 
-  # -------------------------------------------------------------------------
+    # -------------------------------------------------------------------------
+
+  def fix_external_get_force(self, fix_id):
+    """Get access to the array with per-atom forces of a fix external instance with a given fix ID.
+
+    This function is a wrapper around the
+    :py:meth:`lammps.fix_external_get_force() <lammps.lammps.fix_external_get_force()>`
+    method.  It behaves the same as the original method, but returns a NumPy array instead
+    of a ``ctypes`` pointer.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :return: requested data
+    :rtype: numpy.array
+    """
+    import numpy as np
+    nlocal = self.lmp.extract_setting('nlocal')
+    value = self.lmp.fix_external_get_force(fix_id)
+    return self.darray(value,nlocal,3)
+
+    # -------------------------------------------------------------------------
+
+  def fix_external_set_energy_peratom(self, fix_id, eatom):
+    """Set the per-atom energy contribution for a fix external instance with the given ID.
+
+    This function is an alternative to
+    :py:meth:`lammps.fix_external_set_energy_peratom() <lammps.lammps.fix_external_set_energy_peratom()>`
+    method.  It behaves the same as the original method, but accepts a NumPy array
+    instead of a list as argument.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param eatom:   per-atom potential energy
+    :type: numpy.array
+    """
+    import numpy as np
+    nlocal = self.lmp.extract_setting('nlocal')
+    if len(eatom) < nlocal:
+      raise Exception('per-atom energy dimension must be at least nlocal')
+
+    c_double_p = POINTER(c_double)
+    value = eatom.astype(np.double)
+    return self.lmp.lib.lammps_fix_external_set_energy_peratom(self.lmp.lmp, fix_id.encode(),
+                                                               value.ctypes.data_as(c_double_p))
+
+    # -------------------------------------------------------------------------
+
+  def fix_external_set_virial_peratom(self, fix_id, vatom):
+    """Set the per-atom virial contribution for a fix external instance with the given ID.
+
+    This function is an alternative to
+    :py:meth:`lammps.fix_external_set_virial_peratom() <lammps.lammps.fix_external_set_virial_peratom()>`
+    method.  It behaves the same as the original method, but accepts a NumPy array
+    instead of a list as argument.
+
+    :param fix_id:  Fix-ID of a fix external instance
+    :type: string
+    :param eatom:   per-atom potential energy
+    :type: numpy.array
+    """
+    import numpy as np
+    nlocal = self.lmp.extract_setting('nlocal')
+    if len(vatom) < nlocal:
+      raise Exception('per-atom virial first dimension must be at least nlocal')
+    if len(vatom[0]) != 6:
+      raise Exception('per-atom virial second dimension must be 6')
+
+    c_double_pp = np.ctypeslib.ndpointer(dtype=np.uintp, ndim=1, flags='C')
+
+    # recast numpy array to be compatible with library interface
+    value = (vatom.__array_interface__['data'][0]
+                   + np.arange(vatom.shape[0])*vatom.strides[0]).astype(np.uintp)
+
+    # change prototype to our custom type
+    self.lmp.lib.lammps_fix_external_set_virial_peratom.argtypes = [ c_void_p, c_char_p, c_double_pp ]
+
+    self.lmp.lib.lammps_fix_external_set_virial_peratom(self.lmp.lmp, fix_id.encode(), value)
+
+    # -------------------------------------------------------------------------
 
   def get_neighlist(self, idx):
     """Returns an instance of :class:`NumPyNeighList` which wraps access to the neighbor list with the given index
diff --git a/src/MISC/fix_pair_tracker.cpp b/src/MISC/fix_pair_tracker.cpp
index 7939485aff..755025baed 100644
--- a/src/MISC/fix_pair_tracker.cpp
+++ b/src/MISC/fix_pair_tracker.cpp
@@ -246,7 +246,7 @@ void FixPairTracker::reallocate(int n)
 
 double FixPairTracker::memory_usage()
 {
-  double bytes = nmax * nvalues * sizeof(double);
+  double bytes = nmax * (double) nvalues * sizeof(double);
   bytes += nmax * 2 * sizeof(int);
   return bytes;
 }
diff --git a/src/fix_external.cpp b/src/fix_external.cpp
index 8bede4ee19..f95c30c28d 100644
--- a/src/fix_external.cpp
+++ b/src/fix_external.cpp
@@ -15,6 +15,7 @@
 #include "fix_external.h"
 
 #include "atom.h"
+#include "comm.h"
 #include "error.h"
 #include "memory.h"
 #include "update.h"
@@ -80,7 +81,7 @@ FixExternal::~FixExternal()
   atom->delete_callback(id,Atom::GROW);
 
   memory->destroy(fexternal);
-  delete [] caller_vector;
+  delete[] caller_vector;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -163,6 +164,12 @@ void FixExternal::post_force(int vflag)
         f[i][1] += fexternal[i][1];
         f[i][2] += fexternal[i][2];
       }
+
+    // add contribution to global virial from previously stored value
+
+    if (vflag_global)
+      for (int i = 0; i < 6; ++i)
+        virial[i] = user_virial[i];
   }
 }
 
@@ -179,10 +186,13 @@ void FixExternal::min_post_force(int vflag)
 
 /* ----------------------------------------------------------------------
    caller invokes this method to set its contribution to global energy
+   this is the *total* energy across all MPI ranks of the external code
+   and must be set for all MPI ranks.
    unlike other energy/virial set methods:
      do not just return if eflag_global is not set
      b/c input script could access this quantity via compute_scalar()
      even if eflag is not set on a particular timestep
+   this function is compatible with CALLBACK and ARRAY mode
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_energy_global(double caller_energy)
@@ -191,24 +201,32 @@ void FixExternal::set_energy_global(double caller_energy)
 }
 
 /* ----------------------------------------------------------------------
-   caller invokes this method to set its contribution to global virial
+   caller invokes this method to set its contribution to the global virial
+   for all MPI ranks.  the virial value is the *total* contribution across
+   all MPI ranks of the external code and thus we need to divide by the
+   number of MPI ranks since the tallying code expects per MPI rank contributions.
+   this function is compatible with PF_CALLBACK and PF_ARRAY mode
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_virial_global(double *caller_virial)
 {
-  if (!vflag_global) return;
-
+  const double npscale = 1.0/(double)comm->nprocs;
   for (int i = 0; i < 6; i++)
-    virial[i] = caller_virial[i];
+    user_virial[i] = npscale * caller_virial[i];
 }
 
 /* ----------------------------------------------------------------------
-   caller invokes this method to set its contribution to peratom energy
+   caller invokes this method to set its contribution to peratom energy.
+   this is applied to the *local* atoms only.
+   this function is compatible with PF_CALLBACK mode only since it tallies
+   its energy contributions directly into the accumulator arrays.
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_energy_peratom(double *caller_energy)
 {
   if (!eflag_atom) return;
+  if ((mode == PF_ARRAY) && (comm->me == 0))
+    error->warning(FLERR,"Can only set energy/atom for fix external in pf/callback mode");
 
   int nlocal = atom->nlocal;
   for (int i = 0; i < nlocal; i++)
@@ -217,6 +235,9 @@ void FixExternal::set_energy_peratom(double *caller_energy)
 
 /* ----------------------------------------------------------------------
    caller invokes this method to set its contribution to peratom virial
+   this is applied to the *local* atoms only.
+   this function is compatible with PF_CALLBACK mode only since it tallies
+   its virial contributions directly into the accumulator arrays.
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_virial_peratom(double **caller_virial)
@@ -224,6 +245,8 @@ void FixExternal::set_virial_peratom(double **caller_virial)
   int i,j;
 
   if (!vflag_atom) return;
+  if ((mode == PF_ARRAY) && (comm->me == 0))
+    error->warning(FLERR,"Can only set virial/atom for fix external in pf/callback mode");
 
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++)
@@ -232,13 +255,13 @@ void FixExternal::set_virial_peratom(double **caller_virial)
 }
 
 /* ----------------------------------------------------------------------
-   caller invokes this method to set length of vector of values
-   assume all vector values are extensive, could make this an option
+   caller invokes this method to set length of global vector of values
+   assume all vector values are extensive.
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_vector_length(int n)
 {
-  delete [] caller_vector;
+  delete[] caller_vector;
 
   vector_flag = 1;
   size_vector = n;
@@ -248,14 +271,15 @@ void FixExternal::set_vector_length(int n)
 }
 
 /* ----------------------------------------------------------------------
-   caller invokes this method to set Index value in vector
-   index ranges from 1 to N inclusive
+   caller invokes this method to set value for item at "index" in vector
+   index is 1-based, thus index ranges from 1 to N inclusively.
+   Must be called from all MPI ranks.
 ------------------------------------------------------------------------- */
 
 void FixExternal::set_vector(int index, double value)
 {
-  if (index >= size_vector)
-    error->all(FLERR,"Invalid set_vector index in fix external");
+  if (index > size_vector)
+    error->all(FLERR,"Invalid set_vector index ({} of {}) in fix external",index,size_vector);
   caller_vector[index-1] = value;
 }
 
@@ -286,6 +310,7 @@ double FixExternal::compute_vector(int n)
 double FixExternal::memory_usage()
 {
   double bytes = 3*atom->nmax * sizeof(double);
+  bytes += 6*sizeof(double);
   return bytes;
 }
 
@@ -342,3 +367,16 @@ void FixExternal::set_callback(FnPtr caller_callback, void *caller_ptr)
   callback = caller_callback;
   ptr_caller = caller_ptr;
 }
+
+/* ----------------------------------------------------------------------
+   get access to internal data structures
+------------------------------------------------------------------------- */
+
+void *FixExternal::extract(const char *str, int &dim)
+{
+  if (strcmp(str, "fexternal") == 0) {
+    dim = 2;
+    return (void *) fexternal;
+  }
+  return nullptr;
+}
diff --git a/src/fix_external.h b/src/fix_external.h
index 0ace978f99..16db5d5015 100644
--- a/src/fix_external.h
+++ b/src/fix_external.h
@@ -57,11 +57,14 @@ class FixExternal : public Fix {
   typedef void (*FnPtr)(void *, bigint, int, tagint *, double **, double **);
   void set_callback(FnPtr, void *);
 
+  void *extract(const char *, int &);
+
  private:
   int mode, ncall, napply, eflag_caller;
   FnPtr callback;
   void *ptr_caller;
   double user_energy;
+  double user_virial[6];
   double *caller_vector;
 };
 
diff --git a/src/library.cpp b/src/library.cpp
index d4a82489f0..4847c705af 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -4804,15 +4804,49 @@ void lammps_decode_image_flags(imageint image, int *flags)
   flags[2] = (image >> IMG2BITS) - IMGMAX;
 }
 
-/* ----------------------------------------------------------------------
-   find fix external with given ID and set the callback function
-   and caller pointer
-------------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */
 
-void lammps_set_fix_external_callback(void *handle, char *id, FixExternalFnPtr callback_ptr, void * caller)
+/** Set up the callback function for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+Fix :doc:`external <fix_external>` allows programs that are running LAMMPS through
+its library interface to modify certain LAMMPS properties on specific
+timesteps, similar to the way other fixes do.
+
+This function sets the callback function for use with the "pf/callback"
+mode. The function has to have C language bindings with the prototype:
+
+.. code-block:: c
+
+   void func(void *ptr, bigint timestep, int nlocal, tagint *ids, double **x, double **fexternal);
+
+The argument *ptr* to this function will be stored in fix external and
+the passed as the first argument calling the callback function `func()`.
+This would usually be a pointer to the active LAMMPS instance, i.e. the same
+pointer as the *handle* argument.  This would be needed to call
+functions that set the global or per-atom energy or virial contributions
+from within the callback function.
+
+The callback mechanism is one of the two modes of how forces and can be
+applied to a simulation with the help of fix external. The alternative
+is the array mode where you call :cpp:func:`lammps_fix_external_get_force`.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  funcptr  pointer to callback function
+ * \param  ptr      pointer to object in calling code, passed to callback function as first argument */
+
+void lammps_set_fix_external_callback(void *handle, const char *id, FixExternalFnPtr funcptr, void *ptr)
 {
   LAMMPS *lmp = (LAMMPS *) handle;
-  FixExternal::FnPtr callback = (FixExternal::FnPtr) callback_ptr;
+  FixExternal::FnPtr callback = (FixExternal::FnPtr) funcptr;
 
   BEGIN_CAPTURE
   {
@@ -4823,18 +4857,103 @@ void lammps_set_fix_external_callback(void *handle, char *id, FixExternalFnPtr c
     Fix *fix = lmp->modify->fix[ifix];
 
     if (strcmp("external",fix->style) != 0)
-      lmp->error->all(FLERR,"Fix '{}' is not of style "
-                                        "external!", id);
+      lmp->error->all(FLERR,"Fix '{}' is not of style 'external'", id);
 
-    FixExternal * fext = (FixExternal*) fix;
-    fext->set_callback(callback, caller);
+    FixExternal *fext = (FixExternal *) fix;
+    fext->set_callback(callback, ptr);
   }
   END_CAPTURE
 }
 
-/* set global energy contribution from fix external */
-void lammps_fix_external_set_energy_global(void *handle, char *id,
-                                           double energy)
+/** Get pointer to the force array storage in a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+Fix :doc:`external <fix_external>` allows programs that are running LAMMPS through
+its library interface to add or modify certain LAMMPS properties on specific
+timesteps, similar to the way other fixes do.
+
+This function provides access to the per-atom force storage in a fix
+external instance with the given fix-ID to be added to the individual
+atoms when using the "pf/array" mode.  The *fexternal* array can be
+accessed like other "native" per-atom arrays accessible via the
+:cpp:func:`lammps_extract_atom` function.  Please note that the array
+stores holds the forces for *local* atoms for each MPI ranks, in the
+order determined by the neighbor list build.  Because the underlying
+data structures can change as well as the order of atom as they migrate
+between MPI processes because of the domain decomposition
+parallelization, this function should be always called immediately
+before the forces are going to be set to get an up-to-date pointer.
+ You can use e.g. :cpp:func:`lammps_get_natoms` to obtain the number
+of local atoms `nlocal` and then assume the dimensions of the returned
+force array as ``double force[nlocal][3]``.
+
+This is an alternative to the callback mechanism in fix external set up by
+:cpp:func:`lammps_set_fix_external_callback`. The main difference is
+that this mechanism can be used when forces are be pre-computed and the
+control alternates between LAMMPS and the external code, while the
+callback mechanism can call the external code to compute the force when
+the fix is triggered and needs them.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle     pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id         fix ID of fix external instance
+ * \return            a pointer to the per-atom force array allocated by the fix */
+
+double **lammps_fix_external_get_force(void *handle, const char *id)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+  double **fexternal = nullptr;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0)
+      lmp->error->all(FLERR,"Can not find fix with ID '{}'!", id);
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0)
+      lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
+
+    fexternal = (double **)fix->extract("fexternal",ifix);
+  }
+  END_CAPTURE
+  return fexternal;
+}
+
+/** Set the global energy contribution for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback` and
+:cpp:func:`lammps_fix_external_get_force` to also set the contribution
+to the global energy from the external code.  The value of the *eng*
+argument will be stored in the fix and applied on the current and all
+following timesteps until changed by another call to this function.
+The energy is in energy units as determined by the current :doc:`units <units>`
+settings and is the **total** energy of the contribution.  Thus when
+running in parallel all MPI processes have to call this function with
+the **same** value and this will be returned as scalar property of the
+fix external instance when accessed in LAMMPS input commands or from
+variables.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  eng      total energy to be added to the global energy */
+
+void lammps_fix_external_set_energy_global(void *handle, const char *id, double eng)
 {
   LAMMPS *lmp = (LAMMPS *) handle;
 
@@ -4849,15 +4968,41 @@ void lammps_fix_external_set_energy_global(void *handle, char *id,
     if (strcmp("external",fix->style) != 0)
       lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
 
-    FixExternal * fext = (FixExternal*) fix;
-    fext->set_energy_global(energy);
+    FixExternal *fext = (FixExternal*) fix;
+    fext->set_energy_global(eng);
   }
   END_CAPTURE
 }
 
-/* set global virial contribution from fix external */
-void lammps_fix_external_set_virial_global(void *handle, char *id,
-                                           double *virial)
+/** Set the global virial contribution for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback`
+and :cpp:func:`lammps_fix_external_get_force` to set the contribution to
+the global virial from the external code.
+
+The 6 values of the *virial* array will be stored in the fix and applied
+on the current and all following timesteps until changed by another call
+to this function. The components of the virial need to be stored in the
+order: *xx*, *yy*, *zz*, *xy*, *xz*, *yz*.  In LAMMPS the virial is
+stored internally as `stress*volume` in units of `pressure*volume` as
+determined by the current :doc:`units <units>` settings and is the
+**total** contribution.  Thus when running in parallel all MPI processes
+have to call this function with the **same** value and this will then
+be added by fix external.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  virial   the 6 global stress tensor components to be added to the global virial */
+
+void lammps_fix_external_set_virial_global(void *handle, const char *id, double *virial)
 {
   LAMMPS *lmp = (LAMMPS *) handle;
 
@@ -4878,6 +5023,207 @@ void lammps_fix_external_set_virial_global(void *handle, char *id,
   END_CAPTURE
 }
 
+/** Set the per-atom energy contribution for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback`
+to set the per-atom energy contribution due to the fix from the external code
+as part of the callback function.  For this to work, the handle to the
+LAMMPS object must be passed as the *ptr* argument when registering the
+callback function.
+
+.. note::
+
+   This function is fully independent from :cpp:func:`lammps_fix_external_set_energy_global`
+   and will **NOT** add any contributions to the global energy tally
+   and **NOT** check whether the sum of the contributions added here are
+   consistent with the global added energy.
+
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  eng      pointer to array of length nlocal with the energy to be added to the per-atom energy */
+
+void lammps_fix_external_set_energy_peratom(void *handle, const char *id, double *eng)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0)
+      lmp->error->all(FLERR,"Can not find fix with ID '{}'!", id);
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0)
+      lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
+
+    FixExternal *fext = (FixExternal*) fix;
+    fext->set_energy_peratom(eng);
+  }
+  END_CAPTURE
+}
+
+/** Set the per-atom virial contribution for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback`
+to set the per-atom virial contribution due to the fix from the external code
+as part of the callback function.  For this to work, the handle to the
+LAMMPS object must be passed as the *ptr* argument when registering the
+callback function.
+
+.. note::
+
+   This function is fully independent from :cpp:func:`lammps_fix_external_set_virial_global`
+   and will **NOT** add any contributions to the global virial tally
+   and **NOT** check whether the sum of the contributions added here are
+   consistent with the global added virial.
+
+The order and units of the per-atom stress tensor elements are the same
+as for the global virial.  The code in fix external assumes the
+dimensions of the per-atom virial array is ``double virial[nlocal][6]``.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  virial   a list of nlocal entries with the 6 per-atom stress tensor components to be added to the per-atom virial */
+
+void lammps_fix_external_set_virial_peratom(void *handle, const char *id, double **virial)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0)
+      lmp->error->all(FLERR,"Can not find fix with ID '{}'!", id);
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0)
+      lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
+
+    FixExternal * fext = (FixExternal*) fix;
+    fext->set_virial_peratom(virial);
+  }
+  END_CAPTURE
+}
+
+/** Set the vector length for a global vector stored with fix external for analysis
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback` and
+:cpp:func:`lammps_fix_external_get_force` to set the length of a global vector of
+properties that will be stored with the fix via
+:cpp:func:`lammps_fix_external_set_vector`.
+
+This function needs to be called **before** a call to
+:cpp:func:`lammps_fix_external_set_vector` and **before** a run or minimize
+command. When running in parallel it must be called from **all** MPI
+processes and with the same length parameter.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  len      length of the global vector to be stored with the fix */
+
+void lammps_fix_external_set_vector_length(void *handle, const char *id, int len)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0)
+      lmp->error->all(FLERR,"Can not find fix with ID '{}'!", id);
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0)
+      lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
+
+    FixExternal *fext = (FixExternal*) fix;
+    fext->set_vector_length(len);
+  }
+  END_CAPTURE
+}
+
+/** Store a global vector value for a fix external instance with the given ID.
+
+\verbatim embed:rst
+
+This is a companion function to :cpp:func:`lammps_set_fix_external_callback` and
+:cpp:func:`lammps_fix_external_get_force` to set the values of a global vector of
+properties that will be stored with the fix.  And can be accessed from
+within LAMMPS input commands (e.g. fix ave/time or variables) when used
+in a vector context.
+
+This function needs to be called **after** a call to
+:cpp:func:`lammps_fix_external_set_vector_length` and the  and **before** a run or minimize
+command.  When running in parallel it must be called from **all** MPI
+processes and with the **same** index and value parameters.  The value
+is assumed to be extensive.
+
+.. note::
+
+   The index in the *idx* parameter is 1-based, i.e. the first element
+   is set with idx = 1 and the last element of the vector with idx = N,
+   where N is the value of the *len* parameter of the call to
+   :cpp:func:`lammps_fix_external_set_vector_length`.
+
+Please see the documentation for :doc:`fix external <fix_external>` for
+more information about how to use the fix and how to couple it with an
+external code.
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param  id       fix ID of fix external instance
+ * \param  idx      1-based index of in global vector
+ * \param  val      value to be stored in global vector */
+
+void lammps_fix_external_set_vector(void *handle, const char *id, int idx, double val)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0)
+      lmp->error->all(FLERR,"Can not find fix with ID '{}'!", id);
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0)
+      lmp->error->all(FLERR,"Fix '{}' is not of style external!", id);
+
+    FixExternal * fext = (FixExternal*) fix;
+    fext->set_vector(idx, val);
+  }
+  END_CAPTURE
+}
+
 /* ---------------------------------------------------------------------- */
 
 /** Free memory buffer allocated by LAMMPS.
diff --git a/src/library.h b/src/library.h
index a337d7b510..654eda38fa 100644
--- a/src/library.h
+++ b/src/library.h
@@ -226,16 +226,21 @@ void lammps_decode_image_flags(int64_t image, int *flags);
 
 #if defined(LAMMPS_BIGBIG)
 typedef void (*FixExternalFnPtr)(void *, int64_t, int, int64_t *, double **, double **);
-void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void *);
+void lammps_set_fix_external_callback(void *handle, const char *id, FixExternalFnPtr funcptr, void *ptr);
 #elif defined(LAMMPS_SMALLBIG)
 typedef void (*FixExternalFnPtr)(void *, int64_t, int, int *, double **, double **);
-void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void *);
+void lammps_set_fix_external_callback(void *, const char *, FixExternalFnPtr, void *);
 #else
 typedef void (*FixExternalFnPtr)(void *, int, int, int *, double **, double **);
-void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void *);
+void lammps_set_fix_external_callback(void *, const char *, FixExternalFnPtr, void *);
 #endif
-void lammps_fix_external_set_energy_global(void *, char *, double);
-void lammps_fix_external_set_virial_global(void *, char *, double *);
+double **lammps_fix_external_get_force(void *handle, const char *id);
+void lammps_fix_external_set_energy_global(void *handle, const char *id, double eng);
+void lammps_fix_external_set_energy_peratom(void *handle, const char *id, double *eng);
+void lammps_fix_external_set_virial_global(void *handle, const char *id, double *virial);
+void lammps_fix_external_set_virial_peratom(void *handle, const char *id, double **virial);
+void lammps_fix_external_set_vector_length(void *handle, const char *id, int len);
+void lammps_fix_external_set_vector(void *handle, const char *id, int idx, double val);
 
 void lammps_free(void *ptr);
 
diff --git a/tools/swig/lammps.i b/tools/swig/lammps.i
index 56547dda53..5bf47f2463 100644
--- a/tools/swig/lammps.i
+++ b/tools/swig/lammps.i
@@ -128,16 +128,27 @@ extern int    lammps_id_name(void *, const char *, int, char *buffer, int buf_si
 extern int    lammps_plugin_count();
 extern int    lammps_plugin_name(int, char *, char *, int);
 /*
-extern int lammps_encode_image_flags(int ix, int iy, int iz);
-extern void lammps_decode_image_flags(int image, int *flags);
-extern int64_t lammps_encode_image_flags(int ix, int iy, int iz);
-extern void lammps_decode_image_flags(int64_t image, int *flags);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_fix_external_set_energy_global(void *, char *, double);
-extern void lammps_fix_external_set_virial_global(void *, char *, double *);
+ * Have not found a good way to map these functions in a general way.
+ * So some individual customization for the specific use case and compilation is needed.
+ *
+  extern int lammps_encode_image_flags(int ix, int iy, int iz);
+  extern void lammps_decode_image_flags(int image, int *flags);
+  extern int64_t lammps_encode_image_flags(int ix, int iy, int iz);
+  extern void lammps_decode_image_flags(int64_t image, int *flags);
+
+ * Supporting the fix external callback mechanism will require extra code specific to the application.
+  typedef void (*FixExternalFnPtr)(void *, int64_t, int, int64_t *, double **, double **);
+  extern  void lammps_set_fix_external_callback(void *handle, const char *id, FixExternalFnPtr funcptr, void *ptr);
+ * these two functions can only be used from the callback, so we don't support them either
+  extern  void lammps_fix_external_set_energy_peratom(void *handle, const char *id, double *eng);
+  extern void lammps_fix_external_set_virial_peratom(void *handle, const char *id, double **virial);
 */
+extern double **lammps_fix_external_get_force(void *handle, const char *id);
+extern void   lammps_fix_external_set_energy_global(void *handle, const char *id, double eng);
+extern void   lammps_fix_external_set_virial_global(void *handle, const char *id, double *virial);
+extern void   lammps_fix_external_set_vector_length(void *handle, const char *id, int len);
+extern void   lammps_fix_external_set_vector(void *handle, const char *id, int idx, double val);
+
 extern void   lammps_free(void *ptr);
 extern int    lammps_is_running(void *handle);
 extern void   lammps_force_timeout(void *handle);
@@ -255,16 +266,21 @@ extern int lammps_encode_image_flags(int ix, int iy, int iz);
 extern void lammps_decode_image_flags(int image, int *flags);
 extern int64_t lammps_encode_image_flags(int ix, int iy, int iz);
 extern void lammps_decode_image_flags(int64_t image, int *flags);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
-extern void lammps_fix_external_set_energy_global(void *, char *, double);
-extern void lammps_fix_external_set_virial_global(void *, char *, double *);
+typedef void (*FixExternalFnPtr)(void *, int64_t, int, int64_t *, double **, double **);
+extern  void lammps_set_fix_external_callback(void *handle, const char *id, FixExternalFnPtr funcptr, void *ptr);
+extern  void lammps_fix_external_set_energy_peratom(void *handle, const char *id, double *eng);
+extern void lammps_fix_external_set_virial_peratom(void *handle, const char *id, double **virial);
 */
+extern double **lammps_fix_external_get_force(void *handle, const char *id);
+extern void   lammps_fix_external_set_energy_global(void *handle, const char *id, double eng);
+extern void   lammps_fix_external_set_virial_global(void *handle, const char *id, double *virial);
+extern void   lammps_fix_external_set_vector_length(void *handle, const char *id, int len);
+extern void   lammps_fix_external_set_vector(void *handle, const char *id, int idx, double val);
+
 extern void   lammps_free(void *ptr);
 extern int    lammps_is_running(void *handle);
 extern void   lammps_force_timeout(void *handle);
 extern int    lammps_has_error(void *handle);
 extern int    lammps_get_last_error_message(void *handle, char *buffer, int buf_size);
 
-/* last revised for LAMMPS 8 April 2021 */
+/* last revised on 21 July 2021 */
diff --git a/unittest/c-library/CMakeLists.txt b/unittest/c-library/CMakeLists.txt
index 3c24cdcff4..b01cd64677 100644
--- a/unittest/c-library/CMakeLists.txt
+++ b/unittest/c-library/CMakeLists.txt
@@ -7,6 +7,10 @@ add_executable(test_library_commands test_library_commands.cpp test_main.cpp)
 target_link_libraries(test_library_commands PRIVATE lammps GTest::GTest GTest::GMock)
 add_test(LibraryCommands test_library_commands)
 
+add_executable(test_library_external test_library_external.cpp test_main.cpp)
+target_link_libraries(test_library_external PRIVATE lammps GTest::GTest GTest::GMock)
+add_test(LibraryExternal test_library_external)
+
 add_executable(test_library_properties test_library_properties.cpp test_main.cpp)
 target_link_libraries(test_library_properties PRIVATE lammps GTest::GTest GTest::GMock)
 target_compile_definitions(test_library_properties PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
diff --git a/unittest/c-library/test_library_external.cpp b/unittest/c-library/test_library_external.cpp
new file mode 100644
index 0000000000..005b31fcab
--- /dev/null
+++ b/unittest/c-library/test_library_external.cpp
@@ -0,0 +1,200 @@
+// unit tests creating LAMMPS instances via the library interface
+
+#include "library.h"
+
+#include <cinttypes>
+#include <string>
+
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+#include "test_main.h"
+
+using ::testing::HasSubstr;
+using ::testing::StartsWith;
+
+extern "C" {
+#ifdef LAMMPS_SMALLSMALL
+typedef int32_t step_t;
+typedef int32_t tag_t;
+#elif LAMMPS_SMALLBIG
+typedef int64_t step_t;
+typedef int32_t tag_t;
+#else
+typedef int64_t step_t;
+typedef int64_t tag_t;
+#endif
+static void callback(void *handle, step_t timestep, int nlocal, tag_t *, double **, double **f)
+{
+    for (int i = 0; i < nlocal; ++i)
+        f[i][0] = f[i][1] = f[i][2] = (double)timestep;
+
+    double v[6] = {1.0, 1.0, 1.0, 0.0, 0.0, 0.0};
+    lammps_fix_external_set_virial_global(handle, "ext", v);
+    if (timestep < 10) {
+        lammps_fix_external_set_energy_global(handle, "ext", 0.5);
+        lammps_fix_external_set_vector(handle, "ext", 1, timestep);
+        lammps_fix_external_set_vector(handle, "ext", 3, 1.0);
+        lammps_fix_external_set_vector(handle, "ext", 4, -0.25);
+    } else {
+        lammps_fix_external_set_energy_global(handle, "ext", 1.0);
+        lammps_fix_external_set_vector(handle, "ext", 2, timestep);
+        lammps_fix_external_set_vector(handle, "ext", 5, -1.0);
+        lammps_fix_external_set_vector(handle, "ext", 6, 0.25);
+    }
+    double *eatom  = new double[nlocal];
+    double **vatom = new double *[nlocal];
+    vatom[0]       = new double[nlocal * 6];
+    eatom[0]       = 0.0;
+    vatom[0][0] = vatom[0][1] = vatom[0][2] = vatom[0][3] = vatom[0][4] = vatom[0][5] = 0.0;
+
+    for (int i = 1; i < nlocal; ++i) {
+        eatom[i]    = 0.1 * i;
+        vatom[i]    = vatom[0] + 6 * i;
+        vatom[i][0] = vatom[i][1] = vatom[i][2] = 0.1;
+        vatom[i][3] = vatom[i][4] = vatom[i][5] = -0.2;
+    }
+    lammps_fix_external_set_energy_peratom(handle, "ext", eatom);
+    lammps_fix_external_set_virial_peratom(handle, "ext", vatom);
+    delete[] eatom;
+    delete[] vatom[0];
+    delete[] vatom;
+}
+}
+
+TEST(lammps_external, callback)
+{
+    const char *args[] = {"liblammps", "-log", "none", "-nocite"};
+    char **argv        = (char **)args;
+    int argc           = sizeof(args) / sizeof(char *);
+
+    ::testing::internal::CaptureStdout();
+    void *handle       = lammps_open_no_mpi(argc, argv, NULL);
+    std::string output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    ::testing::internal::CaptureStdout();
+    lammps_commands_string(handle, "lattice sc 1.0\n"
+                                   "region box block -1 1 -1 1 -1 1\n"
+                                   "create_box 1 box\n"
+                                   "create_atoms 1 box\n"
+                                   "mass 1 1.0\n"
+                                   "pair_style zero 0.1\n"
+                                   "pair_coeff 1 1\n"
+                                   "velocity all set 0.1 0.0 -0.1\n"
+                                   "fix 1 all nve\n"
+                                   "fix ext all external pf/callback 5 1\n"
+                                   "compute eatm all pe/atom fix\n"
+                                   "compute vatm all stress/atom NULL fix\n"
+                                   "compute sum all reduce sum c_eatm c_vatm[*]\n"
+                                   "thermo_style custom step temp pe ke etotal press c_sum[*]\n"
+                                   "thermo 5\n"
+                                   "fix_modify ext energy yes virial yes\n");
+    lammps_fix_external_set_vector_length(handle, "ext", 6);
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    ::testing::internal::CaptureStdout();
+    lammps_set_fix_external_callback(handle, "ext", &callback, handle);
+    lammps_command(handle, "run 10 post no");
+    double temp  = lammps_get_thermo(handle, "temp");
+    double pe    = lammps_get_thermo(handle, "pe");
+    double press = lammps_get_thermo(handle, "press");
+    double val   = 0.0;
+    double *valp;
+    for (int i = 0; i < 6; ++i) {
+        valp = (double *)lammps_extract_fix(handle, "ext", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR, i, 0);
+        val += *valp;
+        lammps_free(valp);
+    }
+    double *reduce =
+        (double *)lammps_extract_compute(handle, "sum", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR);
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    EXPECT_DOUBLE_EQ(temp, 1.0 / 30.0);
+    EXPECT_DOUBLE_EQ(pe, 1.0 / 8.0);
+    EXPECT_DOUBLE_EQ(press, 0.15416666666666667);
+    EXPECT_DOUBLE_EQ(val, 15);
+    EXPECT_DOUBLE_EQ(reduce[0], 2.8);
+    EXPECT_DOUBLE_EQ(reduce[1], -0.7);
+    EXPECT_DOUBLE_EQ(reduce[2], -0.7);
+    EXPECT_DOUBLE_EQ(reduce[3], -0.7);
+    EXPECT_DOUBLE_EQ(reduce[4], 1.4);
+    EXPECT_DOUBLE_EQ(reduce[5], 1.4);
+    EXPECT_DOUBLE_EQ(reduce[6], 1.4);
+
+    ::testing::internal::CaptureStdout();
+    lammps_close(handle);
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+}
+
+TEST(lammps_external, array)
+{
+    const char *args[] = {"liblammps", "-log", "none", "-nocite"};
+    char **argv        = (char **)args;
+    int argc           = sizeof(args) / sizeof(char *);
+
+    ::testing::internal::CaptureStdout();
+    void *handle       = lammps_open_no_mpi(argc, argv, NULL);
+    std::string output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    ::testing::internal::CaptureStdout();
+    lammps_commands_string(handle, "lattice sc 1.0\n"
+                                   "region box block -1 1 -1 1 -1 1\n"
+                                   "create_box 1 box\n"
+                                   "create_atoms 1 box\n"
+                                   "mass 1 1.0\n"
+                                   "pair_style zero 0.1\n"
+                                   "pair_coeff 1 1\n"
+                                   "velocity all set 0.1 0.0 -0.1\n"
+                                   "fix 1 all nve\n"
+                                   "fix ext all external pf/array 1\n"
+                                   "thermo 5\n");
+
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    ::testing::internal::CaptureStdout();
+    double **force = lammps_fix_external_get_force(handle, "ext");
+    int nlocal     = lammps_extract_setting(handle, "nlocal");
+    for (int i = 0; i < nlocal; ++i)
+        force[i][0] = force[i][1] = force[i][2] = 0.0;
+    lammps_fix_external_set_energy_global(handle, "ext", 0.5);
+    double v[6] = {0.5, 0.5, 0.5, 0.0, 0.0, 0.0};
+    lammps_fix_external_set_virial_global(handle, "ext", v);
+    lammps_command(handle, "run 5 post no");
+    double temp  = lammps_get_thermo(handle, "temp");
+    double pe    = lammps_get_thermo(handle, "pe");
+    double press = lammps_get_thermo(handle, "press");
+    output       = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    EXPECT_DOUBLE_EQ(temp, 4.0 / 525.0);
+    EXPECT_DOUBLE_EQ(pe, 1.0 / 16.0);
+    EXPECT_DOUBLE_EQ(press, 0.069166666666666668);
+
+    ::testing::internal::CaptureStdout();
+    nlocal = lammps_extract_setting(handle, "nlocal");
+    force  = lammps_fix_external_get_force(handle, "ext");
+    for (int i = 0; i < nlocal; ++i)
+        force[i][0] = force[i][1] = force[i][2] = 6.0;
+    lammps_fix_external_set_energy_global(handle, "ext", 1.0);
+    v[0] = v[1] = v[2] = 1.0;
+    v[3] = v[4] = v[5] = 0.0;
+    lammps_fix_external_set_virial_global(handle, "ext", v);
+    lammps_command(handle, "run 5 post no");
+    temp   = lammps_get_thermo(handle, "temp");
+    pe     = lammps_get_thermo(handle, "pe");
+    press  = lammps_get_thermo(handle, "press");
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    EXPECT_DOUBLE_EQ(temp, 1.0 / 30.0);
+    EXPECT_DOUBLE_EQ(pe, 1.0 / 8.0);
+    EXPECT_DOUBLE_EQ(press, 0.15416666666666667);
+
+    ::testing::internal::CaptureStdout();
+    lammps_close(handle);
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+}
diff --git a/unittest/python/CMakeLists.txt b/unittest/python/CMakeLists.txt
index b51d6e340a..9c9b7832ad 100644
--- a/unittest/python/CMakeLists.txt
+++ b/unittest/python/CMakeLists.txt
@@ -85,6 +85,11 @@ if(Python_EXECUTABLE)
            COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-formats.py -v
            WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
   set_tests_properties(PythonFormats PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
+
+  add_test(NAME PythonFixExternal
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-fix-external.py -v
+           WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+  set_tests_properties(PythonFixExternal PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
 else()
   message(STATUS "Skipping Tests for the LAMMPS Python Module: no suitable Python interpreter")
 endif()
diff --git a/unittest/python/python-fix-external.py b/unittest/python/python-fix-external.py
new file mode 100644
index 0000000000..c061b9f466
--- /dev/null
+++ b/unittest/python/python-fix-external.py
@@ -0,0 +1,176 @@
+import sys,os,unittest
+from ctypes import *
+from lammps import lammps, LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR
+
+try:
+    import numpy
+    NUMPY_INSTALLED = True
+except ImportError:
+    NUMPY_INSTALLED = False
+
+# add timestep dependent force
+def callback_one(lmp, ntimestep, nlocal, tag, x, f):
+    lmp.fix_external_set_virial_global("ext",[1.0, 1.0, 1.0, 0.0, 0.0, 0.0])
+    for i in range(nlocal):
+        f[i][0] = float(ntimestep)
+        f[i][1] = float(ntimestep)
+        f[i][2] = float(ntimestep)
+    if ntimestep < 10:
+        lmp.fix_external_set_energy_global("ext", 0.5)
+        lmp.fix_external_set_vector("ext", 1, ntimestep)
+        lmp.fix_external_set_vector("ext", 3, 1.0)
+        lmp.fix_external_set_vector("ext", 4, -0.25)
+    else:
+        lmp.fix_external_set_energy_global("ext", 1.0)
+        lmp.fix_external_set_vector("ext", 2, ntimestep)
+        lmp.fix_external_set_vector("ext", 5, -1.0)
+        lmp.fix_external_set_vector("ext", 6, 0.25)
+
+    eatom = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7]
+    vatom = [ [0.1,0.0,0.0,0.0,0.0,0.0],
+              [0.0,0.2,0.0,0.0,0.0,0.0],
+              [0.0,0.0,0.3,0.0,0.0,0.0],
+              [0.0,0.0,0.0,0.4,0.0,0.0],
+              [0.0,0.0,0.0,0.0,0.5,0.0],
+              [0.0,0.0,0.0,0.0,0.0,0.6],
+              [0.0,0.0,0.0,0.0,-7.0,0.0],
+              [0.0,-8.0,0.0,0.0,0.0,0.0] ]
+    if ntimestep < 5:
+        lmp.fix_external_set_energy_peratom("ext",eatom)
+        lmp.fix_external_set_virial_peratom("ext",vatom)
+    else:
+        import numpy as np
+        eng = np.array(eatom)
+        vir = np.array(vatom)
+
+        lmp.numpy.fix_external_set_energy_peratom("ext",eng)
+        lmp.numpy.fix_external_set_virial_peratom("ext",vir)
+
+    # ------------------------------------------------------------------------
+
+class PythonExternal(unittest.TestCase):
+    @unittest.skipIf(not NUMPY_INSTALLED, "NumPy is not available")
+    def testExternalCallback(self):
+        """Test fix external from Python with pf/callback"""
+
+        machine=None
+        if 'LAMMPS_MACHINE_NAME' in os.environ:
+            machine=os.environ['LAMMPS_MACHINE_NAME']
+        lmp=lammps(name=machine, cmdargs=['-nocite', '-log','none', '-echo', 'screen'])
+
+        # a few commands to set up simple system
+        basic_system="""lattice sc 1.0
+                        region box block -1 1 -1 1 -1 1
+                        create_box 1 box
+                        create_atoms 1 box
+                        mass 1 1.0
+                        pair_style zero 0.1
+                        pair_coeff 1 1
+                        velocity all set 0.1 0.0 -0.1
+                        thermo_style custom step temp pe ke etotal press
+                        thermo 5
+                        fix 1 all nve
+                        fix ext all external pf/callback 5 1
+                        compute eatm all pe/atom fix
+                        compute vatm all stress/atom NULL fix
+                        compute sum all reduce sum c_eatm c_vatm[*]
+                        thermo_style custom step temp pe ke etotal press c_sum[*]
+                        fix_modify ext energy yes virial yes
+"""
+        lmp.commands_string(basic_system)
+        lmp.fix_external_set_vector_length("ext",6);
+        lmp.set_fix_external_callback("ext",callback_one,lmp)
+
+        # check setting per-atom data with python lists
+        lmp.command("run 0 post no")
+        reduce = lmp.extract_compute("sum", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR)
+        self.assertAlmostEqual(reduce[0],2.8,14)
+        self.assertAlmostEqual(reduce[1],-0.1,14)
+        self.assertAlmostEqual(reduce[2],7.8,14)
+        self.assertAlmostEqual(reduce[3],-0.3,14)
+        self.assertAlmostEqual(reduce[4],-0.4,14)
+        self.assertAlmostEqual(reduce[5],6.5,14)
+        self.assertAlmostEqual(reduce[6],-0.6,14)
+
+        lmp.command("run 10 post no")
+        self.assertAlmostEqual(lmp.get_thermo("temp"),1.0/30.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("pe"),1.0/8.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("press"),0.15416666666666667,14)
+        # check setting per-atom data numpy arrays
+        reduce = lmp.extract_compute("sum", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR)
+        self.assertAlmostEqual(reduce[0],2.8,14)
+        self.assertAlmostEqual(reduce[1],-0.1,14)
+        self.assertAlmostEqual(reduce[2],7.8,14)
+        self.assertAlmostEqual(reduce[3],-0.3,14)
+        self.assertAlmostEqual(reduce[4],-0.4,14)
+        self.assertAlmostEqual(reduce[5],6.5,14)
+        self.assertAlmostEqual(reduce[6],-0.6,14)
+        val = 0.0
+        for i in range(0,6):
+            val += lmp.extract_fix("ext",LMP_STYLE_GLOBAL,LMP_TYPE_VECTOR,nrow=i)
+        self.assertAlmostEqual(val,15.0,14)
+
+    def testExternalArray(self):
+        """Test fix external from Python with pf/array"""
+
+        machine=None
+        if 'LAMMPS_MACHINE_NAME' in os.environ:
+            machine=os.environ['LAMMPS_MACHINE_NAME']
+        lmp=lammps(name=machine, cmdargs=['-nocite', '-log','none', '-echo', 'screen'])
+
+        # a few commands to set up simple system
+        basic_system="""lattice sc 1.0
+                        region box block -1 1 -1 1 -1 1
+                        create_box 1 box
+                        create_atoms 1 box
+                        mass 1 1.0
+                        pair_style zero 0.1
+                        pair_coeff 1 1
+                        velocity all set 0.1 0.0 -0.1
+                        fix 1 all nve
+                        fix ext all external pf/array 1
+                        fix_modify ext energy yes virial yes
+                        thermo_style custom step temp pe ke press
+                        thermo 5
+"""
+        lmp.commands_string(basic_system)
+        force = lmp.fix_external_get_force("ext");
+        nlocal = lmp.extract_setting("nlocal");
+        for i in range(nlocal):
+            force[i][0] = 0.0
+            force[i][1] = 0.0
+            force[i][2] = 0.0
+        lmp.fix_external_set_energy_global("ext", 0.5)
+        lmp.fix_external_set_virial_global("ext",[0.5, 0.5, 0.5, 0.0, 0.0, 0.0])
+
+        lmp.command("run 5 post no")
+        self.assertAlmostEqual(lmp.get_thermo("temp"),4.0/525.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("pe"),1.0/16.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("press"),0.06916666666666667,14)
+        if NUMPY_INSTALLED:
+            npforce = lmp.numpy.fix_external_get_force("ext")
+            self.assertEqual(len(npforce),8)
+            self.assertEqual(len(npforce[0]),3)
+            self.assertEqual(npforce[1][1],0.0)
+
+        force = lmp.fix_external_get_force("ext");
+        nlocal = lmp.extract_setting("nlocal");
+        for i in range(nlocal):
+            force[i][0] = 6.0
+            force[i][1] = 6.0
+            force[i][2] = 6.0
+        lmp.fix_external_set_energy_global("ext", 1.0)
+        lmp.fix_external_set_virial_global("ext",[1.0, 1.0, 1.0, 0.0, 0.0, 0.0])
+        lmp.command("run 5 post no")
+        self.assertAlmostEqual(lmp.get_thermo("temp"),1.0/30.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("pe"),1.0/8.0,14)
+        self.assertAlmostEqual(lmp.get_thermo("press"),0.15416666666666667,14)
+        if NUMPY_INSTALLED:
+            npforce = lmp.numpy.fix_external_get_force("ext")
+            self.assertEqual(npforce[0][0],6.0)
+            self.assertEqual(npforce[3][1],6.0)
+            self.assertEqual(npforce[7][2],6.0)
+
+##############################
+if __name__ == "__main__":
+    unittest.main()