Merge branch 'develop' into collected-small-changes

doc/src/Library_utility.rst

@@ -7,6 +7,7 @@ functions.  They do not directly call the LAMMPS library.
 - :cpp:func:`lammps_encode_image_flags`
 - :cpp:func:`lammps_decode_image_flags`
 - :cpp:func:`lammps_set_fix_external_callback`
+- :cpp:func:`lammps_fix_external_get_force`
 - :cpp:func:`lammps_fix_external_set_energy_global`
 - :cpp:func:`lammps_fix_external_set_energy_peratom`
 - :cpp:func:`lammps_fix_external_set_virial_global`
@@ -44,6 +45,11 @@ where such memory buffers were allocated that require the use of
 
 -----------------------
 
+.. doxygenfunction:: lammps_fix_external_get_force
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_fix_external_set_energy_global
    :project: progguide
 

src/library.cpp

@@ -685,7 +685,7 @@ void lammps_commands_string(void *handle, const char *str)
 \verbatim embed:rst
 
 This number may be very large when running large simulations across
-multiple processors.  Depending on compile time choices, LAMMPS may be
+multiple processes.  Depending on compile time choices, LAMMPS may be
 using either 32-bit or a 64-bit integer to store this number. For
 portability this function returns thus a double precision
 floating point number, which can represent up to a 53-bit signed
@@ -2263,7 +2263,7 @@ int lammps_set_variable(void *handle, char *name, char *str)
 // Library functions for scatter/gather operations of data
 // ----------------------------------------------------------------------
 
-/** Gather the named atom-based entity for all atoms across all processors,
+/** Gather the named atom-based entity for all atoms across all processes,
  * in order.
  *
 \verbatim embed:rst
@@ -2282,6 +2282,8 @@ x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], :math:`\dots`);
 *natoms*), as queried by :cpp:func:`lammps_get_natoms`,
 :cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`.
 
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
 \endverbatim
  *
  * \param handle  pointer to a previously created LAMMPS instance
@@ -2301,7 +2303,8 @@ x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], :math:`\dots`);
      Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather_atoms(void *handle, char *name, int type, int count, void *data)
+void lammps_gather_atoms(void *handle, const char *name, int type, int count,
+                         void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -2415,7 +2418,7 @@ void lammps_gather_atoms(void *handle, char *name, int type, int count, void *da
   END_CAPTURE
 }
 
-/** Gather the named atom-based entity for all atoms across all processors,
+/** Gather the named atom-based entity for all atoms across all processes,
  * unordered.
  *
 \verbatim embed:rst
@@ -2432,12 +2435,14 @@ of atoms, use :cpp:func:`lammps_gather_atoms_subset`.
 
 The *data* array will be in groups of *count* values, with *natoms*
 groups total, but not in order by atom ID (e.g., if *name* is *x* and *count*
-is 3, then *data* might be something like = x[10][0], x[10][1], x[10][2],
+is 3, then *data* might be something like x[10][0], x[10][1], x[10][2],
 x[2][0], x[2][1], x[2][2], x[4][0], :math:`\dots`); *data* must be
 pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as
 queried by :cpp:func:`lammps_get_natoms`,
 :cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`.
 
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
 \endverbatim
  *
  * \param handle: pointer to a previously created LAMMPS instance
@@ -2456,7 +2461,8 @@ queried by :cpp:func:`lammps_get_natoms`,
      Allgather Nlocal atoms from each proc into data
 ------------------------------------------------------------------------- */
 
-void lammps_gather_atoms_concat(void *handle, char *name, int type, int count, void *data)
+void lammps_gather_atoms_concat(void *handle, const char *name, int type,
+                                int count, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -2599,6 +2605,8 @@ x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], :math:`\dots`);
 *data* must be pre-allocated by the caller to length
 (*count* :math:`\times` *ndata*).
 
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
 \endverbatim
  *
  * \param handle: pointer to a previously created LAMMPS instance
@@ -2621,8 +2629,8 @@ x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], :math:`\dots`);
      Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather_atoms_subset(void *handle, char *name, int type, int count,
+void lammps_gather_atoms_subset(void *handle, const char *name, int type,
-                                int ndata, int *ids, void *data)
+                                int count, int ndata, int *ids, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -2745,20 +2753,22 @@ void lammps_gather_atoms_subset(void *handle, char *name, int type, int count,
   END_CAPTURE
 }
 
-/** Scatter the named atom-based entities in *data* to all processors.
+/** Scatter the named atom-based entities in *data* to all processes.
  *
 \verbatim embed:rst
 
 This subroutine takes data stored in a one-dimensional array supplied by the
-user and scatters them to all atoms on all processors. The data must be
+user and scatters them to all atoms on all processes. The data must be
 ordered by atom ID, with the requirement that the IDs be consecutive.
 Use :cpp:func:`lammps_scatter_atoms_subset` to scatter data for some (or all)
 atoms, unordered.
 
 The *data* array needs to be ordered in groups of *count* values, sorted by
 atom ID (e.g., if *name* is *x* and *count* = 3, then
-*data* = x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0],
+*data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0],
-:math:`\dots`); *data* must be of length (*count* :math:`\times` *natoms*).
+:math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*).
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
 
 \endverbatim
  *
@@ -2768,7 +2778,7 @@ atom ID (e.g., if *name* is *x* and *count* = 3, then
  * \param count   number of per-atom values (e.g., 1 for *type* or *charge*,
  *                3 for *x* or *f*); use *count* = 3 with *image* if you have
  *                a single image flag packed into (*x*,*y*,*z*) components.
- * \param data    per-atom values packed in a 1-dimensional array of length
+ * \param data    per-atom values packed in a one-dimensional array of length
  *                *natoms* \* *count*.
  *
  */
@@ -2778,7 +2788,8 @@ atom ID (e.g., if *name* is *x* and *count* = 3, then
      loop over Natoms, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
 
-void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *data)
+void lammps_scatter_atoms(void *handle, const char *name, int type, int count,
+                          void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -2879,12 +2890,12 @@ void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *d
 }
 
 /** Scatter the named atom-based entities in *data* from a subset of atoms
- *  to all processors.
+ *  to all processes.
  *
 \verbatim embed:rst
 
 This subroutine takes data stored in a one-dimensional array supplied by the
-user and scatters them to a subset of atoms on all processors. The array
+user and scatters them to a subset of atoms on all processes. The array
 *data* contains data associated with atom IDs, but there is no requirement that
 the IDs be consecutive, as they are provided in a separate array.
 Use :cpp:func:`lammps_scatter_atoms` to scatter data for all atoms, in order.
@@ -2895,6 +2906,8 @@ to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2],
 x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would
 be {1, 100, 57}.
 
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
 \endverbatim
  *
  * \param handle: pointer to a previously created LAMMPS instance
@@ -2928,8 +2941,8 @@ be {1, 100, 57}.
      loop over Ndata, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
 
-void lammps_scatter_atoms_subset(void *handle, char *name, int type, int count,
+void lammps_scatter_atoms_subset(void *handle, const char *name, int type,
-                                 int ndata, int *ids, void *data)
+                                 int count, int ndata, int *ids, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -3146,6 +3159,43 @@ void lammps_gather_bonds(void *handle, void *data)
   END_CAPTURE
 }
 
+/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities
+ *  from all processes, in order by atom ID.
+ *
+\verbatim embed:rst
+
+This subroutine gathers data from all processes and stores them in a one-dimensional array
+allocated by the user. The array *data* will be ordered by atom ID, which requires consecutive IDs
+(1 to *natoms*\ ). If you need a similar array but for non-consecutive atom IDs, see
+:cpp:func:`lammps_gather_concat`; for a similar array but for a subset of atoms, see
+:cpp:func:`lammps_gather_subset`.
+
+The *data* array will be ordered in groups of *count* values, sorted by atom ID (e.g., if *name* is
+*x*, then *data* is {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0],
+:math:`\dots`}); *data* must be pre-allocated by the caller to the correct length
+(*count*\ :math:`{}\times{}`\ *natoms*), as queried by :cpp:func:`lammps_get_natoms`,
+:cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`.
+
+This function will return an error if fix or compute data are requested and the fix or compute ID
+given does not have per-atom data.
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+\endverbatim
+ *
+ * \param handle  pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix
+ *                data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix
+ *                property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix
+ *                property/atom vectors with *count* > 1)
+ * \param type    0 for ``int`` values, 1 for ``double`` values
+ * \param count   number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*);
+ *                use *count* = 3 with *image* if you want the image flags unpacked into
+ *                (*x*,*y*,*z*) components.
+ * \param data    per-atom values packed into a one-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
   Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France)
   gather the named atom-based entity for all atoms
@@ -3172,7 +3222,7 @@ void lammps_gather_bonds(void *handle, void *data)
     Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather(void *handle, char *name, int type, int count, void *data)
+void lammps_gather(void *handle, const char *name, int type, int count, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -3381,6 +3431,44 @@ void lammps_gather(void *handle, char *name, int type, int count, void *data)
   END_CAPTURE
 }
 
+/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities
+ *  from all processes, unordered.
+ *
+\verbatim embed:rst
+
+This subroutine gathers data for all atoms and stores them in a one-dimensional array allocated by
+the user. The data will be a concatenation of chunks from each processor's owned atoms, in
+whatever order the atoms are in on each processor. This process has no requirement that the atom
+IDs be consecutive. If you need the ID of each atom, you can do another call to either
+:cpp:func:`lammps_gather_atoms_concat` or :cpp:func:`lammps_gather_concat` with *name* set to
+``id``. If you have consecutive IDs and want the data to be in order, use
+:cpp:func:`lammps_gather`; for a similar array but for a subset of atoms, use
+:cpp:func:`lammps_gather_subset`.
+
+The *data* array will be in groups of *count* values, with *natoms* groups total, but not in order
+by atom ID (e.g., if *name* is *x* and *count* is 3, then *data* might be something like
+{x[10][0], x[10][1], x[10][2], x[2][0], x[2][1], x[2][2], x[4][0], :math:`\dots`}); *data* must be
+pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as queried by
+:cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`, or
+:cpp:func:`lammps_extract_setting`.
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name:   desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix
+ *                data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix
+ *                property/atom vectors with count = 1, "d2_name" or "i2_name" for fix
+ *                property/atom vectors with count > 1)
+ * \param type:   0 for ``int`` values, 1 for ``double`` values
+ * \param count:  number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*);
+ *                use *count* = 3 with *image* if you want the image flags unpacked into
+ *                (*x*,*y*,*z*) components.
+ * \param data:   per-atom values packed into a one-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
   Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France)
   gather the named atom-based entity for all atoms
@@ -3407,7 +3495,8 @@ void lammps_gather(void *handle, char *name, int type, int count, void *data)
     Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather_concat(void *handle, char *name, int type, int count, void *data)
+void lammps_gather_concat(void *handle, const char *name, int type, int count,
+                          void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -3633,6 +3722,41 @@ void lammps_gather_concat(void *handle, char *name, int type, int count, void *d
   END_CAPTURE
 }
 
+/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities
+ *  from all processes for a subset of atoms.
+ *
+\verbatim embed:rst
+
+This subroutine gathers data for the requested atom IDs and stores them in a one-dimensional array
+allocated by the user. The data will be ordered by atom ID, but there is no requirement that the
+IDs be consecutive. If you wish to return a similar array for *all* the atoms, use
+:cpp:func:`lammps_gather` or :cpp:func:`lammps_gather_concat`.
+
+The *data* array will be in groups of *count* values, sorted by atom ID in the same order as the
+array *ids* (e.g., if *name* is *x*, *count* = 3, and *ids* is {100, 57, 210}, then *data* might
+look like {x[100][0], x[100][1], x[100][2], x[57][0], x[57][1], x[57][2], x[210][0],
+:math:`\dots`}); *ids* must be provided by the user with length *ndata*, and *data* must be
+pre-allocated by the caller to length (*count*\ :math:`{}\times{}`\ *ndata*).
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix
+ *                data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix
+ *                property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix
+ *                property/atom vectors with *count* > 1)
+ * \param type    0 for ``int`` values, 1 for ``double`` values
+ * \param count   number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*);
+ *                use *count* = 3 with *image* if you want the image flags unpacked into
+ *                (*x*,*y*,*z*) components.
+ * \param ndata:  number of atoms for which to return data (can be all of them)
+ * \param ids:    list of *ndata* atom IDs for which to return data
+ * \param data    per-atom values packed into a one-dimensional array of length
+ *                *ndata* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
   Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France)
   gather the named atom-based entity for all atoms
@@ -3659,9 +3783,8 @@ void lammps_gather_concat(void *handle, char *name, int type, int count, void *d
     Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather_subset(void *handle, char *name,
+void lammps_gather_subset(void *handle, const char *name, int type, int count,
-                                int type, int count,
+                          int ndata, int *ids, void *data)
-                                int ndata, int *ids, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -3884,6 +4007,37 @@ void lammps_gather_subset(void *handle, char *name,
   END_CAPTURE
 }
 
+/** Scatter the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based
+ *  entity in *data* to all processes.
+ *
+\verbatim embed:rst
+
+This subroutine takes data stored in a one-dimensional array supplied by the user and scatters
+them to all atoms on all processes. The data must be ordered by atom ID, with the requirement that
+the IDs be consecutive. Use :cpp:func:`lammps_scatter_subset` to scatter data for some (or all)
+atoms, unordered.
+
+The *data* array needs to be ordered in groups of *count* values, sorted by atom ID (e.g., if
+*name* is *x* and *count* = 3, then *data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1],
+x[1][2], x[2][0], :math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*).
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+\endverbatim
+ *
+ * \param handle  pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix
+ *                data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix
+ *                property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix
+ *                property/atom vectors with *count* > 1)
+ * \param type    0 for ``int`` values, 1 for ``double`` values
+ * \param count   number of per-atom values (e.g., 1 for *type* or *charge*,
+ *                3 for *x* or *f*); use *count* = 3 with *image* if you have
+ *                a single image flag packed into (*x*,*y*,*z*) components.
+ * \param data    per-atom values packed in a one-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
   Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France)
   scatter the named atom-based entity in data to all atoms
@@ -3908,7 +4062,8 @@ void lammps_gather_subset(void *handle, char *name,
     Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_scatter(void *handle, char *name, int type, int count, void *data)
+void lammps_scatter(void *handle, const char *name, int type, int count,
+                    void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -4103,6 +4258,42 @@ void lammps_scatter(void *handle, char *name, int type, int count, void *data)
   END_CAPTURE
 }
 
+/** Scatter the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based
+ * entities in *data* from a subset of atoms to all processes.
+ *
+\verbatim embed:rst
+
+This subroutine takes data stored in a one-dimensional array supplied by the
+user and scatters them to a subset of atoms on all processes. The array
+*data* contains data associated with atom IDs, but there is no requirement that
+the IDs be consecutive, as they are provided in a separate array.
+Use :cpp:func:`lammps_scatter` to scatter data for all atoms, in order.
+
+The *data* array needs to be organized in groups of *count* values, with the
+groups in the same order as the array *ids*. For example, if you want *data*
+to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2],
+x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would
+be {1, 100, 57}.
+
+This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix
+ *                data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix
+ *                property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix
+ *                property/atom vectors with *count* > 1)
+ * \param type:   0 for ``int`` values, 1 for ``double`` values
+ * \param count:  number of per-atom values (e.g., 1 for *type* or *charge*,
+ *                3 for *x* or *f*); use *count* = 3 with "image" if you want
+ *                single image flags unpacked into (*x*,*y*,*z*)
+ * \param ndata:  number of atoms listed in *ids* and *data* arrays
+ * \param ids:    list of *ndata* atom IDs to scatter data to
+ * \param data    per-atom values packed in a 1-dimensional array of length
+ *                *ndata* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
   Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France)
    scatter the named atom-based entity in data to a subset of atoms
@@ -4125,7 +4316,7 @@ void lammps_scatter(void *handle, char *name, int type, int count, void *data)
      loop over Ndata, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
 
-void lammps_scatter_subset(void *handle, char *name,int type, int count,
+void lammps_scatter_subset(void *handle, const char *name,int type, int count,
                                  int ndata, int *ids, void *data)
 {
   auto lmp = (LAMMPS *) handle;

src/library.h

@@ -181,23 +181,23 @@ int lammps_set_variable(void *, char *, char *);
  * Library functions for scatter/gather operations of data
  * ---------------------------------------------------------------------- */
 
-void lammps_gather_atoms(void *handle, char *name, int type, int count, void *data);
+void lammps_gather_atoms(void *handle, const char *name, int type, int count, void *data);
-void lammps_gather_atoms_concat(void *handle, char *name, int type, int count, void *data);
+void lammps_gather_atoms_concat(void *handle, const char *name, int type, int count, void *data);
-void lammps_gather_atoms_subset(void *handle, char *name, int type, int count, int ndata, int *ids,
+void lammps_gather_atoms_subset(void *handle, const char *name, int type, int count, int ndata,
-                                void *data);
+                                int *ids, void *data);
-void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *data);
+void lammps_scatter_atoms(void *handle, const char *name, int type, int count, void *data);
-void lammps_scatter_atoms_subset(void *handle, char *name, int type, int count, int ndata, int *ids,
+void lammps_scatter_atoms_subset(void *handle, const char *name, int type, int count, int ndata,
-                                 void *data);
+                                 int *ids, void *data);
 
 void lammps_gather_bonds(void *handle, void *data);
 
-void lammps_gather(void *handle, char *name, int type, int count, void *data);
+void lammps_gather(void *handle, const char *name, int type, int count, void *data);
-void lammps_gather_concat(void *handle, char *name, int type, int count, void *data);
+void lammps_gather_concat(void *handle, const char *name, int type, int count, void *data);
-void lammps_gather_subset(void *handle, char *name, int type, int count, int ndata, int *ids,
+void lammps_gather_subset(void *handle, const char *name, int type, int count, int ndata,
-                          void *data);
+                          int *ids, void *data);
-void lammps_scatter(void *handle, char *name, int type, int count, void *data);
+void lammps_scatter(void *handle, const char *name, int type, int count, void *data);
-void lammps_scatter_subset(void *handle, char *name, int type, int count, int ndata, int *ids,
+void lammps_scatter_subset(void *handle, const char *name, int type, int count, int ndata,
-                           void *data);
+                           int *ids, void *data);
 
 #if !defined(LAMMPS_BIGBIG)
 int lammps_create_atoms(void *handle, int n, const int *id, const int *type, const double *x,

unittest/fortran/CMakeLists.txt

@@ -86,6 +86,10 @@ if(CMAKE_Fortran_COMPILER)
   target_link_libraries(test_fortran_neighlist PRIVATE flammps lammps GTest::GMockMain)
   add_test(NAME FortranNeighlist COMMAND test_fortran_neighlist)
 
+  add_executable(test_fortran_fixexternal wrap_fixexternal.cpp test_fortran_fixexternal.f90)
+  target_link_libraries(test_fortran_fixexternal PRIVATE flammps lammps GTest::GMockMain)
+  add_test(NAME FortranFixExternal COMMAND test_fortran_fixexternal)
+
 else()
   message(STATUS "Skipping Tests for the LAMMPS Fortran Module: no Fortran compiler")
 endif()

unittest/fortran/test_fortran_fixexternal.f90

@@ -0,0 +1,424 @@
+MODULE ext_stuff
+  USE, INTRINSIC :: ISO_Fortran_ENV, ONLY : error_unit
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int, c_int64_t, c_loc
+  USE LIBLAMMPS
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER :: vec_length = 8
+  REAL(c_double), SAVE :: direction = 1.0_c_double
+  REAL(c_double), DIMENSION(:,:), POINTER, SAVE :: f3 => NULL(), f4 => NULL()
+
+CONTAINS
+
+  SUBROUTINE f_lammps_reverse_direction() BIND(C)
+    direction = -direction
+  END SUBROUTINE f_lammps_reverse_direction
+
+  SUBROUTINE f_callback_ss(instance, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT) :: instance
+    INTEGER(c_int) :: timestep
+    INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+    REAL(c_double), DIMENSION(SIZE(id)) :: e
+    REAL(c_double), DIMENSION(6,SIZE(id)) :: v
+
+    WHERE (id == 1)
+      f(1,:) = 1.0_c_double
+      f(2,:) = -1.0_c_double
+      f(3,:) = 1.25_c_double
+      e = 1.0_c_double
+      v(1,:) = 1.0_c_double
+      v(2,:) = 2.0_c_double
+      v(3,:) = -1.0_c_double
+      v(4,:) = -2.0_c_double
+      v(5,:) = 3.0_c_double
+      v(6,:) = -3.0_c_double
+    ELSEWHERE
+      f(1,:) = -1.0_c_double
+      f(2,:) = +1.0_c_double
+      f(3,:) = -1.25_c_double
+      e = 10.0_c_double
+      v(1,:) = 10.0_c_double
+      v(2,:) = 20.0_c_double
+      v(3,:) = -10.0_c_double
+      v(4,:) = -20.0_c_double
+      v(5,:) = 30.0_c_double
+      v(6,:) = -30.0_c_double
+    END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+        CALL instance%fix_external_set_virial_peratom('ext1', v)
+      CLASS DEFAULT
+        WRITE(error_unit,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback_ss
+
+  SUBROUTINE f_callback_sb(instance, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT) :: instance
+    INTEGER(c_int64_t) :: timestep
+    INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+    REAL(c_double), DIMENSION(SIZE(id)) :: e
+    REAL(c_double), DIMENSION(6,SIZE(id)) :: v
+
+    WHERE (id == 1_c_int)
+      f(1,:) = 1.0_c_double
+      f(2,:) = -1.0_c_double
+      f(3,:) = 1.25_c_double
+      e = 1.0_c_double
+      v(1,:) = 1.0_c_double
+      v(2,:) = 2.0_c_double
+      v(3,:) = -1.0_c_double
+      v(4,:) = -2.0_c_double
+      v(5,:) = 3.0_c_double
+      v(6,:) = -3.0_c_double
+    ELSEWHERE
+      f(1,:) = -1.0_c_double
+      f(2,:) = +1.0_c_double
+      f(3,:) = -1.25_c_double
+      e = 10.0_c_double
+      v(1,:) = 10.0_c_double
+      v(2,:) = 20.0_c_double
+      v(3,:) = -10.0_c_double
+      v(4,:) = -20.0_c_double
+      v(5,:) = 30.0_c_double
+      v(6,:) = -30.0_c_double
+    END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+        CALL instance%fix_external_set_virial_peratom('ext1', v)
+      CLASS DEFAULT
+        WRITE(error_unit,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback_sb
+
+  SUBROUTINE f_callback_bb(instance, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT) :: instance
+    INTEGER(c_int64_t) :: timestep
+    INTEGER(c_int64_t), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+    REAL(c_double), DIMENSION(SIZE(id)) :: e
+    REAL(c_double), DIMENSION(6,SIZE(id)) :: v
+
+    WHERE (id == 1_c_int64_t)
+      f(1,:) = 1.0_c_double
+      f(2,:) = -1.0_c_double
+      f(3,:) = 1.25_c_double
+      e = 1.0_c_double
+      v(1,:) = 1.0_c_double
+      v(2,:) = 2.0_c_double
+      v(3,:) = -1.0_c_double
+      v(4,:) = -2.0_c_double
+      v(5,:) = 3.0_c_double
+      v(6,:) = -3.0_c_double
+    ELSEWHERE
+      f(1,:) = -1.0_c_double
+      f(2,:) = +1.0_c_double
+      f(3,:) = -1.25_c_double
+      e = 10.0_c_double
+      v(1,:) = 10.0_c_double
+      v(2,:) = 20.0_c_double
+      v(3,:) = -10.0_c_double
+      v(4,:) = -20.0_c_double
+      v(5,:) = 30.0_c_double
+      v(6,:) = -30.0_c_double
+    END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+        CALL instance%fix_external_set_virial_peratom('ext1', v)
+      CLASS DEFAULT
+        WRITE(error_unit,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback_bb
+
+  SUBROUTINE f_callback2_ss(entity, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT), target :: entity
+    INTEGER(c_int) :: timestep
+    INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+
+    SELECT TYPE (entity)
+      TYPE IS (REAL(c_double))
+        WHERE (id == 1_c_int)
+          f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double
+        ELSEWHERE
+          f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double)
+        END WHERE
+      CLASS DEFAULT
+        WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in&
+          & f_callback2_ss'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback2_ss
+
+  SUBROUTINE f_callback2_sb(entity, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT), target :: entity
+    INTEGER(c_int64_t) :: timestep
+    INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+
+    SELECT TYPE (entity)
+      TYPE IS (REAL(c_double))
+        WHERE (id == 1_c_int)
+          f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double
+        ELSEWHERE
+          f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double)
+        END WHERE
+      CLASS DEFAULT
+        WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in&
+          & f_callback2_sb'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback2_sb
+
+  SUBROUTINE f_callback2_bb(entity, timestep, id, x, f)
+    CLASS(*), INTENT(INOUT), target :: entity
+    INTEGER(c_int64_t) :: timestep
+    INTEGER(c_int64_t), DIMENSION(:), INTENT(IN) :: id
+    REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x
+    REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f
+
+    SELECT TYPE (entity)
+      TYPE IS (REAL(c_double))
+        WHERE (id == 1_c_int64_t)
+          f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double
+        ELSEWHERE
+          f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double)
+          f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double
+          f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double)
+        END WHERE
+      CLASS DEFAULT
+        WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in&
+          & f_callback2_sb'
+        STOP 1
+    END SELECT
+  END SUBROUTINE f_callback2_bb
+END MODULE ext_stuff
+
+FUNCTION f_lammps_with_args() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_ptr
+  USE liblammps
+  USE keepstuff, ONLY: lmp
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_with_args
+
+  CHARACTER(len=12), DIMENSION(12), PARAMETER :: args = &
+      [ CHARACTER(len=12) :: 'liblammps', '-log', 'none', &
+      '-echo','screen','-nocite','-var','zpos','1.5','-var','x','2']
+
+  lmp = lammps(args)
+  f_lammps_with_args = lmp%handle
+END FUNCTION f_lammps_with_args
+
+SUBROUTINE f_lammps_close() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_null_ptr
+  USE liblammps
+  USE keepstuff, ONLY: lmp
+  IMPLICIT NONE
+
+  CALL lmp%close()
+  lmp%handle = c_null_ptr
+END SUBROUTINE f_lammps_close
+
+SUBROUTINE f_lammps_setup_fix_external_callback() BIND(C)
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp, demo_input, cont_input, pair_input
+  USE ext_stuff, ONLY : vec_length
+  IMPLICIT NONE
+
+  CALL lmp%commands_list(demo_input)
+  CALL lmp%commands_list(cont_input)
+  CALL lmp%commands_list(pair_input)
+  CALL lmp%command('neigh_modify exclude group all all')
+  CALL lmp%command('fix ext1 all external pf/callback 1 1')
+  CALL lmp%command('fix ext2 all external pf/callback 1 1')
+  CALL lmp%fix_external_set_vector_length('ext2', vec_length)
+END SUBROUTINE f_lammps_setup_fix_external_callback
+
+SUBROUTINE f_lammps_setup_fix_external_array() BIND(C)
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp, demo_input, cont_input, pair_input
+  USE ext_stuff, ONLY : f3, f4
+  IMPLICIT NONE
+
+  CALL lmp%commands_list(demo_input)
+  CALL lmp%commands_list(cont_input)
+  CALL lmp%commands_list(pair_input)
+  CALL lmp%command('neigh_modify exclude group all all')
+  CALL lmp%command('fix ext3 all external pf/array 1')
+  CALL lmp%command('fix ext4 all external pf/array 1')
+  CALL lmp%command('thermo_style custom step pxx pe etotal')
+  CALL lmp%command('thermo_modify norm no')
+  CALL lmp%command('thermo 100')
+  f3 = lmp%fix_external_get_force('ext3')
+  f4 = lmp%fix_external_get_force('ext4')
+END SUBROUTINE f_lammps_setup_fix_external_array
+
+SUBROUTINE f_lammps_set_fix_external_callbacks() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE ext_stuff
+  IMPLICIT NONE
+  INTEGER :: size_bigint, size_tagint, nlocal
+
+  nlocal = lmp%extract_setting('nlocal')
+
+  size_bigint = lmp%extract_setting('bigint')
+  size_tagint = lmp%extract_setting('tagint')
+  IF (size_bigint == 4_c_int .AND. size_tagint == 4_c_int) THEN
+    CALL lmp%set_fix_external_callback('ext1', f_callback_ss, lmp)
+    CALL lmp%set_fix_external_callback('ext2', f_callback2_ss, direction)
+  ELSE IF (size_bigint == 8_c_int .AND. size_tagint == 8_c_int) THEN
+    CALL lmp%set_fix_external_callback('ext1', f_callback_bb, lmp)
+    CALL lmp%set_fix_external_callback('ext2', f_callback2_bb, direction)
+  ELSE
+    CALL lmp%set_fix_external_callback('ext1', f_callback_sb, lmp)
+    CALL lmp%set_fix_external_callback('ext2', f_callback2_sb, direction)
+  END IF
+END SUBROUTINE f_lammps_set_fix_external_callbacks
+
+SUBROUTINE f_lammps_get_force (i, ptr) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double, c_ptr, C_F_POINTER
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  TYPE(c_ptr), INTENT(IN), VALUE :: ptr
+  REAL(c_double), DIMENSION(:,:), POINTER :: force => NULL()
+  REAL(c_double), DIMENSION(:), POINTER :: f => NULL()
+
+  CALL C_F_POINTER(ptr, f, [3])
+  force = lmp%extract_atom('f')
+  f = force(:,i)
+END SUBROUTINE f_lammps_get_force
+
+SUBROUTINE f_lammps_find_forces() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int, c_int64_t
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE ext_stuff, ONLY : f3, f4
+  IMPLICIT NONE
+  INTEGER(c_int) :: size_tagint
+  INTEGER(c_int), DIMENSION(:), POINTER :: id
+  INTEGER(c_int64_t), DIMENSION(:), POINTER :: tag
+
+  f3(:,:) = 0.0_c_double
+  f4(:,:) = 0.0_c_double
+  size_tagint = lmp%extract_setting('tagint')
+  IF (size_tagint == 4_c_int) THEN
+    id = lmp%extract_atom('id')
+    WHERE (id == 1_c_int)
+      f3(1,:) = 4.0_c_double
+      f3(2,:) = -4.0_c_double
+      f3(3,:) = 6.0_c_double
+      f4(1,:) = 10.0_c_double
+      f4(2,:) = -10.0_c_double
+      f4(3,:) = 12.0_c_double
+    ELSEWHERE
+      f3(1,:) = 5.0_c_double
+      f3(2,:) = -5.0_c_double
+      f3(3,:) = 7.0_c_double
+      f4(1,:) = 11.0_c_double
+      f4(2,:) = -11.0_c_double
+      f4(3,:) = 13.0_c_double
+    END WHERE
+  ELSE
+    tag = lmp%extract_atom('id')
+    WHERE (tag == 1_c_int64_t)
+      f3(1,:) = 4.0_c_double
+      f3(2,:) = -4.0_c_double
+      f3(3,:) = 6.0_c_double
+      f4(1,:) = 10.0_c_double
+      f4(2,:) = -10.0_c_double
+      f4(3,:) = 12.0_c_double
+    ELSEWHERE
+      f3(1,:) = 5.0_c_double
+      f3(2,:) = -5.0_c_double
+      f3(3,:) = 7.0_c_double
+      f4(1,:) = 11.0_c_double
+      f4(2,:) = -11.0_c_double
+      f4(3,:) = 13.0_c_double
+    END WHERE
+  END IF
+END SUBROUTINE f_lammps_find_forces
+
+SUBROUTINE f_lammps_add_energy() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+
+  CALL lmp%fix_external_set_energy_global('ext3', -20.2_c_double);
+END SUBROUTINE f_lammps_add_energy
+
+SUBROUTINE f_lammps_set_virial() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+
+  CALL lmp%fix_external_set_virial_global('ext4', [1.0_c_double, &
+    2.0_c_double, 2.5_c_double, -1.0_c_double, -2.25_c_double, -3.02_c_double])
+END SUBROUTINE f_lammps_set_virial
+
+FUNCTION f_lammps_find_peratom_energy(i) RESULT(energy) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double) :: energy
+  REAL(c_double), DIMENSION(:), POINTER :: e
+
+  e = lmp%extract_compute('peratom', lmp%style%atom, lmp%type%vector)
+  energy = e(i)
+END FUNCTION f_lammps_Find_peratom_energy
+
+SUBROUTINE f_lammps_find_peratom_virial(v, i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double), DIMENSION(6) :: v
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double), DIMENSION(:,:), POINTER :: virial
+
+  virial = lmp%extract_compute('vperatom', lmp%style%atom, lmp%type%array)
+  v = virial(:,i)
+END SUBROUTINE f_lammps_find_peratom_virial
+
+SUBROUTINE f_lammps_fixexternal_set_vector() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE ext_stuff, ONLY : vec_length
+  IMPLICIT NONE
+  REAL(c_double), DIMENSION(vec_length) :: v
+  INTEGER :: i
+  DO i = 1, vec_length
+    v(i) = REAL(i, c_double)
+    CALL lmp%fix_external_set_vector('ext2', i, v(i))
+  END DO
+END SUBROUTINE f_lammps_fixexternal_set_vector

unittest/fortran/test_fortran_gather_scatter.f90

@@ -24,13 +24,16 @@ END SUBROUTINE f_lammps_close
 
 SUBROUTINE f_lammps_setup_gather_scatter() BIND(C)
   USE LIBLAMMPS
-  USE keepstuff, ONLY : lmp, big_input, cont_input, more_input
+  USE keepstuff, ONLY : lmp, big_input, cont_input, more_input, pair_input
   IMPLICIT NONE
 
   CALL lmp%command('atom_modify map array')
   CALL lmp%commands_list(big_input)
   CALL lmp%commands_list(cont_input)
   CALL lmp%commands_list(more_input)
+  CALL lmp%commands_list(pair_input)
+  CALL lmp%command('mass 1 1.0')
+  CALL lmp%command("compute pe all pe/atom")
 END SUBROUTINE f_lammps_setup_gather_scatter
 
 FUNCTION f_lammps_gather_atoms_mask(i) BIND(C)
@@ -262,3 +265,90 @@ FUNCTION f_lammps_test_gather_bonds_big() BIND(C) RESULT(success)
     success = 0_c_int
   END IF
 END FUNCTION f_lammps_test_gather_bonds_big
+
+FUNCTION f_lammps_gather_pe_atom(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double) :: f_lammps_gather_pe_atom
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom
+
+  CALL lmp%gather('c_pe', 1_c_int, pe_atom)
+  f_lammps_gather_pe_atom = pe_atom(i)
+END FUNCTION f_lammps_gather_pe_atom
+
+FUNCTION f_lammps_gather_pe_atom_concat(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double) :: f_lammps_gather_pe_atom_concat
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: tag
+  INTEGER :: j
+
+  CALL lmp%gather_concat('id', 1_c_int, tag)
+  CALL lmp%gather_concat('c_pe', 1_c_int, pe_atom)
+  DO j = 1, SIZE(tag)
+    IF (tag(j) == i) THEN
+       f_lammps_gather_pe_atom_concat = pe_atom(j)
+       EXIT
+    END IF
+  END DO
+  f_lammps_gather_pe_atom_concat = pe_atom(i)
+END FUNCTION f_lammps_gather_pe_atom_concat
+
+SUBROUTINE f_lammps_gather_pe_atom_subset(ids, pe) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN) :: ids(2)
+  REAL(c_double), INTENT(OUT) :: pe(2)
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom
+  INTEGER(c_int) :: natoms
+
+  natoms = NINT(lmp%get_natoms(), c_int)
+  CALL lmp%gather_subset('c_pe', 1, ids, pe_atom)
+  pe(1:natoms) = pe_atom
+END SUBROUTINE f_lammps_gather_pe_atom_subset
+
+SUBROUTINE f_lammps_scatter_compute() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom
+  REAL(c_double) :: swap
+
+  CALL lmp%gather('c_pe', 1_c_int, pe_atom)
+
+  ! swap the computed energy of atoms 1 and 3
+  swap = pe_atom(1)
+  pe_atom(1) = pe_atom(3)
+  pe_atom(3) = swap
+
+  CALL lmp%scatter('c_pe', pe_atom) ! push the swap back to LAMMPS
+END SUBROUTINE f_lammps_scatter_compute
+
+SUBROUTINE f_lammps_scatter_subset_compute() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), PARAMETER :: ids(2) = [3,1]
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom
+  REAL(c_double) :: swap
+
+  CALL lmp%gather_subset('c_pe', 1_c_int, ids, pe_atom)
+
+  ! swap the computed energy of atoms 1 and 3
+  swap = pe_atom(1)
+  pe_atom(1) = pe_atom(2)
+  pe_atom(2) = swap
+
+  CALL lmp%scatter_subset('c_pe', ids, pe_atom) ! push the swap back to LAMMPS
+END SUBROUTINE f_lammps_scatter_subset_compute

unittest/fortran/wrap_fixexternal.cpp

@@ -0,0 +1,194 @@
+
+// unit tests for gathering and scattering data from a LAMMPS instance through
+// the Fortran wrapper
+
+#include "lammps.h"
+#include "library.h"
+#include <cstdint>
+#include <cstdlib>
+#include <mpi.h>
+#include <string>
+
+#include "gtest/gtest.h"
+
+// prototypes for Fortran reverse wrapper functions
+extern "C" {
+void *f_lammps_with_args();
+void f_lammps_close();
+void f_lammps_setup_fix_external_callback();
+void f_lammps_setup_fix_external_array();
+void f_lammps_set_fix_external_callbacks();
+void f_lammps_get_force(int, double*);
+void f_lammps_reverse_direction();
+void f_lammps_find_forces();
+void f_lammps_add_energy();
+void f_lammps_set_virial();
+double f_lammps_find_peratom_energy(int);
+void f_lammps_find_peratom_virial(double[6], int);
+void f_lammps_fixexternal_set_vector();
+}
+
+using namespace LAMMPS_NS;
+
+class LAMMPS_fixexternal : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_fixexternal()           = default;
+    ~LAMMPS_fixexternal() override = default;
+
+    void SetUp() override
+    {
+        ::testing::internal::CaptureStdout();
+        lmp                = (LAMMPS_NS::LAMMPS *)f_lammps_with_args();
+        std::string output = ::testing::internal::GetCapturedStdout();
+        EXPECT_STREQ(output.substr(0, 8).c_str(), "LAMMPS (");
+    }
+    void TearDown() override
+    {
+        ::testing::internal::CaptureStdout();
+        f_lammps_close();
+        std::string output = ::testing::internal::GetCapturedStdout();
+        EXPECT_STREQ(output.substr(0, 16).c_str(), "Total wall time:");
+        lmp = nullptr;
+    }
+};
+
+TEST_F(LAMMPS_fixexternal, callback)
+{
+    f_lammps_setup_fix_external_callback();
+    f_lammps_set_fix_external_callbacks();
+    lammps_command(lmp, "run 0");
+    double f[3];
+    f_lammps_get_force(1,f);
+    EXPECT_DOUBLE_EQ(f[0], 3.0);
+    EXPECT_DOUBLE_EQ(f[1], -3.0);
+    EXPECT_DOUBLE_EQ(f[2], 3.75);
+    f_lammps_get_force(2,f);
+    EXPECT_DOUBLE_EQ(f[0], -3.0);
+    EXPECT_DOUBLE_EQ(f[1], 3.0);
+    EXPECT_DOUBLE_EQ(f[2], -3.75);
+
+    f_lammps_reverse_direction();
+    f_lammps_set_fix_external_callbacks();
+    lammps_command(lmp, "run 0");
+    f_lammps_get_force(1,f);
+    EXPECT_DOUBLE_EQ(f[0], -1.0);
+    EXPECT_DOUBLE_EQ(f[1], 1.0);
+    EXPECT_DOUBLE_EQ(f[2], -1.25);
+    f_lammps_get_force(2,f);
+    EXPECT_DOUBLE_EQ(f[0], 1.0);
+    EXPECT_DOUBLE_EQ(f[1], -1.0);
+    EXPECT_DOUBLE_EQ(f[2], 1.25);
+};
+
+TEST_F(LAMMPS_fixexternal, array)
+{
+    f_lammps_setup_fix_external_array();
+    double **f;
+    f = (double**) lammps_extract_atom(lmp, "f");
+    f_lammps_find_forces();
+    lammps_command(lmp, "run 0");
+    EXPECT_DOUBLE_EQ(f[0][0], 14.0);
+    EXPECT_DOUBLE_EQ(f[0][1], -14.0);
+    EXPECT_DOUBLE_EQ(f[0][2], 18.0);
+    EXPECT_DOUBLE_EQ(f[1][0], 16.0);
+    EXPECT_DOUBLE_EQ(f[1][1], -16.0);
+    EXPECT_DOUBLE_EQ(f[1][2], 20.0);
+};
+
+TEST_F(LAMMPS_fixexternal, energy_global)
+{
+    f_lammps_setup_fix_external_array();
+    double energy;
+    f_lammps_add_energy();
+    lammps_command(lmp, "run 0");
+    energy = lammps_get_thermo(lmp, "etotal");
+    EXPECT_DOUBLE_EQ(energy, -20.2);
+};
+
+TEST_F(LAMMPS_fixexternal, virial_global)
+{
+    f_lammps_setup_fix_external_array();
+    double virial[6], volume;
+    f_lammps_set_virial();
+    lammps_command(lmp, "run 0");
+    volume = lammps_get_thermo(lmp, "vol");
+    virial[0] = lammps_get_thermo(lmp, "pxx");
+    virial[1] = lammps_get_thermo(lmp, "pyy");
+    virial[2] = lammps_get_thermo(lmp, "pzz");
+    virial[3] = lammps_get_thermo(lmp, "pxy");
+    virial[4] = lammps_get_thermo(lmp, "pxz");
+    virial[5] = lammps_get_thermo(lmp, "pyz");
+    EXPECT_DOUBLE_EQ(virial[0], 1.0/volume);
+    EXPECT_DOUBLE_EQ(virial[1], 2.0/volume);
+    EXPECT_DOUBLE_EQ(virial[2], 2.5/volume);
+    EXPECT_DOUBLE_EQ(virial[3], -1.0/volume);
+    EXPECT_DOUBLE_EQ(virial[4], -2.25/volume);
+    EXPECT_DOUBLE_EQ(virial[5], -3.02/volume);
+};
+
+TEST_F(LAMMPS_fixexternal, energy_peratom)
+{
+    f_lammps_setup_fix_external_callback();
+    f_lammps_set_fix_external_callbacks();
+    lammps_command(lmp, "compute peratom all pe/atom");
+    double energy;
+    lammps_command(lmp, "run 0");
+    int nlocal = lammps_extract_setting(lmp, "nlocal");
+    for (int i = 1; i <= nlocal; i++)
+    {
+        energy = f_lammps_find_peratom_energy(i);
+        if (i == 1)
+          EXPECT_DOUBLE_EQ(energy, 1.0);
+        else
+          EXPECT_DOUBLE_EQ(energy, 10.0);
+    }
+};
+
+TEST_F(LAMMPS_fixexternal, virial_peratom)
+{
+    f_lammps_setup_fix_external_callback();
+    f_lammps_set_fix_external_callbacks();
+    lammps_command(lmp, "compute vperatom all stress/atom NULL");
+    double virial[6];
+    lammps_command(lmp, "run 0");
+    int nlocal = lammps_extract_setting(lmp, "nlocal");
+    for (int i = 1; i <= nlocal; i++)
+    {
+        f_lammps_find_peratom_virial(virial, i);
+        if (i == 1)
+        {
+          EXPECT_DOUBLE_EQ(virial[0], -1.0);
+          EXPECT_DOUBLE_EQ(virial[1], -2.0);
+          EXPECT_DOUBLE_EQ(virial[2], 1.0);
+          EXPECT_DOUBLE_EQ(virial[3], 2.0);
+          EXPECT_DOUBLE_EQ(virial[4], -3.0);
+          EXPECT_DOUBLE_EQ(virial[5], 3.0);
+        }
+        else
+        {
+          EXPECT_DOUBLE_EQ(virial[0], -10.0);
+          EXPECT_DOUBLE_EQ(virial[1], -20.0);
+          EXPECT_DOUBLE_EQ(virial[2], 10.0);
+          EXPECT_DOUBLE_EQ(virial[3], 20.0);
+          EXPECT_DOUBLE_EQ(virial[4], -30.0);
+          EXPECT_DOUBLE_EQ(virial[5], 30.0);
+        }
+    }
+};
+
+TEST_F(LAMMPS_fixexternal, vector)
+{
+    f_lammps_setup_fix_external_callback();
+    f_lammps_set_fix_external_callbacks();
+    f_lammps_fixexternal_set_vector();
+    lammps_command(lmp, "run 0");
+    double *v;
+    for (int i = 0; i < 8; i++)
+    {
+      v = (double*) lammps_extract_fix(lmp, "ext2", LMP_STYLE_GLOBAL,
+        LMP_TYPE_VECTOR, i, 1);
+      EXPECT_DOUBLE_EQ(i+1, *v);
+      std::free(v);
+    }
+};

unittest/fortran/wrap_gather_scatter.cpp

@@ -3,6 +3,7 @@
 
 #include "lammps.h"
 #include "library.h"
+#include "atom.h"
 #include <cstdint>
 #include <cstdlib>
 #include <mpi.h>
@@ -26,6 +27,11 @@ void f_lammps_scatter_atoms_positions();
 void f_lammps_setup_gather_bonds();
 int f_lammps_test_gather_bonds_small();
 int f_lammps_test_gather_bonds_big();
+double f_lammps_gather_pe_atom(int);
+double f_lammps_gather_pe_atom_concat(int);
+void f_lammps_gather_pe_atom_subset(int*, double*);
+void f_lammps_scatter_compute();
+void f_lammps_scatter_subset_compute();
 }
 
 using namespace LAMMPS_NS;
@@ -216,3 +222,113 @@ TEST_F(LAMMPS_gather_scatter, gather_bonds)
     EXPECT_EQ(f_lammps_test_gather_bonds_small(), 1);
 #endif
 };
+
+TEST_F(LAMMPS_gather_scatter, gather_compute)
+{
+#ifdef LAMMPS_BIGBIG
+    GTEST_SKIP();
+#else
+    f_lammps_setup_gather_scatter();
+    lammps_command(lmp, "run 0");
+    int natoms = lmp->atom->natoms;
+    int *tag = lmp->atom->tag;
+    double *pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM,
+        LMP_TYPE_VECTOR);
+    for (int i = 0; i < natoms; i++)
+        EXPECT_DOUBLE_EQ(f_lammps_gather_pe_atom(tag[i]), pe[i]);
+#endif
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_compute_concat)
+{
+#ifdef LAMMPS_BIGBIG
+    GTEST_SKIP();
+#else
+    f_lammps_setup_gather_scatter();
+    lammps_command(lmp, "run 0");
+    int natoms = lmp->atom->natoms;
+    int *tag = lmp->atom->tag;
+    double *pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM,
+        LMP_TYPE_VECTOR);
+    for (int i = 0; i < natoms; i++)
+        EXPECT_DOUBLE_EQ(f_lammps_gather_pe_atom(tag[i]), pe[i]);
+#endif
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_compute_subset)
+{
+#ifdef LAMMPS_BIGBIG
+    GTEST_SKIP();
+#else
+    f_lammps_setup_gather_scatter();
+    lammps_command(lmp, "run 0");
+    int ids[2] = {3, 1};
+    int *tag = lmp->atom->tag;
+    double pe[2] = {0.0, 0.0};
+    int nlocal = lammps_extract_setting(lmp, "nlocal");
+    double *pa_pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM,
+      LMP_TYPE_VECTOR);
+
+    for (int i = 0; i < nlocal; i++) {
+      if(tag[i] == ids[0]) pe[0] = pa_pe[i];
+      if(tag[i] == ids[1]) pe[1] = pa_pe[i];
+    }
+
+    double ftn_pe[2];
+    f_lammps_gather_pe_atom_subset(ids, ftn_pe);
+    EXPECT_DOUBLE_EQ(ftn_pe[0], pe[0]);
+    EXPECT_DOUBLE_EQ(ftn_pe[1], pe[1]);
+#endif
+};
+
+TEST_F(LAMMPS_gather_scatter, scatter_compute)
+{
+#ifdef LAMMPS_BIGBIG
+    GTEST_SKIP();
+#else
+    f_lammps_setup_gather_scatter();
+    int natoms = lmp->atom->natoms;
+    double *pe = new double[natoms];
+    lammps_command(lmp, "run 0");
+    lammps_gather(lmp, "c_pe", 1, 1, pe);
+    double *old_pe = new double[natoms];
+    for (int i = 0; i < natoms; i++)
+      old_pe[i] = pe[i];
+    EXPECT_DOUBLE_EQ(pe[0], old_pe[0]);
+    EXPECT_DOUBLE_EQ(pe[1], old_pe[1]);
+    EXPECT_DOUBLE_EQ(pe[2], old_pe[2]);
+    f_lammps_scatter_compute();
+    lammps_gather(lmp, "c_pe", 1, 1, pe);
+    EXPECT_DOUBLE_EQ(pe[0], old_pe[2]);
+    EXPECT_DOUBLE_EQ(pe[1], old_pe[1]);
+    EXPECT_DOUBLE_EQ(pe[2], old_pe[0]);
+    delete[] old_pe;
+    delete[] pe;
+#endif
+};
+
+TEST_F(LAMMPS_gather_scatter, scatter_subset_compute)
+{
+#ifdef LAMMPS_BIGBIG
+    GTEST_SKIP();
+#else
+    f_lammps_setup_gather_scatter();
+    int natoms = lmp->atom->natoms;
+    double *pe = new double[natoms];
+    lammps_command(lmp, "run 0");
+    lammps_gather(lmp, "c_pe", 1, 1, pe);
+    double *old_pe = new double[natoms];
+    for (int i = 0; i < natoms; i++)
+      old_pe[i] = pe[i];
+    EXPECT_DOUBLE_EQ(pe[0], old_pe[0]);
+    EXPECT_DOUBLE_EQ(pe[1], old_pe[1]);
+    EXPECT_DOUBLE_EQ(pe[2], old_pe[2]);
+    f_lammps_scatter_subset_compute();
+    lammps_gather(lmp, "c_pe", 1, 1, pe);
+    EXPECT_DOUBLE_EQ(pe[0], old_pe[2]);
+    EXPECT_DOUBLE_EQ(pe[1], old_pe[1]);
+    EXPECT_DOUBLE_EQ(pe[2], old_pe[0]);
+    delete[] old_pe;
+    delete[] pe;
+#endif
+};