Merge pull request #3470 from hammondkd/fortran-further-tinkering

Further extending the Fortran interface

doc/src/Build_manual.rst

@@ -216,7 +216,7 @@ be multiple tests run automatically:
 - A test that only standard, printable ASCII text characters are used.
   This runs the command ``env LC_ALL=C grep -n '[^ -~]' src/*.rst`` and
   thus prints all offending lines with filename and line number
-  prepended to the screen.  Special characters like greek letters
+  prepended to the screen.  Special characters like Greek letters
   (:math:`\alpha~~\sigma~~\epsilon`), super- or subscripts
   (:math:`x^2~~\mathrm{U}_{LJ}`), mathematical expressions
   (:math:`\frac{1}{2}\mathrm{N}~~x\to\infty`), or the Angstrom symbol

doc/src/Library_objects.rst

@@ -6,6 +6,7 @@ fixes, or variables in LAMMPS using the following functions:
 
 - :cpp:func:`lammps_extract_compute`
 - :cpp:func:`lammps_extract_fix`
+- :cpp:func:`lammps_extract_variable_datatype`
 - :cpp:func:`lammps_extract_variable`
 - :cpp:func:`lammps_set_variable`
 
@@ -21,6 +22,11 @@ fixes, or variables in LAMMPS using the following functions:
 
 -----------------------
 
+.. doxygenfunction:: lammps_extract_variable_datatype
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_extract_variable
    :project: progguide
 
@@ -36,3 +42,5 @@ fixes, or variables in LAMMPS using the following functions:
 .. doxygenenum:: _LMP_STYLE_CONST
 
 .. doxygenenum:: _LMP_TYPE_CONST
+
+.. doxygenenum:: _LMP_VAR_CONST

doc/src/Library_properties.rst

@@ -16,8 +16,8 @@ This section documents the following functions:
 --------------------
 
 The library interface allows the extraction of different kinds of
-information about the active simulation instance and also - in some
+information about the active simulation instance and also---in some
-cases - to apply modifications to it.  This enables combining of a
+cases---to apply modifications to it.  This enables combining of a
 LAMMPS simulation with other processing and simulation methods computed
 by the calling code, or by another code that is coupled to LAMMPS via
 the library interface.  In some cases the data returned is direct
@@ -25,9 +25,9 @@ reference to the original data inside LAMMPS, cast to a void pointer.
 In that case the data needs to be cast to a suitable pointer for the
 calling program to access it, and you may need to know the correct
 dimensions and lengths.  This also means you can directly change those
-value(s) from the calling program, e.g. to modify atom positions.  Of
+value(s) from the calling program (e.g., to modify atom positions).  Of
-course, this should be done with care.  When accessing per-atom data,
+course, changing values should be done with care.  When accessing per-atom
-please note that this data is the per-processor **local** data and is
+data, please note that these data are the per-processor **local** data and are
 indexed accordingly. Per-atom data can change sizes and ordering at
 every neighbor list rebuild or atom sort event as atoms migrate between
 sub-domains and processors.

doc/src/Run_basics.rst

@@ -30,12 +30,13 @@ executable itself can be placed elsewhere.
 
 .. note::
 
-   The redirection operator "<" will not always work when running
+   The redirection operator "<" will not always work when running in
-   in parallel with mpirun or mpiexec; for those systems the -in form is required.
+   parallel with ``mpirun`` or ``mpiexec``; for those systems the -in
+   form is required.
 
 As LAMMPS runs it prints info to the screen and a logfile named
-*log.lammps*\ .  More info about output is given on the
+*log.lammps*\ .  More info about output is given on the :doc:`screen and
-:doc:`screen and logfile output <Run_output>` page.
+logfile output <Run_output>` page.
 
 If LAMMPS encounters errors in the input script or while running a
 simulation it will print an ERROR message and stop or a WARNING

doc/src/Run_options.rst

@@ -93,13 +93,13 @@ switch is not set (the default), LAMMPS will operate as if the KOKKOS
 package were not installed; i.e. you can run standard LAMMPS or with
 the GPU or OPENMP packages, for testing or benchmarking purposes.
 
-Additional optional keyword/value pairs can be specified which
+Additional optional keyword/value pairs can be specified which determine
-determine how Kokkos will use the underlying hardware on your
+how Kokkos will use the underlying hardware on your platform.  These
-platform.  These settings apply to each MPI task you launch via the
+settings apply to each MPI task you launch via the ``mpirun`` or
-"mpirun" or "mpiexec" command.  You may choose to run one or more MPI
+``mpiexec`` command.  You may choose to run one or more MPI tasks per
-tasks per physical node.  Note that if you are running on a desktop
+physical node.  Note that if you are running on a desktop machine, you
-machine, you typically have one physical node.  On a cluster or
+typically have one physical node.  On a cluster or supercomputer there
-supercomputer there may be dozens or 1000s of physical nodes.
+may be dozens or 1000s of physical nodes.
 
 Either the full word or an abbreviation can be used for the keywords.
 Note that the keywords do not use a leading minus sign.  I.e. the
@@ -148,9 +148,9 @@ one of these 4 environment variables
    MV2_COMM_WORLD_LOCAL_RANK (Mvapich)
    OMPI_COMM_WORLD_LOCAL_RANK (OpenMPI)
 
-which are initialized by the "srun", "mpirun" or "mpiexec" commands.
+which are initialized by the ``srun``, ``mpirun``, or ``mpiexec``
-The environment variable setting for each MPI rank is used to assign a
+commands.  The environment variable setting for each MPI rank is used to
-unique GPU ID to the MPI task.
+assign a unique GPU ID to the MPI task.
 
 .. parsed-literal::
 

doc/src/Speed_gpu.rst

@@ -76,10 +76,11 @@ instructions.
 
 **Run with the GPU package from the command line:**
 
-The mpirun or mpiexec command sets the total number of MPI tasks used
+The ``mpirun`` or ``mpiexec`` command sets the total number of MPI tasks
-by LAMMPS (one or multiple per compute node) and the number of MPI
+used by LAMMPS (one or multiple per compute node) and the number of MPI
-tasks used per node.  E.g. the mpirun command in MPICH does this via
+tasks used per node.  E.g. the ``mpirun`` command in MPICH does this via
-its -np and -ppn switches.  Ditto for OpenMPI via -np and -npernode.
+its ``-np`` and ``-ppn`` switches.  Ditto for OpenMPI via ``-np`` and
+``-npernode``.
 
 When using the GPU package, you cannot assign more than one GPU to a
 single MPI task.  However multiple MPI tasks can share the same GPU,
@@ -129,8 +130,8 @@ GPU package pair styles.
 
 **Or run with the GPU package by editing an input script:**
 
-The discussion above for the mpirun/mpiexec command, MPI tasks/node,
+The discussion above for the ``mpirun`` or ``mpiexec`` command, MPI
-and use of multiple MPI tasks/GPU is the same.
+tasks/node, and use of multiple MPI tasks/GPU is the same.
 
 Use the :doc:`suffix gpu <suffix>` command, or you can explicitly add an
 "gpu" suffix to individual styles in your input script, e.g.

doc/src/Speed_kokkos.rst

@@ -72,12 +72,12 @@ See the :ref:`Build extras <kokkos>` page for instructions.
 Running LAMMPS with the KOKKOS package
 """"""""""""""""""""""""""""""""""""""
 
-All Kokkos operations occur within the context of an individual MPI
+All Kokkos operations occur within the context of an individual MPI task
-task running on a single node of the machine. The total number of MPI
+running on a single node of the machine. The total number of MPI tasks
-tasks used by LAMMPS (one or multiple per compute node) is set in the
+used by LAMMPS (one or multiple per compute node) is set in the usual
-usual manner via the mpirun or mpiexec commands, and is independent of
+manner via the ``mpirun`` or ``mpiexec`` commands, and is independent of
-Kokkos. E.g. the mpirun command in OpenMPI does this via its -np and
+Kokkos. E.g. the mpirun command in OpenMPI does this via its ``-np`` and
--npernode switches. Ditto for MPICH via -np and -ppn.
+``-npernode`` switches. Ditto for MPICH via ``-np`` and ``-ppn``.
 
 Running on a multi-core CPU
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -310,7 +310,8 @@ Alternatively the effect of the "-sf" or "-pk" switches can be
 duplicated by adding the :doc:`package kokkos <package>` or :doc:`suffix kk <suffix>` commands to your input script.
 
 The discussion above for building LAMMPS with the KOKKOS package, the
-mpirun/mpiexec command, and setting appropriate thread are the same.
+``mpirun`` or ``mpiexec`` command, and setting appropriate thread
+properties are the same.
 
 You must still use the "-k on" :doc:`command-line switch <Run_options>`
 to enable the KOKKOS package, and specify its additional arguments for

doc/src/Speed_omp.rst

@@ -33,8 +33,8 @@ These examples assume one or more 16-core nodes.
    mpirun -np 4 lmp_omp -sf omp -pk omp 4 -in in.script           # 4 MPI tasks, 4 threads/task
    mpirun -np 32 -ppn 4 lmp_omp -sf omp -pk omp 4 -in in.script   # 8 nodes, 4 MPI tasks/node, 4 threads/task
 
-The mpirun or mpiexec command sets the total number of MPI tasks used
+The ``mpirun`` or ``mpiexec`` command sets the total number of MPI tasks
-by LAMMPS (one or multiple per compute node) and the number of MPI
+used by LAMMPS (one or multiple per compute node) and the number of MPI
 tasks used per node.  E.g. the mpirun command in MPICH does this via
 its -np and -ppn switches.  Ditto for OpenMPI via -np and -npernode.
 
@@ -58,8 +58,8 @@ OMP_NUM_THREADS environment variable.
 Or run with the OPENMP package by editing an input script
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-The discussion above for the mpirun/mpiexec command, MPI tasks/node,
+The discussion above for the ``mpirun`` or ``mpiexec`` command, MPI
-and threads/MPI task is the same.
+tasks/node, and threads/MPI task is the same.
 
 Use the :doc:`suffix omp <suffix>` command, or you can explicitly add an
 "omp" suffix to individual styles in your input script, e.g.

doc/utils/sphinx-config/false_positives.txt

@@ -78,6 +78,7 @@ Alexey
 ali
 aliceblue
 Allinger
+allocatable
 allocator
 allocators
 allosws
@@ -657,6 +658,7 @@ Dcut
 de
 dE
 De
+deallocate
 deallocated
 debye
 Debye
@@ -691,6 +693,7 @@ dequidt
 Dequidt
 der
 dereference
+dereferenced
 derekt
 Deresiewicz
 Derjagin
@@ -1486,6 +1489,7 @@ interfacial
 interial
 interlayer
 intermolecular
+interoperable
 Interparticle
 interstitials
 intertube
@@ -3621,6 +3625,7 @@ Universite
 unix
 unmaintained
 unoptimized
+unordered
 unpadded
 unphysical
 unphysically

fortran/README

@@ -1,9 +1,9 @@
-This directory contains Fortran code which interface LAMMPS as a library
+This directory contains Fortran code that acts as an interface to LAMMPS as a
-and allows the LAMMPS library interface to be invoked from Fortran codes.
+library and allows the LAMMPS library interface to be invoked from Fortran
-It requires a Fortran compiler that supports the Fortran 2003 standard.
+code.  It requires a Fortran compiler that supports the Fortran 2003 standard.
 
 This interface is based on and supersedes the previous Fortran interfaces
-in the examples/COUPLE/fortran* folders, but is fully supported by the
+in the examples/COUPLE/fortran* folders, but it is fully supported by the
 LAMMPS developers and included in the documentation and unit testing.
 
 Details on this Fortran interface and how to build programs using it

python/lammps/constants.py

@@ -22,7 +22,8 @@ LAMMPS_INT64       = 4
 LAMMPS_INT64_2D    = 5
 LAMMPS_STRING      = 6
 
-# these must be kept in sync with the enums in library.h
+# these must be kept in sync with the enums in src/library.h, tools/swig/lammps.i
+# and the constants in fortran/lammps.f90
 LMP_STYLE_GLOBAL   = 0
 LMP_STYLE_ATOM     = 1
 LMP_STYLE_LOCAL    = 2
@@ -42,6 +43,8 @@ LMP_ERROR_UNIVERSE = 8
 
 LMP_VAR_EQUAL      = 0
 LMP_VAR_ATOM       = 1
+LMP_VAR_VECTOR     = 2
+LMP_VAR_STRING     = 3
 
 # -------------------------------------------------------------------------
 

python/lammps/core.py

@@ -301,6 +301,8 @@ class lammps(object):
     self.lib.lammps_extract_fix.argtypes = [c_void_p, c_char_p, c_int, c_int, c_int, c_int]
 
     self.lib.lammps_extract_variable.argtypes = [c_void_p, c_char_p, c_char_p]
+    self.lib.lammps_extract_variable_datatype.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_extract_variable_datatype.restype = c_int
 
     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))
@@ -1083,21 +1085,23 @@ class lammps(object):
   # for vector, must copy nlocal returned values to local c_double vector
   # memory was allocated by library interface function
 
-  def extract_variable(self, name, group=None, vartype=LMP_VAR_EQUAL):
+  def extract_variable(self, name, group=None, vartype=None):
     """ Evaluate a LAMMPS variable and return its data
 
     This function is a wrapper around the function
-    :cpp:func:`lammps_extract_variable` of the C-library interface,
+    :cpp:func:`lammps_extract_variable` of the C library interface,
     evaluates variable name and returns a copy of the computed data.
     The memory temporarily allocated by the C-interface is deleted
     after the data is copied to a Python variable or list.
     The variable must be either an equal-style (or equivalent)
-    variable or an atom-style variable. The variable type has to
+    variable or an atom-style variable. The variable type can be
-    provided as ``vartype`` parameter which may be one of two constants:
+    provided as the ``vartype`` parameter, which may be one of several
-    ``LMP_VAR_EQUAL`` or ``LMP_VAR_ATOM``; it defaults to
+    constants: ``LMP_VAR_EQUAL``, ``LMP_VAR_ATOM``, ``LMP_VAR_VECTOR``,
-    equal-style variables.
+    or ``LMP_VAR_STRING``. If omitted or ``None``, LAMMPS will determine its
-    The group parameter is only used for atom-style variables and
+    value for you based on a call to
-    defaults to the group "all" if set to ``None``, which is the default.
+    :cpp:func:`lammps_extract_variable_datatype` from the C library interface.
+    The group parameter is only used for atom-style variables and defaults to
+    the group "all".
 
     :param name: name of the variable to execute
     :type name: string
@@ -1111,6 +1115,8 @@ class lammps(object):
     if name: name = name.encode()
     else: return None
     if group: group = group.encode()
+    if vartype is None :
+      vartype = self.lib.lammps_extract_variable_datatype(self.lmp, name)
     if vartype == LMP_VAR_EQUAL:
       self.lib.lammps_extract_variable.restype = POINTER(c_double)
       with ExceptionCheck(self):
@@ -1130,6 +1136,31 @@ class lammps(object):
         self.lib.lammps_free(ptr)
       else: return None
       return result
+    elif vartype == LMP_VAR_VECTOR :
+      nvector = 0
+      self.lib.lammps_extract_variable.restype = POINTER(c_int)
+      ptr = self.lib.lammps_extract_variable(self.lmp,name,
+              'LMP_SIZE_VECTOR'.encode())
+      if ptr :
+        nvector = ptr[0]
+        self.lib.lammps_free(ptr)
+      else :
+        return None
+      self.lib.lammps_extract_variable.restype = POINTER(c_double)
+      result = (c_double*nvector)()
+      values = self.lib.lammps_extract_variable(self.lmp,name,group)
+      if values :
+        for i in range(nvector) :
+          result[i] = values[i]
+        # do NOT free the values pointer (points to internal vector data)
+        return result
+      else :
+        return None
+    elif vartype == LMP_VAR_STRING :
+      self.lib.lammps_extract_variable.restype = c_char_p
+      with ExceptionCheck(self) :
+        ptr = self.lib.lammps_extract_variable(self.lmp, name, group)
+        return ptr.decode('utf-8')
     return None
 
   # -------------------------------------------------------------------------

src/library.cpp

@@ -949,7 +949,8 @@ int lammps_get_mpi_comm(void *handle)
 This function will retrieve or compute global properties. In contrast to
 :cpp:func:`lammps_get_thermo` this function returns an ``int``.  The
 following tables list the currently supported keyword.  If a keyword is
-not recognized, the function returns -1.
+not recognized, the function returns -1.  The integer sizes functions may
+be called without a valid LAMMPS object handle (it is ignored).
 
 * :ref:`Integer sizes <extract_integer_sizes>`
 * :ref:`System status <extract_system_status>`
@@ -1145,7 +1146,7 @@ int lammps_extract_setting(void *handle, const char *keyword)
 This function returns an integer that encodes the data type of the global
 property with the specified name. See :cpp:enum:`_LMP_DATATYPE_CONST` for valid
 values. Callers of :cpp:func:`lammps_extract_global` can use this information
-to then decide how to cast the (void*) pointer and access the data.
+to then decide how to cast the ``void *`` pointer and access the data.
 
 .. versionadded:: 18Sep2020
 
@@ -1622,7 +1623,7 @@ void *lammps_extract_global(void *handle, const char *name)
 This function returns an integer that encodes the data type of the per-atom
 property with the specified name. See :cpp:enum:`_LMP_DATATYPE_CONST` for valid
 values. Callers of :cpp:func:`lammps_extract_atom` can use this information
-to then decide how to cast the (void*) pointer and access the data.
+to then decide how to cast the ``void *`` pointer and access the data.
 
 .. versionadded:: 18Sep2020
 
@@ -2043,16 +2044,19 @@ void *lammps_extract_fix(void *handle, const char *id, int style, int type,
 
 This function returns a pointer to data from a LAMMPS :doc:`variable`
 identified by its name.  When the variable is either an *equal*\ -style
-compatible or an *atom*\ -style variable the variable is evaluated and
+compatible variable, a *vector*\ -style variable, or an *atom*\ -style
-the corresponding value(s) returned.  Variables of style *internal*
+variable, the variable is evaluated and the corresponding value(s) returned.
-are compatible with *equal*\ -style variables and so are *python*\
+Variables of style *internal* are compatible with *equal*\ -style variables and
--style variables, if they return a numeric value.  For other
+so are *python*\ -style variables, if they return a numeric value.  For other
-variable styles their string value is returned.  The function returns
+variable styles, their string value is returned.  The function returns
 ``NULL`` when a variable of the provided *name* is not found or of an
 incompatible style.  The *group* argument is only used for *atom*\
--style variables and ignored otherwise.  If set to ``NULL`` when
+-style variables and ignored otherwise, with one exception: for style *vector*,
-extracting data from and *atom*\ -style variable, the group is assumed
+if *group* is "GET_VECTOR_SIZE", the returned pointer will yield the length
-to be "all".
+of the vector to be returned when dereferenced. This pointer must be
+deallocated after the value is read to avoid a memory leak.
+If *group* is set to ``NULL`` when extracting data from an *atom*\ -style
+variable, the group is assumed to be "all".
 
 When requesting data from an *equal*\ -style or compatible variable
 this function allocates storage for a single double value, copies the
@@ -2066,15 +2070,23 @@ use to avoid a memory leak. Example:
    double value = *dptr;
    lammps_free((void *)dptr);
 
-For *atom*\ -style variables the data returned is a pointer to an
+For *atom*\ -style variables, the return value is a pointer to an
 allocated block of storage of double of the length ``atom->nlocal``.
-Since the data is returned a copy, the location will persist, but its
+Since the data returned are a copy, the location will persist, but its
-content will not be updated, in case the variable is re-evaluated.
+content will not be updated in case the variable is re-evaluated.
-To avoid a memory leak this pointer needs to be freed after use in
+To avoid a memory leak, this pointer needs to be freed after use in
 the calling program.
 
+For *vector*\ -style variables, the returned pointer is to actual LAMMPS data.
+The pointer should not be deallocated. Its length depends on the variable,
+compute, or fix data used to construct the *vector*\ -style variable.
+This length can be fetched by calling this function with *group* set to the
+constant "LMP_SIZE_VECTOR", which returns a ``void *`` pointer that can be
+dereferenced to an integer that is the length of the vector. This pointer
+needs to be deallocated when finished with it to avoid memory leaks.
+
 For other variable styles the returned pointer needs to be cast to
-a char pointer.
+a char pointer. It should not be deallocated.
 
 .. code-block:: c
 
@@ -2084,10 +2096,10 @@ a char pointer.
 .. note::
 
    LAMMPS cannot easily check if it is valid to access the data
-   referenced by the variables, e.g. computes or fixes or thermodynamic
+   referenced by the variables (e.g., computes, fixes, or thermodynamic
-   info, so it may fail with an error.  The caller has to make certain,
+   info), so it may fail with an error.  The caller has to make certain
-   that the data is extracted only when it safe to evaluate the variable
+   that the data are extracted only when it safe to evaluate the variable
-   and thus an error and crash is avoided.
+   and thus an error or crash are avoided.
 
 \endverbatim
  *
@@ -2118,6 +2130,15 @@ void *lammps_extract_variable(void *handle, const char *name, const char *group)
       auto vector = (double *) malloc(nlocal*sizeof(double));
       lmp->input->variable->compute_atom(ivar,igroup,vector,1,0);
       return (void *) vector;
+    } else if (lmp->input->variable->vectorstyle(ivar)) {
+      double *values = nullptr;
+      int nvector = lmp->input->variable->compute_vector(ivar, &values);
+      if ( group != nullptr && strcmp(group,"LMP_SIZE_VECTOR") == 0 ) {
+          int* nvecptr = (int *) malloc(sizeof(int));
+          *nvecptr = nvector;
+          return (void *) nvecptr;
+      } else
+        return (void *) values;
     } else {
       return lmp->input->variable->retrieve(name);
     }
@@ -2130,6 +2151,49 @@ void *lammps_extract_variable(void *handle, const char *name, const char *group)
 
 /* ---------------------------------------------------------------------- */
 
+/** Get data type of a LAMMPS variable.
+ *
+\verbatim embed:rst
+
+This function returns an integer that encodes the data type of the variable
+with the specified name. See :cpp:enum:`_LMP_VAR_CONST` for valid values.
+Callers of :cpp:func:`lammps_extract_variable` can use this information to
+decide how to cast the ``void *`` pointer and access the data.
+
+.. versionadded:: TBD
+
+\endverbatim
+ *
+ * \param  handle  pointer to a previously created LAMMPS instance
+ * \param  name    string with the name of the extracted variable
+ * \return         integer constant encoding the data type of the property
+ *                 or -1 if not found.
+ **/
+
+int lammps_extract_variable_datatype(void *handle, const char *name)
+{
+  auto lmp = (LAMMPS*) handle;
+
+  BEGIN_CAPTURE
+  {
+    int ivar = lmp->input->variable->find(name);
+    if ( ivar < 0 ) return -1;
+
+    if (lmp->input->variable->equalstyle(ivar))
+      return LMP_VAR_EQUAL;
+    else if (lmp->input->variable->atomstyle(ivar))
+      return LMP_VAR_ATOM;
+    else if (lmp->input->variable->vectorstyle(ivar))
+      return LMP_VAR_VECTOR;
+    else
+      return LMP_VAR_STRING;
+  }
+  END_CAPTURE
+  return -1;
+}
+
+/* ---------------------------------------------------------------------- */
+
 /** Set the value of a string-style variable.
  *
  * This function assigns a new value from the string str to the
@@ -2159,21 +2223,38 @@ 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,
+ * in order.
+ *
+\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 ordered by
+atom ID, which requires consecutive atom IDs (1 to *natoms*\ ). If you need
+a similar array but have non-consecutive atom IDs, see
+:cpp:func:`lammps_gather_atoms_concat`; for a similar array but for a subset
+of atoms, see :cpp:func:`lammps_gather_atoms_subset`.
+
+The *data* array will 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 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`.
+
+\endverbatim
+ *
+ * \param handle  pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., *x* or *charge*)
+ * \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
+ *                a single image flag unpacked into (*x*,*y*,*z*) components.
+ * \param data    per-atom values packed in a 1-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
 /* ----------------------------------------------------------------------
-   gather the named atom-based entity for all atoms
-     return it in user-allocated data
-   data will be ordered by atom ID
-     requirement for consecutive atom IDs (1 to N)
-   see gather_atoms_concat() to return data for all atoms, unordered
-   see gather_atoms_subset() to return data for only a subset of atoms
-   name = desired quantity, e.g. x or charge
-   type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
-     use count = 3 with "image" if want single image flag unpacked into xyz
-   return atom-based values in 1d data, ordered by count, then by atom ID
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
-     data must be pre-allocated by caller to correct length
-     correct length = count*Natoms, as queried by get_natoms()
    method:
      alloc and zero count*Natom length vector
      loop over Nlocal to fill vector with my values
@@ -2294,23 +2375,43 @@ 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,
+ * 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
+:cpp:func:`lammps_gather_atoms_concat` call with *name* set to ``id``.
+If you have consecutive IDs and want the data to be in order, use
+:cpp:func:`lammps_gather_atoms`; for a similar array but for a subset
+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],
+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`.
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name:   desired quantity (e.g., *x* or *charge*\ )
+ * \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 data:   per-atom values packed in a 1-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
+
 /* ----------------------------------------------------------------------
-   gather the named atom-based entity for all atoms
-     return it in user-allocated data
-   data will be a concatenation of chunks of each proc's atoms,
-     in whatever order the atoms are on each proc
-     no requirement for consecutive atom IDs (1 to N)
-     can do a gather_atoms_concat for "id" if need to know atom IDs
-   see gather_atoms() to return data ordered by consecutive atom IDs
-   see gather_atoms_subset() to return data for only a subset of atoms
-   name = desired quantity, e.g. x or charge
-   type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
-     use count = 3 with "image" if want single image flag unpacked into xyz
-   return atom-based values in 1d data, ordered by count, then by atom
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
-     data must be pre-allocated by caller to correct length
-     correct length = count*Natoms, as queried by get_natoms()
    method:
      Allgather Nlocal atoms from each proc into data
 ------------------------------------------------------------------------- */
@@ -2440,23 +2541,40 @@ void lammps_gather_atoms_concat(void *handle, char *name, int type, int count, v
   END_CAPTURE
 }
 
+/** Gather the named atom-based entity 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_atoms`
+or :cpp:func:`lammps_gather_atoms_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*).
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name:   desired quantity (e.g., *x* or *charge*)
+ * \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 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 in a 1-dimensional array of length
+ *                *ndata* \* *count*.
+ *
+ */
+
 /* ----------------------------------------------------------------------
-   gather the named atom-based entity for a subset of atoms
-     return it in user-allocated data
-   data will be ordered by requested atom IDs
-     no requirement for consecutive atom IDs (1 to N)
-   see gather_atoms() to return data for all atoms, ordered by consecutive IDs
-   see gather_atoms_concat() to return data for all atoms, unordered
-   name = desired quantity, e.g. x or charge
-   type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
-     use count = 3 with "image" if want single image flag unpacked into xyz
-   ndata = # of atoms to return data for (could be all atoms)
-   ids = list of Ndata atom IDs to return data for
-   return atom-based values in 1d data, ordered by count, then by atom
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
-     data must be pre-allocated by caller to correct length
-     correct length = count*Ndata
    method:
      alloc and zero count*Ndata length vector
      loop over Ndata to fill vector with my values
@@ -2477,15 +2595,16 @@ void lammps_gather_atoms_subset(void *handle, char *name, int type, int count,
     int i,j,m,offset;
     tagint id;
 
-    // error if tags are not defined
+    // error if tags are not defined or no atom map
     // NOTE: test that name = image or ids is not a 64-bit int in code?
 
     int flag = 0;
     if (lmp->atom->tag_enable == 0) flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
+    if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
       if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms_subset");
+        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms_subset: atoms must have mappable ids");
       return;
     }
 
@@ -2586,18 +2705,35 @@ 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.
+ *
+\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
+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],
+:math:`\dots`); *data* must be of length (*count* :math:`\times` *natoms*).
+
+\endverbatim
+ *
+ * \param handle  pointer to a previously created LAMMPS instance
+ * \param name    desired quantity (e.g., *x* or *charge*)
+ * \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 1-dimensional array of length
+ *                *natoms* \* *count*.
+ *
+ */
+
 /* ----------------------------------------------------------------------
-   scatter the named atom-based entity in data to all atoms
-   data is ordered by atom ID
-     requirement for consecutive atom IDs (1 to N)
-   see scatter_atoms_subset() to scatter data for some (or all) atoms, unordered
-   name = desired quantity, e.g. x or charge
-   type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
-     use count = 3 with "image" for xyz to be packed into single image flag
-   data = atom-based values in 1d data, ordered by count, then by atom ID
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
-     data must be correct length = count*Natoms, as queried by get_natoms()
    method:
      loop over Natoms, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
@@ -2624,7 +2760,7 @@ void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *d
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
       if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms");
+        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms: ids must exist, be consecutive, and be mapped");
       return;
     }
 
@@ -2702,19 +2838,51 @@ void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *d
   END_CAPTURE
 }
 
+/** Scatter the named atom-based entities in *data* from a subset of atoms
+ *  to all processors.
+ *
+\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
+*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.
+
+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}.
+
+\endverbatim
+ *
+ * \param handle: pointer to a previously created LAMMPS instance
+ * \param name:   desired quantity (e.g., *x* or *charge*)
+ * \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
+ *                all the image flags packed into (*xyz*)
+ * \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*.
+ *
+ */
+
 /* ----------------------------------------------------------------------
    scatter the named atom-based entity in data to a subset of atoms
    data is ordered by provided atom IDs
      no requirement for consecutive atom IDs (1 to N)
    see scatter_atoms() to scatter data for all atoms, ordered by consecutive IDs
-   name = desired quantity, e.g. x or charge
+   name = desired quantity (e.g., x or charge)
    type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+   count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
      use count = 3 with "image" for xyz to be packed into single image flag
    ndata = # of atoms in ids and data (could be all atoms)
    ids = list of Ndata atom IDs to scatter data to
    data = atom-based values in 1d data, ordered by count, then by atom ID
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
      data must be correct length = count*Ndata
    method:
      loop over Ndata, if I own atom ID, set its values from data
@@ -2743,7 +2911,7 @@ void lammps_scatter_atoms_subset(void *handle, char *name, int type, int count,
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
       if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms_subset");
+        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms_subset: atoms must have mapped ids");
       return;
     }
 
@@ -2951,10 +3119,10 @@ void lammps_gather_bonds(void *handle, void *data)
          "d2_name" or "i2_name" for fix property/atom arrays with count > 1
          will return error if fix/compute isn't atom-based
   type = 0 for integer values, 1 for double values
-  count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+  count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
     use count = 3 with "image" if want single image flag unpacked into xyz
   return atom-based values in 1d data, ordered by count, then by atom ID
-    e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+    (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
     data must be pre-allocated by caller to correct length
     correct length = count*Natoms, as queried by get_natoms()
   method:
@@ -3186,10 +3354,10 @@ void lammps_gather(void *handle, char *name, int type, int count, void *data)
          "d2_name" or "i2_name" for fix property/atom arrays with count > 1
          will return error if fix/compute isn't atom-based
   type = 0 for integer values, 1 for double values
-  count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+  count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
     use count = 3 with "image" if want single image flag unpacked into xyz
   return atom-based values in 1d data, ordered by count, then by atom ID
-    e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+    (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
     data must be pre-allocated by caller to correct length
     correct length = count*Natoms, as queried by get_natoms()
   method:
@@ -3438,10 +3606,10 @@ void lammps_gather_concat(void *handle, char *name, int type, int count, void *d
          "d2_name" or "i2_name" for fix property/atom arrays with count > 1
          will return error if fix/compute isn't atom-based
   type = 0 for integer values, 1 for double values
-  count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+  count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
     use count = 3 with "image" if want single image flag unpacked into xyz
   return atom-based values in 1d data, ordered by count, then by atom ID
-    e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+    (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
     data must be pre-allocated by caller to correct length
     correct length = count*Natoms, as queried by get_natoms()
   method:
@@ -3465,11 +3633,12 @@ void lammps_gather_subset(void *handle, char *name,
     int i,j,m,offset,ltype;
     tagint id;
 
-    // error if tags are not defined or not consecutive
+    // error if tags are not defined or no atom map
 
     int flag = 0;
     if (lmp->atom->tag_enable == 0) flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
+    if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
       if (lmp->comm->me == 0)
         lmp->error->warning(FLERR,"Library error in lammps_gather_subset");
@@ -3686,10 +3855,10 @@ void lammps_gather_subset(void *handle, char *name,
          "d2_name" or "i2_name" for fix property/atom arrays with count > 1
          will return error if fix/compute isn't atom-based
   type = 0 for integer values, 1 for double values
-  count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+  count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
     use count = 3 with "image" if want single image flag unpacked into xyz
   return atom-based values in 1d data, ordered by count, then by atom ID
-    e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+    (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
     data must be pre-allocated by caller to correct length
     correct length = count*Natoms, as queried by get_natoms()
   method:
@@ -3904,12 +4073,12 @@ void lammps_scatter(void *handle, char *name, int type, int count, void *data)
           "f_fix", "c_compute" for fixes / computes
           will return error if fix/compute doesn't isn't atom-based
    type = 0 for integer values, 1 for double values
-   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
+   count = # of per-atom values (e.g., 1 for type or charge, 3 for x or f)
      use count = 3 with "image" for xyz to be packed into single image flag
    ndata = # of atoms in ids and data (could be all atoms)
    ids = list of Ndata atom IDs to scatter data to
    data = atom-based values in 1d data, ordered by count, then by atom ID
-     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+     (e.g., x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...)
      data must be correct length = count*Ndata
    method:
      loop over Ndata, if I own atom ID, set its values from data
@@ -4152,8 +4321,8 @@ boundaries atoms will be wrapped back into the simulation cell and its
 image flags adjusted accordingly, unless explicit image flags are
 provided.
 
-The function returns the number of atoms created or -1 on failure, e.g.
+The function returns the number of atoms created or -1 on failure (e.g.,
-when called before as box has been created.
+when called before as box has been created).
 
 Coordinates and velocities have to be given in a 1d-array in the order
 X(1),Y(1),Z(1),X(2),Y(2),Z(2),...,X(N),Y(N),Z(N).
@@ -4511,7 +4680,7 @@ int lammps_config_has_mpi_support()
  * files via a pipe to gzip or similar compression programs
 
 \verbatim embed:rst
-Several LAMMPS commands (e.g. :doc:`read_data`, :doc:`write_data`,
+Several LAMMPS commands (e.g., :doc:`read_data`, :doc:`write_data`,
 :doc:`dump styles atom, custom, and xyz <dump>`) support reading and
 writing compressed files via creating a pipe to the ``gzip`` program.
 This function checks whether this feature was :ref:`enabled at compile
@@ -4531,8 +4700,8 @@ int lammps_config_has_gzip_support() {
 
 \verbatim embed:rst
 The LAMMPS :doc:`dump style image <dump_image>` supports writing multiple
-image file formats.  Most of them need, however, support from an external
+image file formats.  Most of them, however, need support from an external
-library and using that has to be :ref:`enabled at compile time <graphics>`.
+library, and using that has to be :ref:`enabled at compile time <graphics>`.
 This function checks whether support for the `PNG image file format
 <https://en.wikipedia.org/wiki/Portable_Network_Graphics>`_ is available
 in the current LAMMPS library.
@@ -4550,8 +4719,8 @@ int lammps_config_has_png_support() {
 
 \verbatim embed:rst
 The LAMMPS :doc:`dump style image <dump_image>` supports writing multiple
-image file formats.  Most of them need, however, support from an external
+image file formats.  Most of them, however, need support from an external
-library and using that has to be :ref:`enabled at compile time <graphics>`.
+library, and using that has to be :ref:`enabled at compile time <graphics>`.
 This function checks whether support for the `JPEG image file format
 <https://jpeg.org/jpeg/>`_ is available in the current LAMMPS library.
 \endverbatim
@@ -4582,14 +4751,14 @@ int lammps_config_has_ffmpeg_support() {
 
 /* ---------------------------------------------------------------------- */
 
-/** Check whether LAMMPS errors will throw a C++ exception
+/** Check whether LAMMPS errors will throw C++ exceptions.
  *
 \verbatim embed:rst
-In case of errors LAMMPS will either abort or throw a C++ exception.
+In case of an error, LAMMPS will either abort or throw a C++ exception.
 The latter has to be :ref:`enabled at compile time <exceptions>`.
 This function checks if exceptions were enabled.
 
-When using the library interface and C++ exceptions are enabled,
+When using the library interface with C++ exceptions enabled,
 the library interface functions will "catch" them and the
 error status can then be checked by calling
 :cpp:func:`lammps_has_error` and the most recent error message
@@ -4608,10 +4777,10 @@ int lammps_config_has_exceptions() {
 
 /* ---------------------------------------------------------------------- */
 
-/** Check if a specific package has been included in LAMMPS
+/** Check whether a specific package has been included in LAMMPS
  *
 \verbatim embed:rst
-This function checks if the LAMMPS library in use includes the
+This function checks whether the LAMMPS library in use includes the
 specific :doc:`LAMMPS package <Packages>` provided as argument.
 \endverbatim
  *
@@ -5183,8 +5352,8 @@ 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
+You can use, for example, :cpp:func:`lammps_extract_setting` to obtain the
-of local atoms `nlocal` and then assume the dimensions of the returned
+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
@@ -5470,7 +5639,7 @@ void lammps_fix_external_set_vector_length(void *handle, const char *id, int len
 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
+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
@@ -5568,7 +5737,7 @@ int lammps_is_running(void *handle)
   return lmp->update->whichflag;
 }
 
-/** Force a timeout to cleanly stop an ongoing run
+/** Force a timeout to stop an ongoing run cleanly.
  *
  * This function can be used from signal handlers or multi-threaded
  * applications to cleanly terminate an ongoing run.
@@ -5594,9 +5763,9 @@ has thrown a :ref:`C++ exception <exceptions>`.
 .. note::
 
    This function will always report "no error" when the LAMMPS library
-   has been compiled without ``-DLAMMPS_EXCEPTIONS`` which turns fatal
+   has been compiled without ``-DLAMMPS_EXCEPTIONS``, which turns fatal
-   errors aborting LAMMPS into a C++ exceptions. You can use the library
+   errors aborting LAMMPS into C++ exceptions. You can use the library
-   function :cpp:func:`lammps_config_has_exceptions` to check if this is
+   function :cpp:func:`lammps_config_has_exceptions` to check whether this is
    the case.
 \endverbatim
  *
@@ -5626,15 +5795,15 @@ error message is longer, it will be truncated accordingly.  The return
 value of the function corresponds to the kind of error: a "1" indicates
 an error that occurred on all MPI ranks and is often recoverable, while
 a "2" indicates an abort that would happen only in a single MPI rank
-and thus may not be recoverable as other MPI ranks may be waiting on
+and thus may not be recoverable, as other MPI ranks may be waiting on
 the failing MPI ranks to send messages.
 
 .. note::
 
    This function will do nothing when the LAMMPS library has been
-   compiled without ``-DLAMMPS_EXCEPTIONS`` which turns errors aborting
+   compiled without ``-DLAMMPS_EXCEPTIONS``, which turns errors aborting
-   LAMMPS into a C++ exceptions.  You can use the library function
+   LAMMPS into C++ exceptions.  You can use the library function
-   :cpp:func:`lammps_config_has_exceptions` to check if this is the case.
+   :cpp:func:`lammps_config_has_exceptions` to check whether this is the case.
 \endverbatim
  *
  * \param  handle    pointer to a previously created LAMMPS instance cast to ``void *``.

src/library.h

@@ -40,7 +40,8 @@
 
 /** Data type constants for extracting data from atoms, computes and fixes
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in python/lammps/constants.py,
+ * fortran/lammps.f90, and tools/swig/lammps.i */
 
 enum _LMP_DATATYPE_CONST {
   LAMMPS_INT = 0,       /*!< 32-bit integer (array) */
@@ -54,7 +55,8 @@ enum _LMP_DATATYPE_CONST {
 
 /** Style constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in python/lammps/constants.py,
+ * fortran/lammps.f90, and tools/swig/lammps.i */
 
 enum _LMP_STYLE_CONST {
   LMP_STYLE_GLOBAL = 0, /*!< return global data */
@@ -64,7 +66,8 @@ enum _LMP_STYLE_CONST {
 
 /** Type and size constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in python/lammps/constants.py,
+ * fortran/lammps.f90, and tools/swig/lammps.i */
 
 enum _LMP_TYPE_CONST {
   LMP_TYPE_SCALAR = 0, /*!< return scalar */
@@ -77,7 +80,8 @@ enum _LMP_TYPE_CONST {
 
 /** Error codes to select the suitable function in the Error class
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in python/lammps/constants.py,
+ * fortran/lammps.f90, and tools/swig/lammps.i */
 
 enum _LMP_ERROR_CONST {
   LMP_ERROR_WARNING = 0, /*!< call Error::warning() */
@@ -87,6 +91,18 @@ enum _LMP_ERROR_CONST {
   LMP_ERROR_UNIVERSE = 8 /*!< error on Comm::universe */
 };
 
+/** Variable style constants for extracting data from variables.
+ *
+ * Must be kept in sync with the equivalent constants in python/lammps/constants.py,
+ * fortran/lammps.f90, and tools/swig/lammps.i */
+
+enum _LMP_VAR_CONST {
+  LMP_VAR_EQUAL = 0,  /*!< compatible with equal-style variables */
+  LMP_VAR_ATOM = 1,   /*!< compatible with atom-style variables */
+  LMP_VAR_VECTOR = 2, /*!< compatible with vector-style variables */
+  LMP_VAR_STRING = 3  /*!< return value will be a string (catch-all) */
+};
+
 /* Ifdefs to allow this file to be included in C and C++ programs */
 
 #ifdef __cplusplus
@@ -153,6 +169,7 @@ void *lammps_extract_atom(void *handle, const char *name);
 void *lammps_extract_compute(void *handle, const char *, int, int);
 void *lammps_extract_fix(void *handle, const char *, int, int, int, int);
 void *lammps_extract_variable(void *handle, const char *, const char *);
+int lammps_extract_variable_datatype(void *handle, const char *name);
 int lammps_set_variable(void *, char *, char *);
 
 /* ----------------------------------------------------------------------

tools/coding_standard/whitespace.py

@@ -24,11 +24,13 @@ include:
     - cmake/**
     - doc
     - doc/src/**
-    - python
+    - fortran/**
+    - python/**
     - src/**
     - lib/**
     - tools/coding_standard
     - tools/python
+    - unittest/**
 exclude:
     - lib/colvars/Install.py
     - lib/gpu/geryon/file_to_cstr.sh

tools/swig/lammps.i

@@ -22,6 +22,11 @@
 
 %{
 
+/** Data type constants for extracting data from atoms, computes and fixes
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
 enum _LMP_DATATYPE_CONST {
   LAMMPS_INT       = 0,     /*!< 32-bit integer (array) */
   LAMMPS_INT_2D    = 1,     /*!< two-dimensional 32-bit integer array */
@@ -34,7 +39,8 @@ enum _LMP_DATATYPE_CONST {
 
 /** Style constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
 
 enum _LMP_STYLE_CONST {
   LMP_STYLE_GLOBAL = 0,     /*!< return global data */
@@ -44,7 +50,8 @@ enum _LMP_STYLE_CONST {
 
 /** Type and size constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
 
 enum _LMP_TYPE_CONST {
   LMP_TYPE_SCALAR = 0,      /*!< return scalar */
@@ -55,6 +62,31 @@ enum _LMP_TYPE_CONST {
   LMP_SIZE_COLS   = 5       /*!< return number of columns */
 };
 
+/** Error codes to select the suitable function in the Error class
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
+enum _LMP_ERROR_CONST {
+  LMP_ERROR_WARNING  = 0,   /*!< call Error::warning() */
+  LMP_ERROR_ONE      = 1,   /*!< called from one MPI rank */
+  LMP_ERROR_ALL      = 2,   /*!< called from all MPI ranks */
+  LMP_ERROR_WORLD    = 4,   /*!< error on Comm::world */
+  LMP_ERROR_UNIVERSE = 8    /*!< error on Comm::universe */
+};
+
+/** Variable style constants for extracting data from variables.
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
+enum _LMP_VAR_CONST {
+  LMP_VAR_EQUAL  = 0,       /*!< compatible with equal-style variables */
+  LMP_VAR_ATOM   = 1,       /*!< compatible with atom-style variables */
+  LMP_VAR_VECTOR = 2,       /*!< compatible with vector-style variables */
+  LMP_VAR_STRING = 3        /*!< return value will be a string (catch-all) */
+};
+
 /*
  extern void *lammps_open(int argc, char **argv, MPI_Comm comm, void **ptr);
 */
@@ -65,10 +97,13 @@ extern void   lammps_mpi_init();
 extern void   lammps_mpi_finalize();
 extern void   lammps_kokkos_finalize();
 extern void   lammps_python_finalize();
+extern void   lammps_error(void *handle, int error_type, const char *error_text);
+
 extern void   lammps_file(void *handle, const char *file);
 extern char  *lammps_command(void *handle, const char *cmd);
 extern void   lammps_commands_list(void *handle, int ncmd, const char **cmds);
 extern void   lammps_commands_string(void *handle, const char *str);
+
 extern double lammps_get_natoms(void *handle);
 extern double lammps_get_thermo(void *handle, const char *keyword);
 extern void   lammps_extract_box(void *handle, double *boxlo, double *boxhi,
@@ -81,12 +116,16 @@ extern int    lammps_get_mpi_comm(void *handle);
 extern int    lammps_extract_setting(void *handle, const char *keyword);
 extern int    lammps_extract_global_datatype(void *handle, const char *name);
 extern void  *lammps_extract_global(void *handle, const char *name);
+
 extern int    lammps_extract_atom_datatype(void *handle, const char *name);
 extern void  *lammps_extract_atom(void *handle, const char *name);
+
 extern void  *lammps_extract_compute(void *handle, char *id, int, int);
 extern void  *lammps_extract_fix(void *handle, char *, int, int, int, int);
 extern void  *lammps_extract_variable(void *handle, char *, char *);
+extern int    lammps_extract_variable_datatype(void *handle, const char *name);
 extern int    lammps_set_variable(void *, char *, char *);
+
 extern void   lammps_gather_atoms(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_concat(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
@@ -107,6 +146,7 @@ extern int    lammps_find_fix_neighlist(void*, char *, int);
 extern int    lammps_find_compute_neighlist(void*, char *, int);
 extern int    lammps_neighlist_num_elements(void*, int);
 extern void   lammps_neighlist_element_neighbors(void *, int, int, int *, int *, int ** );
+
 extern int    lammps_version(void *handle);
 extern void   lammps_get_os_info(char *buffer, int buf_size);
 extern int    lammps_config_has_mpi_support();
@@ -151,15 +191,21 @@ extern void   lammps_fix_external_set_virial_global(void *handle, const char *id
 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_flush_buffers(void *ptr);
 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);
+extern int    lammps_python_api_version();
 
-extern void   lammps_flush_buffers(void *ptr);
 %}
 
+/** Data type constants for extracting data from atoms, computes and fixes
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
 enum _LMP_DATATYPE_CONST {
   LAMMPS_INT       = 0,     /*!< 32-bit integer (array) */
   LAMMPS_INT_2D    = 1,     /*!< two-dimensional 32-bit integer array */
@@ -172,7 +218,8 @@ enum _LMP_DATATYPE_CONST {
 
 /** Style constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
 
 enum _LMP_STYLE_CONST {
   LMP_STYLE_GLOBAL = 0,     /*!< return global data */
@@ -182,7 +229,8 @@ enum _LMP_STYLE_CONST {
 
 /** Type and size constants for extracting data from computes and fixes.
  *
- * Must be kept in sync with the equivalent constants in lammps/constants.py */
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
 
 enum _LMP_TYPE_CONST {
   LMP_TYPE_SCALAR = 0,      /*!< return scalar */
@@ -193,6 +241,31 @@ enum _LMP_TYPE_CONST {
   LMP_SIZE_COLS   = 5       /*!< return number of columns */
 };
 
+/** Error codes to select the suitable function in the Error class
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
+enum _LMP_ERROR_CONST {
+  LMP_ERROR_WARNING  = 0,   /*!< call Error::warning() */
+  LMP_ERROR_ONE      = 1,   /*!< called from one MPI rank */
+  LMP_ERROR_ALL      = 2,   /*!< called from all MPI ranks */
+  LMP_ERROR_WORLD    = 4,   /*!< error on Comm::world */
+  LMP_ERROR_UNIVERSE = 8    /*!< error on Comm::universe */
+};
+
+/** Variable style constants for extracting data from variables.
+ *
+ * Must be kept in sync with the equivalent constants in src/library.h,
+ * python/lammps/constants.py, and fortran/lammps.f90 */
+
+enum _LMP_VAR_CONST {
+  LMP_VAR_EQUAL  = 0,       /*!< compatible with equal-style variables */
+  LMP_VAR_ATOM   = 1,       /*!< compatible with atom-style variables */
+  LMP_VAR_VECTOR = 2,       /*!< compatible with vector-style variables */
+  LMP_VAR_STRING = 3        /*!< return value will be a string (catch-all) */
+};
+
 /*
  extern void *lammps_open(int argc, char **argv, MPI_Comm comm, void **ptr);
 */
@@ -203,10 +276,13 @@ extern void   lammps_mpi_init();
 extern void   lammps_mpi_finalize();
 extern void   lammps_kokkos_finalize();
 extern void   lammps_python_finalize();
+extern void   lammps_error(void *handle, int error_type, const char *error_text);
+
 extern void   lammps_file(void *handle, const char *file);
 extern char  *lammps_command(void *handle, const char *cmd);
 extern void   lammps_commands_list(void *handle, int ncmd, const char **cmds);
 extern void   lammps_commands_string(void *handle, const char *str);
+
 extern double lammps_get_natoms(void *handle);
 extern double lammps_get_thermo(void *handle, const char *keyword);
 extern void   lammps_extract_box(void *handle, double *boxlo, double *boxhi,
@@ -219,12 +295,16 @@ extern int    lammps_get_mpi_comm(void *handle);
 extern int    lammps_extract_setting(void *handle, const char *keyword);
 extern int    lammps_extract_global_datatype(void *handle, const char *name);
 extern void  *lammps_extract_global(void *handle, const char *name);
+
 extern int    lammps_extract_atom_datatype(void *handle, const char *name);
 extern void  *lammps_extract_atom(void *handle, const char *name);
+
 extern void  *lammps_extract_compute(void *handle, char *id, int, int);
 extern void  *lammps_extract_fix(void *handle, char *, int, int, int, int);
 extern void  *lammps_extract_variable(void *handle, char *, char *);
+extern int    lammps_extract_variable_datatype(void *handle, const char *name);
 extern int    lammps_set_variable(void *, char *, char *);
+
 extern void   lammps_gather_atoms(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_concat(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
@@ -245,6 +325,7 @@ extern int    lammps_find_fix_neighlist(void*, char *, int);
 extern int    lammps_find_compute_neighlist(void*, char *, int);
 extern int    lammps_neighlist_num_elements(void*, int);
 extern void   lammps_neighlist_element_neighbors(void *, int, int, int *, int *, int ** );
+
 extern int    lammps_version(void *handle);
 extern void   lammps_get_os_info(char *buffer, int buf_size);
 extern int    lammps_config_has_mpi_support();
@@ -268,12 +349,18 @@ 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);
 /*
+ * 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);
 */
@@ -283,12 +370,12 @@ extern void   lammps_fix_external_set_virial_global(void *handle, const char *id
 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_flush_buffers(void *ptr);
 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);
+extern int    lammps_python_api_version();
 
-extern void   lammps_flush_buffers(void *ptr);
+/* last revised on 3 October 2022 */
-
-/* last revised on 4 February 2022 */

unittest/c-library/test_library_properties.cpp

@@ -434,6 +434,33 @@ TEST_F(LibraryProperties, neighlist)
     }
 };
 
+TEST_F(LibraryProperties, has_error)
+{
+    // need errors to throw exceptions to be able to intercept them.
+    if (!lammps_config_has_exceptions()) GTEST_SKIP();
+
+    EXPECT_EQ(lammps_has_error(lmp), 0);
+
+    // trigger an error, but hide output
+    ::testing::internal::CaptureStdout();
+    lammps_command(lmp, "this_is_not_a_known_command");
+    ::testing::internal::GetCapturedStdout();
+
+    EXPECT_EQ(lammps_has_error(lmp), 1);
+
+    // retrieve error message
+    char errmsg[1024];
+    int err = lammps_get_last_error_message(lmp, errmsg, 1024);
+    EXPECT_EQ(err, 1);
+    EXPECT_THAT(errmsg, HasSubstr("ERROR: Unknown command: this_is_not_a_known_command"));
+
+    // retrieving the error message clear the error status
+    EXPECT_EQ(lammps_has_error(lmp), 0);
+    err = lammps_get_last_error_message(lmp, errmsg, 1024);
+    EXPECT_EQ(err, 0);
+    EXPECT_THAT(errmsg, StrEq(""));
+};
+
 class AtomProperties : public ::testing::Test {
 protected:
     void *lmp;

unittest/fortran/CMakeLists.txt

@@ -50,13 +50,34 @@ if(CMAKE_Fortran_COMPILER)
   add_test(NAME FortranBox COMMAND test_fortran_box)
 
   add_executable(test_fortran_properties wrap_properties.cpp test_fortran_properties.f90 test_fortran_commands.f90)
-  target_link_libraries(test_fortran_properties PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  target_link_libraries(test_fortran_properties PRIVATE flammps lammps MPI::MPI_Fortran GTest::GMockMain)
   add_test(NAME FortranProperties COMMAND test_fortran_properties)
 
   add_executable(test_fortran_extract_global wrap_extract_global.cpp test_fortran_extract_global.f90)
   target_link_libraries(test_fortran_extract_global PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
   add_test(NAME FortranExtractGlobal COMMAND test_fortran_extract_global)
 
+  add_executable(test_fortran_extract_atom wrap_extract_atom.cpp test_fortran_extract_atom.f90)
+  target_link_libraries(test_fortran_extract_atom PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  add_test(NAME FortranExtractAtom COMMAND test_fortran_extract_atom)
+
+  add_executable(test_fortran_extract_compute wrap_extract_compute.cpp test_fortran_extract_compute.f90)
+  target_link_libraries(test_fortran_extract_compute PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  add_test(NAME FortranExtractCompute COMMAND test_fortran_extract_compute)
+
+  add_executable(test_fortran_extract_fix wrap_extract_fix.cpp test_fortran_extract_fix.f90)
+  target_link_libraries(test_fortran_extract_fix PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  add_test(NAME FortranExtractFix COMMAND test_fortran_extract_fix)
+
+  add_executable(test_fortran_extract_variable wrap_extract_variable.cpp test_fortran_extract_variable.f90)
+  target_compile_definitions(test_fortran_extract_variable PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
+  target_link_libraries(test_fortran_extract_variable PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  add_test(NAME FortranExtractVariable COMMAND test_fortran_extract_variable)
+
+  add_executable(test_fortran_gather_scatter wrap_gather_scatter.cpp test_fortran_gather_scatter.f90)
+  target_link_libraries(test_fortran_gather_scatter PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
+  add_test(NAME FortranGatherScatter COMMAND test_fortran_gather_scatter)
+
 else()
   message(STATUS "Skipping Tests for the LAMMPS Fortran Module: no Fortran compiler")
 endif()

unittest/fortran/atomdata.txt

@@ -0,0 +1,8 @@
+3
+2 1.6
+1 5.2
+3 -1.4
+
+2
+3 2.5
+1 -1.1

unittest/fortran/greetings.txt

@@ -0,0 +1,9 @@
+hello
+god_dag
+hola
+bonjour
+guten_Tag
+konnichiwa
+shalom
+salve
+goedendag

unittest/fortran/keepstuff.f90

@@ -3,24 +3,26 @@ MODULE keepstuff
   IMPLICIT NONE
   TYPE(LAMMPS) :: lmp
   INTEGER :: mycomm
-  CHARACTER(len=40), DIMENSION(3), PARAMETER :: demo_input = &
+  CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: demo_input = &
-      [ CHARACTER(len=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                &
       'region       box block 0 $x 0 2 0 2',                &
       'create_box 1 box',                                   &
       'create_atoms 1 single 1.0 1.0 ${zpos}' ]
   CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: big_input = &
-      [ CHARACTER(len=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                &
       'region       box block 0 $x 0 3 0 4',                &
       'create_box 1 box',                                   &
       'create_atoms 1 single 1.0 1.0 ${zpos}' ]
-  CHARACTER(len=40), DIMENSION(2), PARAMETER :: cont_input = &
+  CHARACTER(LEN=40), DIMENSION(2), PARAMETER :: cont_input = &
-      [ CHARACTER(len=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                &
       'create_atoms 1 single &',                            &
       ' 0.2 0.1 0.1' ]
+  CHARACTER(LEN=40), DIMENSION(1), PARAMETER :: more_input = &
+      [ CHARACTER(LEN=40) :: 'create_atoms 1 single 0.5 0.5 0.5' ]
   CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: pair_input = &
       [ CHARACTER(LEN=40) ::                                &
       'pair_style lj/cut 2.5',                              &
       'pair_coeff 1 1 1.0 1.0',                             &
-      'mass 1 1.0' ]
+      'mass 1 2.0' ]
 END MODULE keepstuff
 

unittest/fortran/test_fortran_extract_atom.f90

@@ -0,0 +1,123 @@
+FUNCTION f_lammps_with_args() BIND(C, name="f_lammps_with_args")
+  USE 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, name="f_lammps_close")
+  USE 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_extract_atom() BIND(C)
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp, big_input, cont_input, pair_input
+   IMPLICIT NONE
+
+   CALL lmp%commands_list(big_input)
+   CALL lmp%commands_list(cont_input)
+   CALL lmp%commands_list(pair_input)
+END SUBROUTINE f_lammps_setup_extract_atom
+
+FUNCTION f_lammps_extract_atom_mass() BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   REAL(c_double) :: f_lammps_extract_atom_mass
+   REAL(c_double), DIMENSION(:), POINTER :: mass => NULL()
+
+   mass = lmp%extract_atom('mass')
+   f_lammps_extract_atom_mass = mass(1)
+END FUNCTION f_lammps_extract_atom_mass
+
+FUNCTION f_lammps_extract_atom_tag_int(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   INTEGER(c_int) :: f_lammps_extract_atom_tag_int
+   INTEGER(c_int), DIMENSION(:), POINTER :: tag => NULL()
+
+   tag = lmp%extract_atom('id')
+   f_lammps_extract_atom_tag_int = tag(i)
+END FUNCTION f_lammps_extract_atom_tag_int
+
+FUNCTION f_lammps_extract_atom_tag_int64(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int64_t
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int64_t), INTENT(IN), VALUE :: i
+   INTEGER(c_int64_t) :: f_lammps_extract_atom_tag_int64
+   INTEGER(c_int64_t), DIMENSION(:), POINTER :: tag => NULL()
+
+   tag = lmp%extract_atom('id')
+   f_lammps_extract_atom_tag_int64 = tag(i)
+END FUNCTION f_lammps_extract_atom_tag_int64
+
+FUNCTION f_lammps_extract_atom_type(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   INTEGER(c_int) :: f_lammps_extract_atom_type
+   INTEGER(c_int), DIMENSION(:), POINTER :: atype => NULL()
+
+   atype = lmp%extract_atom('type')
+   f_lammps_extract_atom_type = atype(i)
+END FUNCTION f_lammps_extract_atom_type
+
+FUNCTION f_lammps_extract_atom_mask(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   INTEGER(c_int) :: f_lammps_extract_atom_mask
+   INTEGER(c_int), DIMENSION(:), POINTER :: mask => NULL()
+
+   mask = lmp%extract_atom('mask')
+   f_lammps_extract_atom_mask = mask(i)
+END FUNCTION f_lammps_extract_atom_mask
+
+SUBROUTINE f_lammps_extract_atom_x(i, x) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double), DIMENSION(3) :: x
+   REAL(c_double), DIMENSION(:,:), POINTER :: xptr => NULL()
+
+   xptr = lmp%extract_atom('x')
+   x = xptr(:,i)
+END SUBROUTINE f_lammps_extract_atom_x
+
+SUBROUTINE f_lammps_extract_atom_v(i, v) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double), DIMENSION(3) :: v
+   REAL(c_double), DIMENSION(:,:), POINTER :: vptr => NULL()
+
+   vptr = lmp%extract_atom('v')
+   v = vptr(:,i)
+END SUBROUTINE f_lammps_extract_atom_v

unittest/fortran/test_fortran_extract_compute.f90

@@ -0,0 +1,133 @@
+FUNCTION f_lammps_with_args() BIND(C)
+  USE 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 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_extract_compute() BIND(C)
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp, big_input, cont_input, more_input, pair_input
+   IMPLICIT NONE
+
+   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("compute peratompe all pe/atom") ! per-atom vector
+   call lmp%command("compute stress all stress/atom thermo_temp") ! per-atom array
+   CALL lmp%command("compute totalpe all reduce sum c_peratompe") ! global scalar
+   CALL lmp%command("compute COM all com") ! global vector
+   CALL lmp%command("compute RDF all rdf 100") ! global array
+   CALL lmp%command("compute pairdist all pair/local dist") ! local vector
+   CALL lmp%command("compute pairlocal all pair/local dist dx dy dz") ! local array
+   CALL lmp%command("thermo_style custom step pe c_totalpe c_COM[1]")
+   CALL lmp%command("run 0") ! must be here, otherwise will SEGFAULT
+END SUBROUTINE f_lammps_setup_extract_compute
+
+FUNCTION f_lammps_extract_compute_peratom_vector(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double) :: f_lammps_extract_compute_peratom_vector
+   REAL(c_double), DIMENSION(:), POINTER :: vector => NULL()
+
+   vector = lmp%extract_compute('peratompe', lmp%style%atom, lmp%type%vector)
+   f_lammps_extract_compute_peratom_vector = vector(i)
+END FUNCTION f_lammps_extract_compute_peratom_vector
+
+FUNCTION f_lammps_extract_compute_peratom_array(i,j) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i, j
+   REAL(c_double) :: f_lammps_extract_compute_peratom_array
+   REAL(c_double), DIMENSION(:,:), POINTER :: array => NULL()
+
+   array = lmp%extract_compute('stress', lmp%style%atom, lmp%type%array)
+   f_lammps_extract_compute_peratom_array = array(i,j)
+END FUNCTION f_lammps_extract_compute_peratom_array
+
+FUNCTION f_lammps_extract_compute_global_scalar() BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   REAL(c_double) :: f_lammps_extract_compute_global_scalar
+   REAL(c_double), POINTER :: scalar
+
+   scalar = lmp%extract_compute('totalpe', lmp%style%global, lmp%type%scalar)
+   f_lammps_extract_compute_global_scalar = scalar
+END FUNCTION f_lammps_extract_compute_global_scalar
+
+FUNCTION f_lammps_extract_compute_global_vector(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double) :: f_lammps_extract_compute_global_vector
+   REAL(c_double), DIMENSION(:), POINTER :: vector
+
+   vector = lmp%extract_compute('COM', lmp%style%global, lmp%type%vector)
+   f_lammps_extract_compute_global_vector = vector(i)
+END FUNCTION f_lammps_extract_compute_global_vector
+
+FUNCTION f_lammps_extract_compute_global_array(i,j) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i, j
+   REAL(c_double) :: f_lammps_extract_compute_global_array
+   REAL(c_double), DIMENSION(:,:), POINTER :: array
+
+   array = lmp%extract_compute('RDF', lmp%style%global, lmp%type%array)
+   f_lammps_extract_compute_global_array = array(i,j)
+END FUNCTION f_lammps_extract_compute_global_array
+
+FUNCTION f_lammps_extract_compute_local_vector(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double) :: f_lammps_extract_compute_local_vector
+   REAL(c_double), DIMENSION(:), POINTER :: vector
+
+   vector = lmp%extract_compute('pairdist', lmp%style%local, lmp%type%vector)
+   f_lammps_extract_compute_local_vector = vector(i)
+END FUNCTION f_lammps_extract_compute_local_vector
+
+FUNCTION f_lammps_extract_compute_local_array(i, j) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i, j
+   REAL(c_double) :: f_lammps_extract_compute_local_array
+   REAL(c_double), DIMENSION(:,:), POINTER :: array
+
+   array = lmp%extract_compute('pairlocal', lmp%style%local, lmp%type%array)
+   f_lammps_extract_compute_local_array = array(i,j)
+END FUNCTION f_lammps_extract_compute_local_array

unittest/fortran/test_fortran_extract_fix.f90

@@ -0,0 +1,99 @@
+FUNCTION f_lammps_with_args() BIND(C)
+  USE 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 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_extract_fix() BIND(C)
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp, big_input, cont_input, pair_input, more_input
+   IMPLICIT NONE
+
+   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("fix state all store/state 0 z") ! per-atom vector
+   CALL lmp%command("fix move all move linear 0 0 0") ! for per-atom array
+   CALL lmp%command("fix recenter all recenter 1 1 1") ! global scalar, vector
+   CALL lmp%command("variable natoms equal count(all)")
+   CALL lmp%command("variable ts equal step")
+   CALL lmp%command("fix vec all vector 1 v_natoms v_ts") ! global array
+   CALL lmp%command("run 1") ! must be 1, otherwise move/recenter won't happen
+END SUBROUTINE f_lammps_setup_extract_fix
+
+FUNCTION f_lammps_extract_fix_global_scalar() BIND(C) RESULT(scalar)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   REAL(c_double) :: scalar
+
+   scalar = lmp%extract_fix("recenter", lmp%style%global, lmp%type%scalar)
+END FUNCTION f_lammps_extract_fix_global_scalar
+
+FUNCTION f_lammps_extract_fix_global_vector(i) BIND(C) RESULT(element)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double) :: element
+
+   element = lmp%extract_fix("recenter", lmp%style%global, lmp%type%vector, i)
+END FUNCTION f_lammps_extract_fix_global_vector
+
+FUNCTION f_lammps_extract_fix_global_array(i,j) BIND(C) RESULT(element)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i, j
+   REAL(c_double) :: element
+
+   element = lmp%extract_fix("vec", lmp%style%global, lmp%type%array, i, j)
+END FUNCTION f_lammps_extract_fix_global_array
+
+FUNCTION f_lammps_extract_fix_peratom_vector(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   REAL(c_double) :: f_lammps_extract_fix_peratom_vector
+   REAL(c_double), DIMENSION(:), POINTER :: vector
+
+   vector = lmp%extract_fix("state", lmp%style%atom, lmp%type%vector, -1, -1)
+   f_lammps_extract_fix_peratom_vector = vector(i)
+END FUNCTION f_lammps_extract_fix_peratom_vector
+
+FUNCTION f_lammps_extract_fix_peratom_array(i,j) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i, j
+   REAL(c_double) :: f_lammps_extract_fix_peratom_array
+   REAL(c_double), DIMENSION(:,:), POINTER :: array
+
+   array = lmp%extract_fix("move", lmp%style%atom, lmp%type%array, -1, -1)
+   f_lammps_extract_fix_peratom_array = array(i,j)
+END FUNCTION f_lammps_extract_fix_peratom_array

unittest/fortran/test_fortran_extract_global.f90

@@ -51,324 +51,330 @@ SUBROUTINE f_lammps_setup_full_extract_global () BIND(C)
 END SUBROUTINE f_lammps_setup_full_extract_global
 
 FUNCTION f_lammps_extract_global_units() BIND(C) RESULT(success)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE LIBLAMMPS
    USE keepstuff, ONLY : lmp
    IMPLICIT NONE
-   INTEGER (C_int) :: success
+   INTEGER(c_int) :: success
    CHARACTER(LEN=16) :: units
 
    ! passing strings from Fortran to C is icky, so we do the test here and
    ! report the result instead
    units = lmp%extract_global('units')
    IF (TRIM(units) == 'lj') THEN
-      success = 1_C_int
+      success = 1_c_int
    ELSE
-      success = 0_C_int
+      success = 0_c_int
    END IF
 END FUNCTION f_lammps_extract_global_units
 
 FUNCTION f_lammps_extract_global_ntimestep() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: ntimestep
+   INTEGER(c_int), POINTER :: ntimestep
-   INTEGER (C_int) :: f_lammps_extract_global_ntimestep
+   INTEGER(c_int) :: f_lammps_extract_global_ntimestep
 
    ntimestep = lmp%extract_global("ntimestep")
    f_lammps_extract_global_ntimestep = ntimestep
 END FUNCTION f_lammps_extract_global_ntimestep
+
 FUNCTION f_lammps_extract_global_ntimestep_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: ntimestep
+   INTEGER(c_int64_t), POINTER :: ntimestep
-   INTEGER (C_int64_t) :: f_lammps_extract_global_ntimestep_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_ntimestep_big
 
    ntimestep = lmp%extract_global("ntimestep")
    f_lammps_extract_global_ntimestep_big = ntimestep
 END FUNCTION f_lammps_extract_global_ntimestep_big
 
 FUNCTION f_lammps_extract_global_dt() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double), POINTER :: dt
+   REAL(c_double), POINTER :: dt
-   REAL (C_double) :: f_lammps_extract_global_dt
+   REAL(c_double) :: f_lammps_extract_global_dt
 
    dt = lmp%extract_global("dt")
    f_lammps_extract_global_dt = dt
 END FUNCTION f_lammps_extract_global_dt
 
 SUBROUTINE f_lammps_extract_global_boxlo(C_boxlo) BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double), DIMENSION(3) :: C_boxlo
+   REAL(c_double), DIMENSION(3) :: C_boxlo
-   REAL (C_double), DIMENSION(:), POINTER :: boxlo
+   REAL(c_double), DIMENSION(:), POINTER :: boxlo
 
    boxlo = lmp%extract_global("boxlo")
    C_boxlo = boxlo
 END SUBROUTINE f_lammps_extract_global_boxlo
 
 SUBROUTINE f_lammps_extract_global_boxhi(C_boxhi) BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double), DIMENSION(3) :: C_boxhi
+   REAL(c_double), DIMENSION(3) :: C_boxhi
-   REAL (C_double), DIMENSION(:), POINTER :: boxhi
+   REAL(c_double), DIMENSION(:), POINTER :: boxhi
 
    boxhi = lmp%extract_global("boxhi")
    C_boxhi = boxhi
 END SUBROUTINE f_lammps_extract_global_boxhi
 
 FUNCTION f_lammps_extract_global_boxxlo() BIND(C) RESULT(C_boxxlo)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxxlo
+   REAL(c_double) :: C_boxxlo
-   REAL (C_double), POINTER :: boxxlo
+   REAL(c_double), POINTER :: boxxlo
 
    boxxlo = lmp%extract_global("boxxlo")
    C_boxxlo = boxxlo
 END FUNCTION f_lammps_extract_global_boxxlo
 
 FUNCTION f_lammps_extract_global_boxxhi() BIND(C) RESULT(C_boxxhi)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxxhi
+   REAL(c_double) :: C_boxxhi
-   REAL (C_double), POINTER :: boxxhi
+   REAL(c_double), POINTER :: boxxhi
 
    boxxhi = lmp%extract_global("boxxhi")
    C_boxxhi = boxxhi
 END FUNCTION f_lammps_extract_global_boxxhi
 
 FUNCTION f_lammps_extract_global_boxylo() BIND(C) RESULT(C_boxylo)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxylo
+   REAL(c_double) :: C_boxylo
-   REAL (C_double), POINTER :: boxylo
+   REAL(c_double), POINTER :: boxylo
 
    boxylo = lmp%extract_global("boxylo")
    C_boxylo = boxylo
 END FUNCTION f_lammps_extract_global_boxylo
 
 FUNCTION f_lammps_extract_global_boxyhi() BIND(C) RESULT(C_boxyhi)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxyhi
+   REAL(c_double) :: C_boxyhi
-   REAL (C_double), POINTER :: boxyhi
+   REAL(c_double), POINTER :: boxyhi
 
    boxyhi = lmp%extract_global("boxyhi")
    C_boxyhi = boxyhi
 END FUNCTION f_lammps_extract_global_boxyhi
 
 FUNCTION f_lammps_extract_global_boxzlo() BIND(C) RESULT(C_boxzlo)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxzlo
+   REAL(c_double) :: C_boxzlo
-   REAL (C_double), POINTER :: boxzlo
+   REAL(c_double), POINTER :: boxzlo
 
    boxzlo = lmp%extract_global("boxzlo")
    C_boxzlo = boxzlo
 END FUNCTION f_lammps_extract_global_boxzlo
 
 FUNCTION f_lammps_extract_global_boxzhi() BIND(C) RESULT(C_boxzhi)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_boxzhi
+   REAL(c_double) :: C_boxzhi
-   REAL (C_double), POINTER :: boxzhi
+   REAL(c_double), POINTER :: boxzhi
 
    boxzhi = lmp%extract_global("boxzhi")
    C_boxzhi = boxzhi
 END FUNCTION f_lammps_extract_global_boxzhi
 
 SUBROUTINE f_lammps_extract_global_periodicity(C_periodicity) BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), DIMENSION(3) :: C_periodicity
+   INTEGER(c_int), DIMENSION(3) :: C_periodicity
-   INTEGER (C_int), DIMENSION(:), POINTER :: periodicity
+   INTEGER(c_int), DIMENSION(:), POINTER :: periodicity
 
    periodicity = lmp%extract_global("periodicity")
    C_periodicity = periodicity
 END SUBROUTINE f_lammps_extract_global_periodicity
 
 FUNCTION f_lammps_extract_global_triclinic() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: triclinic
+   INTEGER(c_int), POINTER :: triclinic
-   INTEGER (C_int) :: f_lammps_extract_global_triclinic
+   INTEGER(c_int) :: f_lammps_extract_global_triclinic
 
    triclinic = lmp%extract_global("triclinic")
    f_lammps_extract_global_triclinic = triclinic
 END FUNCTION f_lammps_extract_global_triclinic
 
 FUNCTION f_lammps_extract_global_xy() BIND(C) RESULT(C_xy)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_xy
+   REAL(c_double) :: C_xy
-   REAL (C_double), POINTER :: xy
+   REAL(c_double), POINTER :: xy
 
    xy = lmp%extract_global("xy")
    C_xy = xy
 END FUNCTION f_lammps_extract_global_xy
 
 FUNCTION f_lammps_extract_global_xz() BIND(C) RESULT(C_xz)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_xz
+   REAL(c_double) :: C_xz
-   REAL (C_double), POINTER :: xz
+   REAL(c_double), POINTER :: xz
 
    xz = lmp%extract_global("xz")
    C_xz = xz
 END FUNCTION f_lammps_extract_global_xz
 
 FUNCTION f_lammps_extract_global_yz() BIND(C) RESULT(C_yz)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_yz
+   REAL(c_double) :: C_yz
-   REAL (C_double), POINTER :: yz
+   REAL(c_double), POINTER :: yz
 
    yz = lmp%extract_global("yz")
    C_yz = yz
 END FUNCTION f_lammps_extract_global_yz
 
 FUNCTION f_lammps_extract_global_natoms() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: natoms
+   INTEGER(c_int), POINTER :: natoms
-   INTEGER (C_int) :: f_lammps_extract_global_natoms
+   INTEGER(c_int) :: f_lammps_extract_global_natoms
 
    natoms = lmp%extract_global("natoms")
    f_lammps_extract_global_natoms = natoms
 END FUNCTION f_lammps_extract_global_natoms
+
 FUNCTION f_lammps_extract_global_natoms_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: natoms
+   INTEGER(c_int64_t), POINTER :: natoms
-   INTEGER (C_int64_t) :: f_lammps_extract_global_natoms_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_natoms_big
 
    natoms = lmp%extract_global("natoms")
    f_lammps_extract_global_natoms_big = natoms
 END FUNCTION f_lammps_extract_global_natoms_big
 
 FUNCTION f_lammps_extract_global_nbonds() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nbonds
+   INTEGER(c_int), POINTER :: nbonds
-   INTEGER (C_int) :: f_lammps_extract_global_nbonds
+   INTEGER(c_int) :: f_lammps_extract_global_nbonds
 
    nbonds = lmp%extract_global("nbonds")
    f_lammps_extract_global_nbonds = nbonds
 END FUNCTION f_lammps_extract_global_nbonds
+
 FUNCTION f_lammps_extract_global_nbonds_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: nbonds
+   INTEGER(c_int64_t), POINTER :: nbonds
-   INTEGER (C_int64_t) :: f_lammps_extract_global_nbonds_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_nbonds_big
 
    nbonds = lmp%extract_global("nbonds")
    f_lammps_extract_global_nbonds_big = nbonds
 END FUNCTION f_lammps_extract_global_nbonds_big
 
 FUNCTION f_lammps_extract_global_nangles() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nangles
+   INTEGER(c_int), POINTER :: nangles
-   INTEGER (C_int) :: f_lammps_extract_global_nangles
+   INTEGER(c_int) :: f_lammps_extract_global_nangles
 
    nangles = lmp%extract_global("nangles")
    f_lammps_extract_global_nangles = nangles
 END FUNCTION f_lammps_extract_global_nangles
+
 FUNCTION f_lammps_extract_global_nangles_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: nangles
+   INTEGER(c_int64_t), POINTER :: nangles
-   INTEGER (C_int64_t) :: f_lammps_extract_global_nangles_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_nangles_big
 
    nangles = lmp%extract_global("nangles")
    f_lammps_extract_global_nangles_big = nangles
 END FUNCTION f_lammps_extract_global_nangles_big
 
 FUNCTION f_lammps_extract_global_ndihedrals() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: ndihedrals
+   INTEGER(c_int), POINTER :: ndihedrals
-   INTEGER (C_int) :: f_lammps_extract_global_ndihedrals
+   INTEGER(c_int) :: f_lammps_extract_global_ndihedrals
 
    ndihedrals = lmp%extract_global("ndihedrals")
    f_lammps_extract_global_ndihedrals = ndihedrals
 END FUNCTION f_lammps_extract_global_ndihedrals
+
 FUNCTION f_lammps_extract_global_ndihedrals_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: ndihedrals
+   INTEGER(c_int64_t), POINTER :: ndihedrals
-   INTEGER (C_int64_t) :: f_lammps_extract_global_ndihedrals_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_ndihedrals_big
 
    ndihedrals = lmp%extract_global("ndihedrals")
    f_lammps_extract_global_ndihedrals_big = ndihedrals
 END FUNCTION f_lammps_extract_global_ndihedrals_big
 
 FUNCTION f_lammps_extract_global_nimpropers() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nimpropers
+   INTEGER(c_int), POINTER :: nimpropers
-   INTEGER (C_int) :: f_lammps_extract_global_nimpropers
+   INTEGER(c_int) :: f_lammps_extract_global_nimpropers
 
    nimpropers = lmp%extract_global("nimpropers")
    f_lammps_extract_global_nimpropers = nimpropers
 END FUNCTION f_lammps_extract_global_nimpropers
+
 FUNCTION f_lammps_extract_global_nimpropers_big() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int64_t), POINTER :: nimpropers
+   INTEGER(c_int64_t), POINTER :: nimpropers
-   INTEGER (C_int64_t) :: f_lammps_extract_global_nimpropers_big
+   INTEGER(c_int64_t) :: f_lammps_extract_global_nimpropers_big
 
    nimpropers = lmp%extract_global("nimpropers")
    f_lammps_extract_global_nimpropers_big = nimpropers
@@ -376,96 +382,96 @@ END FUNCTION f_lammps_extract_global_nimpropers_big
 
 
 FUNCTION f_lammps_extract_global_ntypes() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: ntypes
+   INTEGER(c_int), POINTER :: ntypes
-   INTEGER (C_int) :: f_lammps_extract_global_ntypes
+   INTEGER(c_int) :: f_lammps_extract_global_ntypes
 
    ntypes = lmp%extract_global("ntypes")
    f_lammps_extract_global_ntypes = ntypes
 END FUNCTION f_lammps_extract_global_ntypes
 
 FUNCTION f_lammps_extract_global_nlocal() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nlocal
+   INTEGER(c_int), POINTER :: nlocal
-   INTEGER (C_int) :: f_lammps_extract_global_nlocal
+   INTEGER(c_int) :: f_lammps_extract_global_nlocal
 
    nlocal = lmp%extract_global("nlocal")
    f_lammps_extract_global_nlocal = nlocal
 END FUNCTION f_lammps_extract_global_nlocal
 
 FUNCTION f_lammps_extract_global_nghost() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nghost
+   INTEGER(c_int), POINTER :: nghost
-   INTEGER (C_int) :: f_lammps_extract_global_nghost
+   INTEGER(c_int) :: f_lammps_extract_global_nghost
 
    nghost = lmp%extract_global("nghost")
    f_lammps_extract_global_nghost = nghost
 END FUNCTION f_lammps_extract_global_nghost
 
 FUNCTION f_lammps_extract_global_nmax() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int), POINTER :: nmax
+   INTEGER(c_int), POINTER :: nmax
-   INTEGER (C_int) :: f_lammps_extract_global_nmax
+   INTEGER(c_int) :: f_lammps_extract_global_nmax
 
    nmax = lmp%extract_global("nmax")
    f_lammps_extract_global_nmax = nmax
 END FUNCTION f_lammps_extract_global_nmax
 
 FUNCTION f_lammps_extract_global_boltz() BIND(C) RESULT(C_k_B)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_k_B
+   REAL(c_double) :: C_k_B
-   REAL (C_double), POINTER :: k_B
+   REAL(c_double), POINTER :: k_B
 
    k_B = lmp%extract_global("boltz")
    C_k_B = k_B
 END FUNCTION f_lammps_extract_global_boltz
 
 FUNCTION f_lammps_extract_global_hplanck() BIND(C) RESULT(C_h)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: C_h
+   REAL(c_double) :: C_h
-   REAL (C_double), POINTER :: h
+   REAL(c_double), POINTER :: h
 
    h = lmp%extract_global("boltz")
    C_h = h
 END FUNCTION f_lammps_extract_global_hplanck
 
 FUNCTION f_lammps_extract_global_angstrom() BIND(C) RESULT(Angstrom)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: Angstrom
+   REAL(c_double) :: Angstrom
-   REAL (C_double), POINTER :: A
+   REAL(c_double), POINTER :: A
 
    A = lmp%extract_global("angstrom")
    Angstrom = A
 END FUNCTION f_lammps_extract_global_angstrom
 
 FUNCTION f_lammps_extract_global_femtosecond() BIND(C) RESULT(fs)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   REAL (C_double) :: fs
+   REAL(c_double) :: fs
-   REAL (C_double), POINTER :: femtosecond
+   REAL(c_double), POINTER :: femtosecond
 
    femtosecond = lmp%extract_global("femtosecond")
    fs = femtosecond

unittest/fortran/test_fortran_extract_variable.f90

@@ -0,0 +1,426 @@
+MODULE keepvar
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr, c_size_t, c_char
+  USE liblammps
+  IMPLICIT NONE
+
+  INTERFACE
+    FUNCTION c_path_join(a, b) BIND(C)
+      IMPORT :: c_ptr
+      TYPE(c_ptr), VALUE :: a, b
+      TYPE(c_ptr) :: c_path_join
+    END FUNCTION c_path_join
+
+    FUNCTION c_strlen(str) BIND(C,name='strlen')
+      IMPORT :: c_ptr, c_size_t
+      IMPLICIT NONE
+      TYPE(c_ptr), INTENT(IN), VALUE :: str
+      INTEGER(c_size_t) :: c_strlen
+    END FUNCTION c_strlen
+
+    SUBROUTINE c_free(ptr) BIND(C,name='free')
+      IMPORT :: c_ptr
+      TYPE(c_ptr), VALUE :: ptr
+    END SUBROUTINE c_free
+  END INTERFACE
+
+CONTAINS
+
+  FUNCTION absolute_path(filename)
+    USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr, c_size_t, c_char, C_F_POINTER
+    USE keepstuff, ONLY : lmp
+    CHARACTER(LEN=:), ALLOCATABLE :: absolute_path
+    CHARACTER(LEN=*), INTENT(IN) :: filename
+    CHARACTER(LEN=256) :: test_input_directory
+    TYPE(c_ptr) :: c_test_input_directory, c_absolute_path, c_filename
+    CHARACTER(LEN=1,KIND=c_char), DIMENSION(:), POINTER :: F_absolute_path
+    INTEGER(c_size_t) :: i, length
+
+    test_input_directory = lmp%extract_variable('input_dir')
+    c_test_input_directory = f2c_string(test_input_directory)
+    c_filename = f2c_string(filename)
+    c_absolute_path = c_path_join(c_test_input_directory, c_filename)
+    length = c_strlen(c_absolute_path)
+    CALL C_F_POINTER(c_absolute_path, F_absolute_path, [length])
+    ALLOCATE(CHARACTER(LEN=length) :: absolute_path)
+    DO i = 1, length
+      absolute_path(i:i) = F_absolute_path(i)
+    END DO
+    CALL c_free(c_filename)
+    CALL c_free(c_test_input_directory)
+    CALL c_free(c_absolute_path)
+  END FUNCTION absolute_path
+
+  FUNCTION f2c_string(f_string) RESULT(ptr)
+    USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_char, c_ptr, c_size_t, &
+      c_null_char, C_F_POINTER
+    CHARACTER(LEN=*), INTENT(IN)           :: f_string
+    CHARACTER(LEN=1, KIND=c_char), POINTER :: c_string(:)
+    TYPE(c_ptr) :: ptr
+    INTEGER(c_size_t) :: i, n
+
+    INTERFACE
+      FUNCTION lammps_malloc(size) BIND(C, name='malloc')
+        IMPORT :: c_ptr, c_size_t
+        IMPLICIT NONE
+        INTEGER(c_size_t), VALUE :: size
+        TYPE(c_ptr) :: lammps_malloc
+      END FUNCTION lammps_malloc
+    END INTERFACE
+
+    n = LEN_TRIM(f_string)
+    ptr = lammps_malloc(n+1)
+    CALL C_F_POINTER(ptr, c_string, [1])
+    DO i=1, n
+        c_string(i) = f_string(i:i)
+    END DO
+    c_string(n+1) = c_null_char
+  END FUNCTION f2c_string
+
+END MODULE keepvar
+
+FUNCTION f_lammps_with_C_args(argc, argv) BIND(C)
+  USE ISO_C_BINDING, ONLY: c_ptr, c_char, c_int, c_size_t, C_F_POINTER
+  USE liblammps
+  USE keepstuff, ONLY: lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: argc
+  TYPE(c_ptr), VALUE :: argv
+  TYPE(c_ptr), DIMENSION(:), POINTER :: Fargv
+  INTEGER, PARAMETER :: ARG_LENGTH = 256
+  TYPE(c_ptr) :: f_lammps_with_C_args
+  CHARACTER(LEN=ARG_LENGTH), DIMENSION(argc) :: args
+  CHARACTER(LEN=1,KIND=c_char), DIMENSION(:), POINTER :: Cstr
+  INTEGER(c_size_t):: i, length, j
+
+  INTERFACE
+    FUNCTION c_strlen(str) BIND(C,name='strlen')
+      IMPORT :: c_ptr, c_size_t
+      IMPLICIT NONE
+      TYPE(c_ptr), INTENT(IN), VALUE :: str
+      INTEGER(c_size_t) :: c_strlen
+    END FUNCTION c_strlen
+  END INTERFACE
+
+  CALL C_F_POINTER(argv, Fargv, [argc])
+  DO i = 1, argc
+    args(i) = ''
+    length = c_strlen(Fargv(i))
+    CALL C_F_POINTER(Fargv(i), Cstr, [length])
+    DO j = 1, length
+      args(i)(j:j) = Cstr(j)
+    END DO
+  END DO
+
+  lmp = lammps(args)
+  f_lammps_with_C_args = lmp%handle
+END FUNCTION f_lammps_with_C_args
+
+SUBROUTINE f_lammps_close() BIND(C)
+  USE 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_extract_variable() BIND(C)
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp, big_input, cont_input, more_input, pair_input
+  USE keepvar, ONLY : absolute_path
+  IMPLICIT NONE
+
+  ! Had to do this one as one string because lammps_commands_list and
+  ! lammps_commands_string do not play well with triple quotes
+  CHARACTER(LEN=256), PARAMETER :: py_input = &
+      'python square_it input 1 v_lp return v_py format ff here """' &
+        // NEW_LINE(' ') // 'def square_it(N) :' &
+        // NEW_LINE(' ') // '  return N*N' &
+        // NEW_LINE(' ') // '"""'
+
+  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('variable idx index "hello" "goodbye"')
+  CALL lmp%command('variable lp loop 10')
+  CALL lmp%command('variable lp_pad loop 10 pad')
+  CALL lmp%command('variable wld world "group1"')
+  CALL lmp%command('variable uni universe "universe1" "universeA"')
+  CALL lmp%command('variable ulp uloop 2')
+  CALL lmp%command('variable str string "this is a string"')
+  CALL lmp%command('variable ex equal exp(v_lp)')
+  CALL lmp%command('variable fmt format ex %.6G')
+  CALL lmp%command('variable fmt_pad format ex %08.6g')
+  CALL lmp%command('variable username getenv FORTRAN_USER')
+  CALL lmp%command('variable greeting file ' // absolute_path('greetings.txt'))
+  CALL lmp%command('variable atfile atomfile ' &
+    // absolute_path('atomdata.txt'))
+  IF (lmp%config_has_package('PYTHON')) THEN
+    CALL lmp%command(py_input)
+    CALL lmp%command('variable py python square_it')
+  END IF
+  CALL lmp%command('variable time timer')
+  CALL lmp%command('variable int internal 4')
+  CALL lmp%command('variable at_z atom z')
+  CALL lmp%command("compute COM all com") ! defines a global vector
+  CALL lmp%command("variable center vector c_COM")
+  ! make sure COM is computable...
+  CALL lmp%command("thermo_style custom step pe c_COM[1] v_center[1]")
+  CALL lmp%command("run 0") ! so c_COM and v_center have values
+END SUBROUTINE f_lammps_setup_extract_variable
+
+FUNCTION f_lammps_extract_variable_index_1() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int) :: f_lammps_extract_variable_index_1
+  CHARACTER(LEN=256) :: str
+
+  str = lmp%extract_variable("idx")
+  IF (trim(str) == 'hello') THEN
+     f_lammps_extract_variable_index_1 = 1_c_int
+  ELSE
+     f_lammps_extract_variable_index_1 = 0_c_int
+  END IF
+END FUNCTION f_lammps_extract_variable_index_1
+
+FUNCTION f_lammps_extract_variable_index_2() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int) :: f_lammps_extract_variable_index_2
+  CHARACTER(LEN=256) :: str
+
+  str = lmp%extract_variable("idx")
+  IF (trim(str) == 'goodbye') THEN
+     f_lammps_extract_variable_index_2 = 1_c_int
+  ELSE
+     f_lammps_extract_variable_index_2 = 0_c_int
+  END IF
+END FUNCTION f_lammps_extract_variable_index_2
+
+FUNCTION f_lammps_extract_variable_loop() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int) :: f_lammps_extract_variable_loop
+  CHARACTER(LEN=256) :: loop
+
+  loop = lmp%extract_variable('lp')
+  READ(loop,*) f_lammps_extract_variable_loop
+END FUNCTION f_lammps_extract_variable_loop
+
+FUNCTION f_lammps_extract_variable_loop_pad() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_loop_pad
+  CHARACTER(LEN=20) :: loop
+
+  loop = lmp%extract_variable('lp_pad')
+  f_lammps_extract_variable_loop_pad = f2c_string(loop)
+END FUNCTION f_lammps_extract_variable_loop_pad
+
+FUNCTION f_lammps_extract_variable_world() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_world
+  CHARACTER(LEN=20) :: world
+
+  world = lmp%extract_variable('wld')
+  f_lammps_extract_variable_world = f2c_string(world)
+END FUNCTION f_lammps_extract_variable_world
+
+FUNCTION f_lammps_extract_variable_universe() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_universe
+  CHARACTER(LEN=20) :: universe
+
+  universe = lmp%extract_variable('uni')
+  f_lammps_extract_variable_universe = f2c_string(universe)
+END FUNCTION f_lammps_extract_variable_universe
+
+FUNCTION f_lammps_extract_variable_uloop() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int) :: f_lammps_extract_variable_uloop
+  CHARACTER(LEN=256) :: uloop
+
+  uloop = lmp%extract_variable('ulp')
+  READ(uloop,*) f_lammps_extract_variable_uloop
+END FUNCTION f_lammps_extract_variable_uloop
+
+FUNCTION f_lammps_extract_variable_string() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_string
+  CHARACTER(LEN=256) :: string
+
+  string = lmp%extract_variable('str')
+  f_lammps_extract_variable_string = f2c_string(string)
+END FUNCTION f_lammps_extract_variable_string
+
+FUNCTION f_lammps_extract_variable_format() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_format
+  CHARACTER(LEN=20) :: form
+
+  form = lmp%extract_variable('fmt')
+  f_lammps_extract_variable_format = f2c_string(form)
+END FUNCTION f_lammps_extract_variable_format
+
+FUNCTION f_lammps_extract_variable_format_pad() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_format_pad
+  CHARACTER(LEN=20) :: form
+
+  form = lmp%extract_variable('fmt_pad')
+  f_lammps_extract_variable_format_pad = f2c_string(form)
+END FUNCTION f_lammps_extract_variable_format_pad
+
+FUNCTION f_lammps_extract_variable_getenv() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_getenv
+  CHARACTER(LEN=40) :: string
+
+  string = lmp%extract_variable('username')
+  f_lammps_extract_variable_getenv = f2c_string(string)
+END FUNCTION f_lammps_extract_variable_getenv
+
+FUNCTION f_lammps_extract_variable_file() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_ptr
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  TYPE(c_ptr) :: f_lammps_extract_variable_file
+  CHARACTER(LEN=40) :: string
+
+  string = lmp%extract_variable('greeting')
+  f_lammps_extract_variable_file = f2c_string(string)
+END FUNCTION f_lammps_extract_variable_file
+
+FUNCTION f_lammps_extract_variable_atomfile(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_extract_variable_atomfile
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: atom_data
+
+  atom_data = lmp%extract_variable('atfile')
+  f_lammps_extract_variable_atomfile = atom_data(i)
+END FUNCTION f_lammps_extract_variable_atomfile
+
+FUNCTION f_lammps_extract_variable_python() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double) :: f_lammps_extract_variable_python
+
+  f_lammps_extract_variable_python = lmp%extract_variable('py')
+END FUNCTION f_lammps_extract_variable_python
+
+FUNCTION f_lammps_extract_variable_timer() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double) :: f_lammps_extract_variable_timer
+
+  f_lammps_extract_variable_timer = lmp%extract_variable('time')
+END FUNCTION f_lammps_extract_variable_timer
+
+FUNCTION f_lammps_extract_variable_internal() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double) :: f_lammps_extract_variable_internal
+
+  f_lammps_extract_variable_internal = lmp%extract_variable('int')
+END FUNCTION f_lammps_extract_variable_internal
+
+FUNCTION f_lammps_extract_variable_equal() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  REAL(c_double) :: f_lammps_extract_variable_equal
+
+  f_lammps_extract_variable_equal = lmp%extract_variable('ex')
+END FUNCTION f_lammps_extract_variable_equal
+
+FUNCTION f_lammps_extract_variable_atom(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double) :: f_lammps_extract_variable_atom
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: atom
+
+  atom = lmp%extract_variable('at_z') ! z-coordinates
+  f_lammps_extract_variable_atom = atom(i)
+END FUNCTION f_lammps_extract_variable_atom
+
+FUNCTION f_lammps_extract_variable_vector(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  REAL(c_double) :: f_lammps_extract_variable_vector
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: vector
+
+  vector = lmp%extract_variable('center') ! z-coordinates
+  f_lammps_extract_variable_vector = vector(i)
+END FUNCTION f_lammps_extract_variable_vector
+
+SUBROUTINE f_lammps_set_variable_string() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  USE keepvar, ONLY : f2c_string
+  IMPLICIT NONE
+  CHARACTER(LEN=40) :: string
+
+  string = "this is the new string"
+  CALL lmp%set_variable('str', string)
+END SUBROUTINE f_lammps_set_variable_string
+
+! vim: sts=2 ts=2 sw=2 et

unittest/fortran/test_fortran_gather_scatter.f90

@@ -0,0 +1,201 @@
+FUNCTION f_lammps_with_args() BIND(C)
+  USE 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 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_gather_scatter() BIND(C)
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp, big_input, cont_input, more_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)
+END SUBROUTINE f_lammps_setup_gather_scatter
+
+FUNCTION f_lammps_gather_atoms_mask(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  INTEGER(c_int) :: f_lammps_gather_atoms_mask
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: mask
+
+  CALL lmp%gather_atoms('mask', 1_c_int, mask)
+  f_lammps_gather_atoms_mask = mask(i)
+END FUNCTION f_lammps_gather_atoms_mask
+
+FUNCTION f_lammps_gather_atoms_position(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_atoms_position
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: positions
+
+  CALL lmp%gather_atoms('x', 3_c_int, positions)
+  f_lammps_gather_atoms_position = positions(i)
+END FUNCTION f_lammps_gather_atoms_position
+
+FUNCTION f_lammps_gather_atoms_concat_mask(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  INTEGER(c_int) :: f_lammps_gather_atoms_concat_mask
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: mask, tag
+  INTEGER :: j
+
+  CALL lmp%gather_atoms_concat('mask', 1_c_int, mask)
+  CALL lmp%gather_atoms_concat('id', 1_c_int, tag)
+  DO j = 1, SIZE(tag)
+    IF (tag(j) == i) THEN
+      f_lammps_gather_atoms_concat_mask = mask(j)
+      RETURN
+    END IF
+  END DO
+  f_lammps_gather_atoms_concat_mask = -1
+END FUNCTION f_lammps_gather_atoms_concat_mask
+
+FUNCTION f_lammps_gather_atoms_concat_position(xyz, id) 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 :: id, xyz
+  REAL(c_double) :: f_lammps_gather_atoms_concat_position
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: positions
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: tag
+  INTEGER :: j
+
+  CALL lmp%gather_atoms_concat('x', 3_c_int, positions)
+  CALL lmp%gather_atoms_concat('id', 1_c_int, tag)
+  DO j = 1, SIZE(tag)
+    IF (tag(j) == id) THEN
+       f_lammps_gather_atoms_concat_position = positions((j-1)*3 + xyz)
+    END IF
+  END DO
+END FUNCTION f_lammps_gather_atoms_concat_position
+
+FUNCTION f_lammps_gather_atoms_subset_mask(i) BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), INTENT(IN), VALUE :: i
+  INTEGER(c_int) :: f_lammps_gather_atoms_subset_mask
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: mask
+  INTEGER :: j
+  INTEGER(c_int), DIMENSION(*), PARAMETER :: tag = [3,2]
+
+  CALL lmp%gather_atoms_subset('mask', 1_c_int, tag, mask)
+  DO j = 1, SIZE(tag)
+    IF (tag(j) == i) THEN
+      f_lammps_gather_atoms_subset_mask = mask(j)
+      RETURN
+    END IF
+  END DO
+  f_lammps_gather_atoms_subset_mask = -1
+END FUNCTION f_lammps_gather_atoms_subset_mask
+
+FUNCTION f_lammps_gather_atoms_subset_position(xyz,id) 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 :: id, xyz
+  REAL(c_double) :: f_lammps_gather_atoms_subset_position
+  REAL(c_double), DIMENSION(:), ALLOCATABLE :: positions
+  INTEGER(c_int), DIMENSION(*), PARAMETER :: tag = [3,2]
+  INTEGER :: j
+
+  CALL lmp%gather_atoms_subset('x', 3_c_int, tag, positions)
+  DO j = 1, SIZE(tag)
+    IF (tag(j) == id) THEN
+      f_lammps_gather_atoms_subset_position = positions((j-1)*3 + xyz)
+      RETURN
+    END IF
+  END DO
+  f_lammps_gather_atoms_subset_position = -1.0D0
+END FUNCTION f_lammps_gather_atoms_subset_position
+
+SUBROUTINE f_lammps_scatter_atoms_masks() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: masks
+  INTEGER(c_int) :: swap
+
+  CALL lmp%gather_atoms('mask', 1_c_int, masks)
+
+  ! swap masks of atoms 1 and 3
+  swap=masks(1)
+  masks(1) = masks(3)
+  masks(3) = swap
+
+  CALL lmp%scatter_atoms('mask', masks) ! push the swap back to LAMMPS
+END SUBROUTINE f_lammps_scatter_atoms_masks
+
+SUBROUTINE f_lammps_scatter_atoms_positions() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: tags
+  REAL(c_double), DIMENSION(:), ALLOCATABLE, TARGET :: xvec
+  REAL(c_double), DIMENSION(:,:), POINTER :: x
+  REAL(c_double) :: swap(3)
+
+  CALL lmp%gather_atoms('id',1_c_int,tags)
+  CALL lmp%gather_atoms('x',3_c_int,xvec)
+  x(1:3,1:SIZE(xvec)/3) => xvec
+
+  ! swap positions of atoms 1 and 3
+  swap=x(:,1)
+  x(:,1) = x(:,3)
+  x(:,3) = swap
+
+  CALL lmp%scatter_atoms('x', xvec) ! push the swap back to LAMMPS
+END SUBROUTINE f_lammps_scatter_atoms_positions
+
+SUBROUTINE f_lammps_scatter_atoms_subset_mask() BIND(C)
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double
+  USE LIBLAMMPS
+  USE keepstuff, ONLY : lmp
+  IMPLICIT NONE
+  INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: all_masks
+  INTEGER(c_int), DIMENSION(*), PARAMETER :: tags = [3,1]
+  INTEGER(c_int), DIMENSION(2) :: masks
+
+  CALL lmp%gather_atoms('mask', 1_c_int, all_masks)
+
+  ! swap masks of atoms 1 and 3 in the new array (because 'tags' is reversed)
+  masks(1) = all_masks(1)
+  masks(2) = all_masks(3)
+
+  CALL lmp%scatter_atoms_subset('mask', tags, masks) ! push the swap to LAMMPS
+END SUBROUTINE f_lammps_scatter_atoms_subset_mask

unittest/fortran/test_fortran_properties.f90

@@ -1,40 +1,40 @@
 FUNCTION f_lammps_version() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE liblammps
    USE keepstuff, ONLY : lmp
    IMPLICIT NONE
-   INTEGER (C_int) :: f_lammps_version
+   INTEGER(c_int) :: f_lammps_version
 
    f_lammps_version = lmp%version()
 END FUNCTION f_lammps_version
 
 SUBROUTINE f_lammps_memory_usage(meminfo) BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double
    USE liblammps
    USE keepstuff, ONLY : lmp
    IMPLICIT NONE
-   REAL (C_double), DIMENSION(3), INTENT(OUT) :: meminfo
+   REAL(c_double), DIMENSION(3), INTENT(OUT) :: meminfo
 
    CALL lmp%memory_usage(meminfo)
 END SUBROUTINE f_lammps_memory_usage
 
 FUNCTION f_lammps_get_mpi_comm() BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE liblammps
    USE keepstuff, ONLY : lmp
    IMPLICIT NONE
-   INTEGER (C_int) :: f_lammps_get_mpi_comm
+   INTEGER(c_int) :: f_lammps_get_mpi_comm
 
    f_lammps_get_mpi_comm = lmp%get_mpi_comm()
 END FUNCTION f_lammps_get_mpi_comm
 
 FUNCTION f_lammps_extract_setting(Cstr) BIND(C)
-   USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int, C_char
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_char
    USE keepstuff, ONLY : lmp
    USE LIBLAMMPS
    IMPLICIT NONE
-   INTEGER (C_int) :: f_lammps_extract_setting
+   INTEGER(c_int) :: f_lammps_extract_setting
-   CHARACTER (KIND=C_char, LEN=1), DIMENSION(*), INTENT(IN) :: Cstr
+   CHARACTER(KIND=c_char, LEN=1), DIMENSION(*), INTENT(IN) :: Cstr
    INTEGER :: strlen, i
    CHARACTER(LEN=:), ALLOCATABLE :: Fstr
 
@@ -43,10 +43,47 @@ FUNCTION f_lammps_extract_setting (Cstr) BIND(C)
       i = i + 1
    END DO
    strlen = i
-   allocate ( CHARACTER(LEN=strlen) :: Fstr)
+   ALLOCATE(CHARACTER(LEN=strlen) :: Fstr)
-   FORALL (i=1:strlen)
+   DO i = 1, strlen
       Fstr(i:i) = Cstr(i)
-   END FORALL
+   END DO
    f_lammps_extract_setting = lmp%extract_setting(Fstr)
-   deallocate (Fstr)
+   DEALLOCATE(Fstr)
 END FUNCTION f_lammps_extract_setting
+
+FUNCTION f_lammps_has_error() BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+   USE keepstuff, ONLY : lmp
+   USE LIBLAMMPS
+   IMPLICIT NONE
+   INTEGER(c_int) :: f_lammps_has_error
+
+   IF (lmp%has_error()) THEN
+      f_lammps_has_error = 1_c_int
+   ELSE
+      f_lammps_has_error = 0_c_int
+   END IF
+END FUNCTION f_lammps_has_error
+
+FUNCTION f_lammps_get_last_error_message(errmesg, errlen) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_char, c_ptr, C_F_POINTER
+   USE keepstuff, ONLY : lmp
+   USE LIBLAMMPS
+   IMPLICIT NONE
+   INTEGER(c_int) :: f_lammps_get_last_error_message
+   CHARACTER(KIND=c_char), DIMENSION(*) :: errmesg
+   INTEGER(c_int), VALUE, INTENT(IN) :: errlen
+   CHARACTER(LEN=:), ALLOCATABLE :: buffer
+   INTEGER :: status, i
+
+   ! copy error message to buffer
+   ALLOCATE(CHARACTER(errlen) :: buffer)
+   CALL lmp%get_last_error_message(buffer, status)
+   f_lammps_get_last_error_message = status
+   ! and copy to C style string
+   DO i=1, errlen
+      errmesg(i) = buffer(i:i)
+      IF (buffer(i:i) == ACHAR(0)) EXIT
+   END DO
+   DEALLOCATE(buffer)
+END FUNCTION f_lammps_get_last_error_message

unittest/fortran/wrap_extract_atom.cpp

@@ -0,0 +1,120 @@
+// unit tests for extracting Atom class 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_extract_atom();
+double f_lammps_extract_atom_mass();
+int f_lammps_extract_atom_tag_int(int);
+int64_t f_lammps_extract_atom_tag_int64(int64_t);
+int f_lammps_extract_atom_type(int);
+int f_lammps_extract_atom_mask(int);
+void f_lammps_extract_atom_x(int, double *);
+void f_lammps_extract_atom_v(int, double *);
+}
+
+class LAMMPS_extract_atom : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_extract_atom()           = default;
+    ~LAMMPS_extract_atom() 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_extract_atom, mass)
+{
+    f_lammps_setup_extract_atom();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_atom_mass(), 2.0);
+};
+
+TEST_F(LAMMPS_extract_atom, tag)
+{
+    f_lammps_setup_extract_atom();
+#if defined(LAMMPS_BIGBIG)
+    EXPECT_EQ(f_lammps_extract_atom_tag_int64(1l), 1l);
+    EXPECT_EQ(f_lammps_extract_atom_tag_int64(2l), 2l);
+#else
+    EXPECT_EQ(f_lammps_extract_atom_tag_int(1), 1);
+    EXPECT_EQ(f_lammps_extract_atom_tag_int(2), 2);
+#endif
+};
+
+TEST_F(LAMMPS_extract_atom, type)
+{
+    f_lammps_setup_extract_atom();
+    EXPECT_EQ(f_lammps_extract_atom_type(1), 1);
+    EXPECT_EQ(f_lammps_extract_atom_type(2), 1);
+};
+
+TEST_F(LAMMPS_extract_atom, mask)
+{
+    f_lammps_setup_extract_atom();
+    EXPECT_EQ(f_lammps_extract_atom_mask(1), 1);
+    EXPECT_EQ(f_lammps_extract_atom_mask(2), 1);
+    lammps_command(lmp, "group 1 id 1");
+    lammps_command(lmp, "group 2 id 2");
+    EXPECT_EQ(f_lammps_extract_atom_mask(1), 3);
+    EXPECT_EQ(f_lammps_extract_atom_mask(2), 5);
+};
+
+TEST_F(LAMMPS_extract_atom, x)
+{
+    f_lammps_setup_extract_atom();
+    double x1[3];
+    double x2[3];
+    f_lammps_extract_atom_x(1, x1);
+    EXPECT_DOUBLE_EQ(x1[0], 1.0);
+    EXPECT_DOUBLE_EQ(x1[1], 1.0);
+    EXPECT_DOUBLE_EQ(x1[2], 1.5);
+    f_lammps_extract_atom_x(2, x2);
+    EXPECT_DOUBLE_EQ(x2[0], 0.2);
+    EXPECT_DOUBLE_EQ(x2[1], 0.1);
+    EXPECT_DOUBLE_EQ(x2[2], 0.1);
+}
+
+TEST_F(LAMMPS_extract_atom, v)
+{
+    f_lammps_setup_extract_atom();
+    double v1[3];
+    double v2[3];
+    f_lammps_extract_atom_v(1, v1);
+    EXPECT_DOUBLE_EQ(v1[0], 0.0);
+    EXPECT_DOUBLE_EQ(v1[1], 0.0);
+    EXPECT_DOUBLE_EQ(v1[2], 0.0);
+    f_lammps_extract_atom_v(2, v2);
+    EXPECT_DOUBLE_EQ(v2[0], 0.0);
+    EXPECT_DOUBLE_EQ(v2[1], 0.0);
+    EXPECT_DOUBLE_EQ(v2[2], 0.0);
+    lammps_command(lmp, "group one id 1");
+    lammps_command(lmp, "velocity one set 1 2 3");
+    f_lammps_extract_atom_v(1, v1);
+    EXPECT_DOUBLE_EQ(v1[0], 1.0);
+    EXPECT_DOUBLE_EQ(v1[1], 2.0);
+    EXPECT_DOUBLE_EQ(v1[2], 3.0);
+}

unittest/fortran/wrap_extract_compute.cpp

@@ -0,0 +1,168 @@
+// unit tests for extracting compute 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_extract_compute();
+double f_lammps_extract_compute_peratom_vector(int);
+double f_lammps_extract_compute_peratom_array(int, int);
+double f_lammps_extract_compute_global_scalar();
+double f_lammps_extract_compute_global_vector(int);
+double f_lammps_extract_compute_global_array(int, int);
+double f_lammps_extract_compute_local_vector(int);
+double f_lammps_extract_compute_local_array(int, int);
+}
+
+class LAMMPS_extract_compute : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_extract_compute()           = default;
+    ~LAMMPS_extract_compute() 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_extract_compute, peratom_vector)
+{
+    f_lammps_setup_extract_compute();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_vector(1), -0.599703102447981);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_vector(2), 391.817623795857);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_vector(3), 391.430665759871);
+};
+
+TEST_F(LAMMPS_extract_compute, peratom_array)
+{
+    f_lammps_setup_extract_compute();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(1, 1), 0.8837067009319107);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(2, 1), 0.3588584939803668);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(3, 1), 1.2799807127711049);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(4, 1), 0.20477632346642258);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(5, 1), 0.400429511840588);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(6, 1), 0.673995757699694);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(1, 2), -1070.0291234709418);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(2, 2), -1903.651817128683);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(3, 2), -1903.5121520875714);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(4, 2), -1427.867483013);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(5, 2), -1427.8560790941347);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_peratom_array(6, 2), -1903.5971655908565);
+};
+
+TEST_F(LAMMPS_extract_compute, global_scalar)
+{
+    f_lammps_setup_extract_compute();
+    double *scalar;
+    scalar = (double *)lammps_extract_compute(lmp, "totalpe", LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR);
+    // EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_scalar(), 782.64858645328);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_scalar(), *scalar);
+};
+
+TEST_F(LAMMPS_extract_compute, global_vector)
+{
+    f_lammps_setup_extract_compute();
+    double *vector;
+    vector = (double *)lammps_extract_compute(lmp, "COM", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_vector(1), vector[0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_vector(2), vector[1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_vector(3), vector[2]);
+};
+
+TEST_F(LAMMPS_extract_compute, global_array)
+{
+    f_lammps_setup_extract_compute();
+    double **array;
+    array = (double **)lammps_extract_compute(lmp, "RDF", LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(1, 1), array[0][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(2, 1), array[0][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(1, 2), array[1][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(2, 2), array[1][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(1, 3), array[2][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_global_array(1, 4), array[3][0]);
+};
+TEST_F(LAMMPS_extract_compute, local_vector)
+{
+    f_lammps_setup_extract_compute();
+    double *vector;
+    vector = (double *)lammps_extract_compute(lmp, "pairdist", LMP_STYLE_LOCAL, LMP_TYPE_VECTOR);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(1), vector[0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(2), vector[1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(3), vector[2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(4), vector[3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(5), vector[4]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(6), vector[5]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(7), vector[6]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(8), vector[7]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(9), vector[8]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_vector(10), vector[9]);
+};
+
+TEST_F(LAMMPS_extract_compute, local_array)
+{
+    f_lammps_setup_extract_compute();
+    double **array;
+    array = (double **)lammps_extract_compute(lmp, "pairlocal", LMP_STYLE_LOCAL, LMP_TYPE_ARRAY);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 1), array[0][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 1), array[0][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 1), array[0][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 1), array[0][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 2), array[1][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 2), array[1][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 2), array[1][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 2), array[1][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 3), array[2][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 3), array[2][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 3), array[2][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 3), array[2][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 4), array[3][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 4), array[3][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 4), array[3][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 4), array[3][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 5), array[4][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 5), array[4][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 5), array[4][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 5), array[4][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 6), array[5][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 6), array[5][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 6), array[5][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 6), array[5][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 7), array[6][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 7), array[6][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 7), array[6][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 7), array[6][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 8), array[7][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 8), array[7][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 8), array[7][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 8), array[7][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 9), array[8][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 9), array[8][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 9), array[8][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 9), array[8][3]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(1, 10), array[9][0]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(2, 10), array[9][1]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(3, 10), array[9][2]);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_compute_local_array(4, 10), array[9][3]);
+};

unittest/fortran/wrap_extract_fix.cpp

@@ -0,0 +1,107 @@
+// unit tests for extracting compute data from a LAMMPS instance through the
+// Fortran wrapper
+#include <cstdio>
+
+#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_extract_fix();
+double f_lammps_extract_fix_global_scalar();
+double f_lammps_extract_fix_global_vector(int);
+double f_lammps_extract_fix_global_array(int, int);
+double f_lammps_extract_fix_peratom_vector(int);
+double f_lammps_extract_fix_peratom_array(int, int);
+double f_lammps_extract_fix_local_vector(int);
+double f_lammps_extract_fix_local_array(int, int);
+}
+
+class LAMMPS_extract_fix : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_extract_fix()           = default;
+    ~LAMMPS_extract_fix() 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_extract_fix, global_scalar)
+{
+    f_lammps_setup_extract_fix();
+    double *scalar =
+        (double *)lammps_extract_fix(lmp, "recenter", LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR, -1, -1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_scalar(), *scalar);
+    lammps_free(scalar);
+};
+
+TEST_F(LAMMPS_extract_fix, global_vector)
+{
+    f_lammps_setup_extract_fix();
+    double *x =
+        (double *)lammps_extract_fix(lmp, "recenter", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR, 0, -1);
+    double *y =
+        (double *)lammps_extract_fix(lmp, "recenter", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR, 1, -1);
+    double *z =
+        (double *)lammps_extract_fix(lmp, "recenter", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR, 2, -1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_vector(1), *x);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_vector(2), *y);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_vector(3), *z);
+    lammps_free(x);
+    lammps_free(y);
+    lammps_free(z);
+};
+
+TEST_F(LAMMPS_extract_fix, global_array)
+{
+    f_lammps_setup_extract_fix();
+    double natoms = lammps_get_natoms(lmp);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_array(1, 1), natoms);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_array(1, 2), natoms);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_array(2, 1), 0.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_global_array(2, 2), 1.0);
+};
+
+TEST_F(LAMMPS_extract_fix, peratom_vector)
+{
+    f_lammps_setup_extract_fix();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_vector(1), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_vector(2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_vector(3), 0.5);
+};
+
+TEST_F(LAMMPS_extract_fix, peratom_array)
+{
+    f_lammps_setup_extract_fix();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(1, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(2, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(3, 1), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(1, 2), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(2, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(3, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(1, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(2, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_fix_peratom_array(3, 3), 0.5);
+};

unittest/fortran/wrap_extract_global.cpp

@@ -3,10 +3,10 @@
 
 #include "lammps.h"
 #include "library.h"
+#include <cstdint>
+#include <cstdlib>
 #include <mpi.h>
 #include <string>
-#include <cstdlib>
-#include <cstdint>
 
 #include "gtest/gtest.h"
 

unittest/fortran/wrap_extract_variable.cpp

@@ -0,0 +1,283 @@
+// unit tests for extracting compute data from a LAMMPS instance through the
+// Fortran wrapper
+
+#include "lammps.h"
+#include "library.h"
+#include "platform.h"
+
+#include <chrono>
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <mpi.h>
+#include <string>
+#include <thread>
+
+#include "gtest/gtest.h"
+
+#define STRINGIFY(val) XSTR(val)
+#define XSTR(val) #val
+
+// prototypes for Fortran reverse wrapper functions
+extern "C" {
+void *f_lammps_with_c_args(int, char **);
+void f_lammps_close();
+void f_lammps_setup_extract_variable();
+int f_lammps_extract_variable_index_1();
+int f_lammps_extract_variable_index_2();
+int f_lammps_extract_variable_loop();
+char *f_lammps_extract_variable_loop_pad();
+char *f_lammps_extract_variable_world();
+char *f_lammps_extract_variable_universe();
+int f_lammps_extract_variable_uloop();
+char *f_lammps_extract_variable_string();
+char *f_lammps_extract_variable_format();
+char *f_lammps_extract_variable_format_pad();
+char *f_lammps_extract_variable_getenv();
+char *f_lammps_extract_variable_file();
+double f_lammps_extract_variable_atomfile(int);
+double f_lammps_extract_variable_python();
+double f_lammps_extract_variable_timer();
+double f_lammps_extract_variable_internal();
+double f_lammps_extract_variable_equal();
+double f_lammps_extract_variable_atom(int);
+double f_lammps_extract_variable_vector(int);
+void f_lammps_set_variable_string();
+char *c_path_join(const char *, const char *);
+}
+
+char *c_path_join(const char *a, const char *b)
+{
+    std::string A = a;
+    std::string B = b;
+    std::string C = LAMMPS_NS::platform::path_join(A, B);
+    size_t length = C.length() + 1;
+    char *retval  = (char *)malloc(length * sizeof(char));
+    C.copy(retval, length);
+    retval[length - 1] = '\0';
+    return retval;
+}
+
+constexpr char input_dir[] = STRINGIFY(TEST_INPUT_FOLDER);
+class LAMMPS_extract_variable : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_extract_variable()           = default;
+    ~LAMMPS_extract_variable() override = default;
+
+    void SetUp() override
+    {
+        // clang-format off
+        const char *args[] =
+            { "LAMMPS_Fortran_test", "-l", "none", "-echo", "screen", "-nocite",
+              "-var", "input_dir", input_dir, "-var", "zpos", "1.5", "-var", "x", "2" };
+        // clang-format on
+        char **argv = (char **)args;
+        int argc    = sizeof(args) / sizeof(const char *);
+        ::testing::internal::CaptureStdout();
+        lmp = (LAMMPS_NS::LAMMPS *)f_lammps_with_c_args(argc, argv);
+
+        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_extract_variable, index)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_EQ(f_lammps_extract_variable_index_1(), 1);
+    EXPECT_EQ(f_lammps_extract_variable_index_2(), 0);
+    lammps_command(lmp, "next idx");
+    EXPECT_EQ(f_lammps_extract_variable_index_1(), 0);
+    EXPECT_EQ(f_lammps_extract_variable_index_2(), 1);
+};
+
+TEST_F(LAMMPS_extract_variable, loop)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    for (i = 1; i <= 10; i++) {
+        EXPECT_EQ(f_lammps_extract_variable_loop(), i);
+        lammps_command(lmp, "next lp");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, loop_pad)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    char str[10];
+    char *fstr;
+    for (i = 1; i <= 10; i++) {
+        std::sprintf(str, "%02d", i);
+        fstr = f_lammps_extract_variable_loop_pad();
+        EXPECT_STREQ(fstr, str);
+        std::free(fstr);
+        lammps_command(lmp, "next lp_pad");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, world)
+{
+    f_lammps_setup_extract_variable();
+    char *fstr = f_lammps_extract_variable_world();
+    EXPECT_STREQ(fstr, "group1");
+    std::free(fstr);
+};
+
+TEST_F(LAMMPS_extract_variable, universe)
+{
+    f_lammps_setup_extract_variable();
+    char *fstr = f_lammps_extract_variable_universe();
+    EXPECT_STREQ(fstr, "universe1");
+    std::free(fstr);
+};
+
+TEST_F(LAMMPS_extract_variable, uloop)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_EQ(f_lammps_extract_variable_uloop(), 1);
+};
+
+TEST_F(LAMMPS_extract_variable, string)
+{
+    f_lammps_setup_extract_variable();
+    char *fstr = f_lammps_extract_variable_string();
+    EXPECT_STREQ(fstr, "this is a string");
+    std::free(fstr);
+    f_lammps_set_variable_string();
+    fstr = f_lammps_extract_variable_string();
+    EXPECT_STREQ(fstr, "this is the new string");
+    std::free(fstr);
+};
+
+TEST_F(LAMMPS_extract_variable, format)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    char str[16];
+    char *fstr;
+    for (i = 1; i <= 10; i++) {
+        std::sprintf(str, "%.6G", std::exp(i));
+        fstr = f_lammps_extract_variable_format();
+        EXPECT_STREQ(fstr, str);
+        std::free(fstr);
+        lammps_command(lmp, "next lp");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, format_pad)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    char str[16];
+    char *fstr;
+    for (i = 1; i <= 10; i++) {
+        std::sprintf(str, "%08.6G", std::exp(i));
+        fstr = f_lammps_extract_variable_format_pad();
+        EXPECT_STREQ(fstr, str);
+        std::free(fstr);
+        lammps_command(lmp, "next lp");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, getenv)
+{
+    LAMMPS_NS::platform::putenv("FORTRAN_USER=myuser");
+    f_lammps_setup_extract_variable();
+    char *env  = std::getenv("FORTRAN_USER");
+    char *fenv = f_lammps_extract_variable_getenv();
+    EXPECT_STREQ(fenv, env);
+    std::free(fenv);
+};
+
+TEST_F(LAMMPS_extract_variable, file)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    const char *str[9] = {"hello",      "god_dag", "hola",  "bonjour",  "guten_Tag",
+                          "konnichiwa", "shalom",  "salve", "goedendag"};
+    char *fstr;
+    for (i = 0; i < 9; i++) {
+        fstr = f_lammps_extract_variable_file();
+        EXPECT_STREQ(fstr, str[i]);
+        std::free(fstr);
+        lammps_command(lmp, "next greeting");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, atomfile)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(1), 5.2);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(2), 1.6);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(3), -1.4);
+    lammps_command(lmp, "next atfile");
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(1), -1.1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(2), 0.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atomfile(3), 2.5);
+};
+
+TEST_F(LAMMPS_extract_variable, python)
+{
+    if (lammps_config_has_package("PYTHON")) {
+        f_lammps_setup_extract_variable();
+        for (int i = 1; i <= 10; i++) {
+            EXPECT_DOUBLE_EQ(f_lammps_extract_variable_python(), static_cast<double>(i * i));
+            lammps_command(lmp, "next lp");
+        }
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, timer)
+{
+    f_lammps_setup_extract_variable();
+    double initial_t, final_t;
+    initial_t = f_lammps_extract_variable_timer();
+    std::this_thread::sleep_for(std::chrono::milliseconds(100));
+    lammps_command(lmp, "variable time timer"); // update the time
+    final_t = f_lammps_extract_variable_timer();
+    EXPECT_GT(final_t, initial_t + 0.1);
+};
+
+TEST_F(LAMMPS_extract_variable, internal)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_internal(), 4.0);
+};
+
+TEST_F(LAMMPS_extract_variable, equal)
+{
+    f_lammps_setup_extract_variable();
+    int i;
+    for (i = 1; i <= 10; i++) {
+        EXPECT_DOUBLE_EQ(f_lammps_extract_variable_equal(), std::exp(static_cast<double>(i)));
+        lammps_command(lmp, "next lp");
+    }
+};
+
+TEST_F(LAMMPS_extract_variable, atom)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atom(1), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atom(2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_atom(3), 0.5);
+};
+
+TEST_F(LAMMPS_extract_variable, vector)
+{
+    f_lammps_setup_extract_variable();
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_vector(1), (1 + 0.2 + 0.5) / 3.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_vector(2), (1 + 0.1 + 0.5) / 3.0);
+    EXPECT_DOUBLE_EQ(f_lammps_extract_variable_vector(3), (1.5 + 0.1 + 0.5) / 3.0);
+};

unittest/fortran/wrap_gather_scatter.cpp

@@ -0,0 +1,202 @@
+// 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_gather_scatter();
+int f_lammps_gather_atoms_mask(int);
+double f_lammps_gather_atoms_position(int);
+int f_lammps_gather_atoms_concat_mask(int);
+double f_lammps_gather_atoms_concat_position(int, int);
+int f_lammps_gather_atoms_subset_mask(int);
+double f_lammps_gather_atoms_subset_position(int, int);
+void f_lammps_scatter_atoms_masks();
+void f_lammps_scatter_atoms_positions();
+}
+
+class LAMMPS_gather_scatter : public ::testing::Test {
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS_gather_scatter()           = default;
+    ~LAMMPS_gather_scatter() 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_gather_scatter, gather_atoms_masks)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 1);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(2), 1);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 1);
+    lammps_command(lmp, "group special id 1");
+    lammps_command(lmp, "group other id 2");
+    lammps_command(lmp, "group spiffy id 3");
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 3);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(2), 5);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 9);
+    lammps_command(lmp, "group other id 1");
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 7);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(2), 5);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 9);
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_atoms_positions)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(2), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(3), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(4), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(5), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(6), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(7), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(8), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_position(9), 0.5);
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_atoms_concat_masks)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(1), 1);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(2), 1);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(3), 1);
+    lammps_command(lmp, "group special id 1");
+    lammps_command(lmp, "group other id 2");
+    lammps_command(lmp, "group spiffy id 3");
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(1), 3);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(2), 5);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(3), 9);
+    lammps_command(lmp, "group other id 1");
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(1), 7);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(2), 5);
+    EXPECT_EQ(f_lammps_gather_atoms_concat_mask(3), 9);
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_atoms_concat_positions)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 1), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 2), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 3), 0.5);
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_atoms_subset_masks)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(2), 1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(3), 1);
+    lammps_command(lmp, "group special id 1");
+    lammps_command(lmp, "group other id 2");
+    lammps_command(lmp, "group spiffy id 3");
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(2), 5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(3), 9);
+    lammps_command(lmp, "group other id 3");
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(2), 5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_mask(3), 13);
+};
+
+TEST_F(LAMMPS_gather_scatter, gather_atoms_subset_positions)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(1, 2), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(2, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(3, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(1, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(2, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_subset_position(3, 3), 0.5);
+};
+
+TEST_F(LAMMPS_gather_scatter, scatter_atoms_masks)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    lammps_command(lmp, "group special id 1");
+    lammps_command(lmp, "group other id 2");
+    lammps_command(lmp, "group spiffy id 3");
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(1), 3);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(2), 5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(3), 9);
+    f_lammps_scatter_atoms_masks();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(1), 9);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(2), 5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_mask(3), 3);
+};
+
+TEST_F(LAMMPS_gather_scatter, scatter_atoms_positions)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 1), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 1), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 2), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 3), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 3), 0.5);
+    f_lammps_scatter_atoms_positions();
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 3), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 3), 1.0);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 3), 1.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 2), 0.2);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 2), 0.1);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(1, 1), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(2, 1), 0.5);
+    EXPECT_DOUBLE_EQ(f_lammps_gather_atoms_concat_position(3, 1), 0.5);
+};
+
+TEST_F(LAMMPS_gather_scatter, scatter_atoms_subset_mask)
+{
+    if (lammps_extract_setting(nullptr, "tagint") == 8) GTEST_SKIP();
+    f_lammps_setup_gather_scatter();
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 1);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 1);
+    lammps_command(lmp, "group special id 1");
+    lammps_command(lmp, "group other id 2");
+    lammps_command(lmp, "group spiffy id 3");
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 3);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 9);
+    f_lammps_scatter_atoms_masks();
+    EXPECT_EQ(f_lammps_gather_atoms_mask(1), 9);
+    EXPECT_EQ(f_lammps_gather_atoms_mask(3), 3);
+};

unittest/fortran/wrap_properties.cpp

@@ -2,9 +2,10 @@
 
 #include "lammps.h"
 #include "library.h"
-#include <mpi.h>
+
 #include <string>
 
+#include "gmock/gmock.h"
 #include "gtest/gtest.h"
 
 // prototypes for fortran reverse wrapper functions
@@ -15,21 +16,27 @@ int f_lammps_version();
 void f_lammps_memory_usage(double *);
 int f_lammps_get_mpi_comm();
 int f_lammps_extract_setting(const char *);
+int f_lammps_has_error();
+int f_lammps_get_last_error_message(char *, int);
 }
 
+namespace LAMMPS_NS {
+
+using ::testing::ContainsRegex;
+
 class LAMMPS_properties : public ::testing::Test {
 protected:
-    LAMMPS_NS::LAMMPS *lmp;
+    LAMMPS *lmp;
-    LAMMPS_properties()           = default;
-    ~LAMMPS_properties() override = default;
 
     void SetUp() override
     {
         ::testing::internal::CaptureStdout();
-        lmp                = (LAMMPS_NS::LAMMPS *)f_lammps_with_args();
+        lmp = (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();
@@ -104,5 +111,34 @@ TEST_F(LAMMPS_properties, extract_setting)
     EXPECT_EQ(f_lammps_extract_setting("mu_flag"), 0);
     EXPECT_EQ(f_lammps_extract_setting("rmass_flag"), 0);
     EXPECT_EQ(f_lammps_extract_setting("UNKNOWN"), -1);
-
 };
+
+TEST_F(LAMMPS_properties, has_error)
+{
+    // need errors to throw exceptions to be able to intercept them.
+    if (!lammps_config_has_exceptions()) GTEST_SKIP();
+
+    EXPECT_EQ(f_lammps_has_error(), lammps_has_error(lmp));
+    EXPECT_EQ(f_lammps_has_error(), 0);
+
+    // trigger an error, but hide output
+    ::testing::internal::CaptureStdout();
+    lammps_command(lmp, "this_is_not_a_known_command");
+    ::testing::internal::GetCapturedStdout();
+
+    EXPECT_EQ(f_lammps_has_error(), lammps_has_error(lmp));
+    EXPECT_EQ(f_lammps_has_error(), 1);
+
+    // retrieve error message
+    char errmsg[1024];
+    int err = f_lammps_get_last_error_message(errmsg, 1023);
+    EXPECT_EQ(err, 1);
+    EXPECT_THAT(errmsg, ContainsRegex(".*ERROR: Unknown command: this_is_not_a_known_command.*"));
+
+    // retrieving the error message clear the error status
+    EXPECT_EQ(f_lammps_has_error(), 0);
+    err = f_lammps_get_last_error_message(errmsg, 1023);
+    EXPECT_EQ(err, 0);
+    EXPECT_THAT(errmsg, ContainsRegex("                                                  "));
+};
+} // namespace LAMMPS_NS