diff --git a/cmake/Modules/CodeCoverage.cmake b/cmake/Modules/CodeCoverage.cmake
index 3b323e37ff..054e08fc1a 100644
--- a/cmake/Modules/CodeCoverage.cmake
+++ b/cmake/Modules/CodeCoverage.cmake
@@ -3,11 +3,16 @@
 #
 # Requires latest gcovr (for GCC 8.1 support):#
 # pip install git+https://github.com/gcovr/gcovr.git
+#
+# For Python coverage the coverage package needs to be installed
 ###############################################################################
 if(ENABLE_COVERAGE)
     find_program(GCOVR_BINARY gcovr)
     find_package_handle_standard_args(GCOVR DEFAULT_MSG GCOVR_BINARY)
 
+    find_program(COVERAGE_BINARY coverage)
+    find_package_handle_standard_args(COVERAGE DEFAULT_MSG COVERAGE_BINARY)
+
     if(GCOVR_FOUND)
         get_filename_component(ABSOLUTE_LAMMPS_SOURCE_DIR ${LAMMPS_SOURCE_DIR} ABSOLUTE)
 
@@ -46,4 +51,30 @@ if(ENABLE_COVERAGE)
         )
         add_dependencies(reset_coverage clean_coverage_html)
     endif()
+
+    if(COVERAGE_FOUND)
+        set(PYTHON_COVERAGE_HTML_DIR ${CMAKE_BINARY_DIR}/python_coverage_html)
+        add_custom_command(
+            OUTPUT ${CMAKE_BINARY_DIR}/unittest/python/.coverage
+            COMMAND ${COVERAGE_BINARY} combine
+            WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/unittest/python
+            COMMENT "Combine Python coverage files..."
+        )
+
+        add_custom_target(
+            gen_python_coverage_html
+            COMMAND ${COVERAGE_BINARY} html -d ${PYTHON_COVERAGE_HTML_DIR}
+            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage
+            WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/unittest/python
+            COMMENT "Generating HTML Python coverage report..."
+        )
+
+        add_custom_target(
+            gen_python_coverage_xml
+            COMMAND ${COVERAGE_BINARY} xml -o ${CMAKE_BINARY_DIR}/python_coverage.xml
+            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage
+            WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/unittest/python
+            COMMENT "Generating XML Python coverage report..."
+        )
+    endif()
 endif()
diff --git a/doc/src/pg_lib_objects.rst b/doc/src/pg_lib_objects.rst
index 5da858db0c..301511b848 100644
--- a/doc/src/pg_lib_objects.rst
+++ b/doc/src/pg_lib_objects.rst
@@ -26,6 +26,8 @@ computes, fixes, or variables in LAMMPS.
 
 -----------------------
 
+.. doxygenenum:: _LMP_DATATYPE_CONST
+
 .. doxygenenum:: _LMP_STYLE_CONST
 
 .. doxygenenum:: _LMP_TYPE_CONST
diff --git a/doc/src/pg_lib_properties.rst b/doc/src/pg_lib_properties.rst
index 256ed1cf82..bf36dcb8b1 100644
--- a/doc/src/pg_lib_properties.rst
+++ b/doc/src/pg_lib_properties.rst
@@ -84,11 +84,22 @@ event as atoms are migrating between sub-domains.
 
 -----------------------
 
+.. doxygenfunction:: lammps_extract_global_datatype
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_extract_global
    :project: progguide
 
 -----------------------
 
+.. doxygenfunction:: lammps_extract_atom_datatype
+   :project: progguide
+
+
+-----------------------
+
 .. doxygenfunction:: lammps_extract_atom
    :project: progguide
 
diff --git a/doc/src/pg_python.rst b/doc/src/pg_python.rst
index 80e83d4638..87310d43e7 100644
--- a/doc/src/pg_python.rst
+++ b/doc/src/pg_python.rst
@@ -28,16 +28,58 @@ There are multiple Python interface classes in the :py:mod:`lammps` module:
 
 ----------
 
+Setting up a Python virtual environment
+***************************************
+
+LAMMPS and its Python module can be installed together into a Python virtual
+environment. This lets you isolate your customized Python environment from
+your user or system installation. The following is a minimal working example:
+
+.. code-block:: bash
+
+   # create and change into build directory
+   mkdir build
+   cd build
+
+   # create virtual environment
+   virtualenv myenv
+
+   # Add venv lib folder to LD_LIBRARY_PATH when activating it
+   echo 'export LD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$LD_LIBRARY_PATH' >> myenv/bin/activate
+
+   # Add LAMMPS_POTENTIALS path when activating venv
+   echo 'export LAMMPS_POTENTIALS=$VIRTUAL_ENV/share/lammps/potentials' >> myenv/bin/activate
+
+   # activate environment
+   source myenv/bin/activate
+
+   # configure LAMMPS compilation
+   # compiles as shared library with PYTHON package and C++ exceptions
+   # and installs into myvenv
+   (myenv)$ cmake -C ../cmake/presets/minimal.cmake    \
+                  -D BUILD_SHARED_LIBS=on              \
+                  -D PKG_PYTHON=on                     \
+                  -D LAMMPS_EXCEPTIONS=on              \
+                  -D CMAKE_INSTALL_PREFIX=$VIRTUAL_ENV \
+                  ../cmake
+
+    # compile LAMMPS
+    (myenv)$ cmake --build . --parallel
+
+    # install LAMMPS into myvenv
+    (myenv)$ cmake --install .
+
 Creating or deleting a LAMMPS object
 ************************************
 
 With the Python interface the creation of a :cpp:class:`LAMMPS
-<LAMMPS_NS::LAMMPS>` instance is included in the constructor for the
-:py:func:`lammps <lammps.lammps>` class.  Internally it will call either
-:cpp:func:`lammps_open` or :cpp:func:`lammps_open_no_mpi` from the C
+<LAMMPS_NS::LAMMPS>` instance is included in the constructors for the
+:py:meth:`lammps <lammps.lammps.__init__()>`, :py:meth:`PyLammps <lammps.PyLammps.__init__()>`,
+and :py:meth:`PyLammps <lammps.IPyLammps.__init__()>` classes.
+Internally it will call either :cpp:func:`lammps_open` or :cpp:func:`lammps_open_no_mpi` from the C
 library API to create the class instance.
 
-All arguments are optional.  The *name* argument is to allow loading a
+All arguments are optional.  The *name* argument allows loading a
 LAMMPS shared library that is named ``liblammps_machine.so`` instead of
 the default name of ``liblammps.so``.  In most cases the latter will be
 installed or used.  The *ptr* argument is for use of the
@@ -48,22 +90,99 @@ to the Python class and used instead of creating a new instance.  The
 *comm* argument may be used in combination with the `mpi4py <mpi4py_url_>`_
 module to pass an MPI communicator to LAMMPS and thus it is possible
 to run the Python module like the library interface on a subset of the
-MPI ranks after splitting the communicator. Here is a simple example:
+MPI ranks after splitting the communicator.
 
-.. code-block:: python
 
-   from lammps import lammps
+Here are simple examples using all three Python interfaces:
 
-   # NOTE: argv[0] is set by the Python module
-   args = ["-log", "none"]
-   # create LAMMPS instance
-   lmp = lammps(cmdargs=args)
-   # get and print numerical version code
-   print("LAMMPS Version: ", lmp.version())
-   # explicitly close and delete LAMMPS instance (optional)
-   lmp.close()
+.. tabs::
 
-Same as with the :ref:`C library API <lammps_c_api>` this will use the
+   .. tab:: lammps API
+
+      .. code-block:: python
+
+         from lammps import lammps
+
+         # NOTE: argv[0] is set by the lammps class constructor
+         args = ["-log", "none"]
+         # create LAMMPS instance
+         lmp = lammps(cmdargs=args)
+         # get and print numerical version code
+         print("LAMMPS Version: ", lmp.version())
+         # explicitly close and delete LAMMPS instance (optional)
+         lmp.close()
+
+   .. tab:: PyLammps API
+
+      The :py:class:`PyLammps` class is a wrapper around the
+      :py:class:`lammps` class and all of its lower level functions.
+      By default, it will create a new instance of :py:class:`lammps` passing
+      along all arguments to the constructor of :py:class:`lammps`.
+
+      .. code-block:: python
+
+         from lammps import PyLammps
+
+         # NOTE: argv[0] is set by the lammps class constructor
+         args = ["-log", "none"]
+         # create LAMMPS instance
+         L = PyLammps(cmdargs=args)
+         # get and print numerical version code
+         print("LAMMPS Version: ", L.version())
+         # explicitly close and delete LAMMPS instance (optional)
+         L.close()
+
+      :py:class:`PyLammps` objects can also be created on top of an existing :py:class:`lammps` object:
+
+      .. code-block:: Python
+
+         from lammps import lammps, PyLammps
+         ...
+         # create LAMMPS instance
+         lmp = lammps(cmdargs=args)
+         # create PyLammps instance using previously created LAMMPS instance
+         L = PyLammps(ptr=lmp)
+
+      This is useful if you have to create the :py:class:`lammps <lammps.lammps>`
+      instance is a specific way, but want to take advantage of the
+      :py:class:`PyLammps <lammps.PyLammps>` interface.
+
+   .. tab:: IPyLammps API
+
+      The :py:class:`IPyLammps` class is an extension of the
+      :py:class:`PyLammps` class. It has the same construction behavior. By
+      default, it will create a new instance of :py:class:`lammps` passing
+      along all arguments to the constructor of :py:class:`lammps`.
+
+      .. code-block:: python
+
+         from lammps import IPyLammps
+
+         # NOTE: argv[0] is set by the lammps class constructor
+         args = ["-log", "none"]
+         # create LAMMPS instance
+         L = IPyLammps(cmdargs=args)
+         # get and print numerical version code
+         print("LAMMPS Version: ", L.version())
+         # explicitly close and delete LAMMPS instance (optional)
+         L.close()
+
+      You can also initialize IPyLammps on top of an existing :py:class:`lammps` or :py:class:`PyLammps` object:
+
+      .. code-block:: Python
+
+         from lammps import lammps, IPyLammps
+         ...
+         # create LAMMPS instance
+         lmp = lammps(cmdargs=args)
+         # create PyLammps instance using previously created LAMMPS instance
+         L = PyLammps(ptr=lmp)
+
+      This is useful if you have to create the :py:class:`lammps <lammps.lammps>`
+      instance is a specific way, but want to take advantage of the
+      :py:class:`IPyLammps <lammps.IPyLammps>` interface.
+
+In all of the above cases, same as with the :ref:`C library API <lammps_c_api>`, this will use the
 ``MPI_COMM_WORLD`` communicator for the MPI library that LAMMPS was
 compiled with.  The :py:func:`lmp.close() <lammps.lammps.close>` call is
 optional since the LAMMPS class instance will also be deleted
@@ -73,39 +192,109 @@ destructor.
 Executing LAMMPS commands
 *************************
 
-Once an instance of the :py:class:`lammps <lammps.lammps>` class is
-created, there are multiple ways to "feed" it commands. In a way that is
-not very different from running a LAMMPS input script, except that
-Python has many more facilities for structured programming than the
-LAMMPS input script syntax.  Furthermore it is possible to "compute"
-what the next LAMMPS command should be. Same as in the equivalent `C
-library functions <pg_lib_execute>`, commands can be read from a file, a
-single string, a list of strings and a block of commands in a single
-multi-line string. They are processed under the same boundary conditions
-as the C library counterparts.  The example below demonstrates the use
-of :py:func:`lammps.file`, :py:func:`lammps.command`,
-:py:func:`lammps.commands_list`, and :py:func:`lammps.commands_string`:
+Once an instance of the :py:class:`lammps`, :py:class:`PyLammps`, or
+:py:class:`IPyLammps` class is created, there are multiple ways to "feed" it
+commands. In a way that is not very different from running a LAMMPS input
+script, except that Python has many more facilities for structured
+programming than the LAMMPS input script syntax. Furthermore it is possible
+to "compute" what the next LAMMPS command should be.
 
-.. code-block:: python
+.. tabs::
 
-   from lammps import lammps
+   .. tab:: lammps API
 
-   lmp = lammps()
-   # read commands from file 'in.melt'
-   lmp.file('in.melt')
-   # issue a single command
-   lmp.command('variable zpos index 1.0')
-   # create 10 groups with 10 atoms each
-   cmds = ["group g{} id {}:{}".format(i,10*i+1,10*(i+1)) for i in range(10)]
-   lmp.commands_list(cmds)
-   # run commands from a multi-line string
-   block = """
-   clear
-   region  box block 0 2 0 2 0 2
-   create_box 1 box
-   create_atoms 1 single 1.0 1.0 ${zpos}
-   """
-   lmp.commands_string(block)
+      Same as in the equivalent
+      :doc:`C library functions <pg_lib_execute>`, commands can be read from a file, a
+      single string, a list of strings and a block of commands in a single
+      multi-line string. They are processed under the same boundary conditions
+      as the C library counterparts.  The example below demonstrates the use
+      of :py:func:`lammps.file`, :py:func:`lammps.command`,
+      :py:func:`lammps.commands_list`, and :py:func:`lammps.commands_string`:
+
+      .. code-block:: python
+
+         from lammps import lammps
+         lmp = lammps()
+         # read commands from file 'in.melt'
+         lmp.file('in.melt')
+         # issue a single command
+         lmp.command('variable zpos index 1.0')
+         # create 10 groups with 10 atoms each
+         cmds = ["group g{} id {}:{}".format(i,10*i+1,10*(i+1)) for i in range(10)]
+         lmp.commands_list(cmds)
+         # run commands from a multi-line string
+         block = """
+         clear
+         region  box block 0 2 0 2 0 2
+         create_box 1 box
+         create_atoms 1 single 1.0 1.0 ${zpos}
+         """
+         lmp.commands_string(block)
+
+   .. tab:: PyLammps/IPyLammps API
+
+      Unlike the lammps API, the PyLammps/IPyLammps APIs allow running LAMMPS
+      commands by calling equivalent member functions.
+
+      For instance, the following LAMMPS command
+
+      .. code-block:: LAMMPS
+
+         region box block 0 10 0 5 -0.5 0.5
+
+      can be executed using the following Python code if *L* is a :py:class:`lammps` instance:
+
+      .. code-block:: Python
+
+         L.command("region box block 0 10 0 5 -0.5 0.5")
+
+      With the PyLammps interface, any LAMMPS command can be split up into arbitrary parts. 
+      These parts are then passed to a member function with the name of the command.
+      For the ``region`` command that means the :code:`region` method can be called.
+      The arguments of the command can be passed as one string, or
+      individually.
+
+      .. code-block:: Python
+
+         L.region("box block", 0, 10, 0, 5, -0.5, 0.5)
+
+      In this example all parameters except the first are Python floating-point literals. The
+      PyLammps interface takes the entire parameter list and transparently
+      merges it to a single command string.
+
+      The benefit of this approach is avoiding redundant command calls and easier
+      parameterization. In the original interface parameterization this needed to be done
+      manually by creating formatted strings.
+
+      .. code-block:: Python
+
+         L.command("region box block %f %f %f %f %f %f" % (xlo, xhi, ylo, yhi, zlo, zhi))
+
+      In contrast, methods of PyLammps accept parameters directly and will convert
+      them automatically to a final command string.
+
+      .. code-block:: Python
+
+         L.region("box block", xlo, xhi, ylo, yhi, zlo, zhi)
+
+      Using these facilities, the example shown for the lammps API can be rewritten as follows:
+
+      .. code-block:: python
+
+         from lammps import PyLammps
+         L = PyLammps()
+         # read commands from file 'in.melt'
+         L.file('in.melt')
+         # issue a single command
+         L.variable('zpos', 'index', 1.0)
+         # create 10 groups with 10 atoms each
+         for i in range(10):
+            L.group(f"g{i}", "id", f"{10*i+1}:{10*(i+1)}")
+
+         L.clear()
+         L.region("box block", 0, 2, 0, 2, 0, 2)
+         L.create_box(1, "box")
+         L.create_atoms(1, "single", 1.0, 1.0, "${zpos}")
 
 ----------
 
@@ -130,14 +319,34 @@ functions. Below is a detailed documentation of the API.
 The ``PyLammps`` class API
 **************************
 
+The :py:class:`PyLammps <lammps.PyLammps>` class is a wrapper that creates a
+simpler, more "Pythonic" interface to common LAMMPS functionality. LAMMPS
+data structures are exposed through objects and properties. This makes Python
+scripts shorter and more concise. See the :doc:`PyLammps Tutorial
+<Howto_pylammps>` for an introduction on how to use this interface.
+
 .. autoclass:: lammps.PyLammps
    :members:
 
+.. autoclass:: lammps.AtomList
+   :members:
+
+.. autoclass:: lammps.Atom
+   :members:
+
+.. autoclass:: lammps.Atom2D
+   :members:
+
 ----------
 
 The ``IPyLammps`` class API
 ***************************
 
+The :py:class:`IPyLammps <lammps.PyLammps>` class is an extension of
+:py:class:`PyLammps <lammps.PyLammps>`, adding additional functions to
+quickly display visualizations such as images and videos inside of IPython.
+See the :doc:`PyLammps Tutorial <Howto_pylammps>` for examples.
+
 .. autoclass:: lammps.IPyLammps
    :members:
 
@@ -150,14 +359,24 @@ The :py:mod:`lammps` module additionally contains several constants
 and the :py:class:`NeighList <lammps.NeighList>` class:
 
 .. _py_data_constants:
-.. py:data:: LAMMPS_INT, LAMMPS_DOUBLE, LAMMPS_BIGINT, LAMMPS_TAGINT, LAMMPS_STRING
+
+Data Types
+----------
+
+.. py:data:: LAMMPS_INT, LAMMPS_INT_2D, LAMMPS_DOUBLE, LAMMPS_DOUBLE_2D, LAMMPS_INT64, LAMMPS_INT64_2D, LAMMPS_STRING
    :type: int
 
    Constants in the :py:mod:`lammps` module to indicate how to
    cast data when the C library function returns a void pointer.
-   Used in :py:func:`lammps.extract_global`.
+   Used in :py:func:`lammps.extract_global` and :py:func:`lammps.extract_atom`.
+   See :cpp:enum:`_LMP_DATATYPE_CONST` for the equivalent constants in the
+   C library interface.
 
 .. _py_style_constants:
+
+Style Constants
+---------------
+
 .. py:data:: LMP_STYLE_GLOBAL, LMP_STYLE_ATOM, LMP_STYLE_LOCAL
    :type: int
 
@@ -167,6 +386,10 @@ and the :py:class:`NeighList <lammps.NeighList>` class:
    :py:func:`lammps.extract_compute` and :py:func:`lammps.extract_fix`.
 
 .. _py_type_constants:
+
+Type Constants
+--------------
+
 .. py:data:: LMP_TYPE_SCALAR, LMP_TYLE_VECTOR, LMP_TYPE_ARRAY, LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS
    :type: int
 
@@ -176,13 +399,36 @@ and the :py:class:`NeighList <lammps.NeighList>` class:
    :py:func:`lammps.extract_compute` and :py:func:`lammps.extract_fix`.
 
 .. _py_var_constants:
+
+Variable Style Constants
+------------------------
+
 .. py:data:: LMP_VAR_EQUAL, LMP_VAR_ATOM
    :type: int
 
    Constants in the :py:mod:`lammps` module to select what style of
    variable to query when calling :py:func:`lammps.extract_variable`.
 
+Classes representing internal objects
+-------------------------------------
+
 .. autoclass:: lammps.NeighList
    :members:
    :no-undoc-members:
 
+
+LAMMPS error handling in Python
+*******************************
+
+Compiling the shared library with :ref:`C++ exception support <exceptions>` provides a better error
+handling experience. Without exceptions the LAMMPS code will terminate the
+current Python process with an error message.  C++ exceptions allow capturing
+them on the C++ side and rethrowing them on the Python side. This way
+LAMMPS errors can be handled through the Python exception handling mechanism.
+
+.. warning::
+
+   Capturing a LAMMPS exception in Python can still mean that the
+   current LAMMPS process is in an illegal state and must be terminated. It is
+   advised to save your data and terminate the Python instance as quickly as
+   possible.
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 17e4ed4fd9..73041d8c02 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -2018,6 +2018,7 @@ Nakano
 nall
 namespace
 namespaces
+namedtuple
 nan
 NaN
 Nandor
diff --git a/examples/python/in.fix_python_invoke b/examples/python/in.fix_python_invoke
index f77cc15c90..ac435762c5 100644
--- a/examples/python/in.fix_python_invoke
+++ b/examples/python/in.fix_python_invoke
@@ -23,12 +23,12 @@ from lammps import lammps
 
 def end_of_step_callback(lmp):
   L = lammps(ptr=lmp)
-  t = L.extract_global("ntimestep", 0)
+  t = L.extract_global("ntimestep")
   print("### END OF STEP ###", t)
 
 def post_force_callback(lmp, v):
   L = lammps(ptr=lmp)
-  t = L.extract_global("ntimestep", 0)
+  t = L.extract_global("ntimestep")
   print("### POST_FORCE ###", t)
 """
 
diff --git a/examples/python/in.fix_python_invoke_neighlist b/examples/python/in.fix_python_invoke_neighlist
index 50f1d52c33..e5445227b1 100644
--- a/examples/python/in.fix_python_invoke_neighlist
+++ b/examples/python/in.fix_python_invoke_neighlist
@@ -35,14 +35,13 @@ def post_force_callback(lmp, v):
     #mylist = L.get_neighlist(0)
     mylist = L.find_pair_neighlist("lj/cut", request=0)
     print(pid_prefix, mylist)
-    nlocal = L.extract_global("nlocal", 0)
-    nghost = L.extract_global("nghost", 0)
-    ntypes = L.extract_global("ntypes", 0)
-    mass = L.numpy.extract_atom_darray("mass", ntypes+1)
-    atype = L.numpy.extract_atom_iarray("type", nlocal+nghost)
-    x = L.numpy.extract_atom_darray("x", nlocal+nghost, dim=3)
-    v = L.numpy.extract_atom_darray("v", nlocal+nghost, dim=3)
-    f = L.numpy.extract_atom_darray("f", nlocal+nghost, dim=3)
+    nlocal = L.extract_global("nlocal")
+    nghost = L.extract_global("nghost")
+    mass = L.numpy.extract_atom("mass")
+    atype = L.numpy.extract_atom("type", nelem=nlocal+nghost)
+    x = L.numpy.extract_atom("x", nelem=nlocal+nghost, dim=3)
+    v = L.numpy.extract_atom("v", nelem=nlocal+nghost, dim=3)
+    f = L.numpy.extract_atom("f", nelem=nlocal+nghost, dim=3)
 
     for iatom, numneigh, neighs in mylist:
       print(pid_prefix, "- {}".format(iatom), x[iatom], v[iatom], f[iatom], " : ",  numneigh, "Neighbors")
diff --git a/examples/python/py_nve.py b/examples/python/py_nve.py
index 79331528b1..9ff7c0978b 100644
--- a/examples/python/py_nve.py
+++ b/examples/python/py_nve.py
@@ -1,12 +1,12 @@
 from __future__ import print_function
-import lammps
+from lammps import lammps, LAMMPS_INT, LAMMPS_DOUBLE
 import ctypes
 import traceback
 import numpy as np
 
 class LAMMPSFix(object):
     def __init__(self, ptr, group_name="all"):
-        self.lmp = lammps.lammps(ptr=ptr)
+        self.lmp = lammps(ptr=ptr)
         self.group_name = group_name
 
 class LAMMPSFixMove(LAMMPSFix):
@@ -39,19 +39,18 @@ class NVE(LAMMPSFixMove):
         assert(self.group_name == "all")
 
     def init(self):
-        dt = self.lmp.extract_global("dt", 1)
-        ftm2v = self.lmp.extract_global("ftm2v", 1)
-        self.ntypes = self.lmp.extract_global("ntypes", 0)
+        dt = self.lmp.extract_global("dt")
+        ftm2v = self.lmp.extract_global("ftm2v")
+        self.ntypes = self.lmp.extract_global("ntypes")
         self.dtv = dt
         self.dtf = 0.5 * dt * ftm2v
 
     def initial_integrate(self, vflag):
-        nlocal = self.lmp.extract_global("nlocal", 0)
-        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
-        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
-        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
-        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
-        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        mass = self.lmp.numpy.extract_atom("mass")
+        atype = self.lmp.numpy.extract_atom("type")
+        x = self.lmp.numpy.extract_atom("x")
+        v = self.lmp.numpy.extract_atom("v")
+        f = self.lmp.numpy.extract_atom("f")
 
         for i in range(x.shape[0]):
             dtfm = self.dtf / mass[int(atype[i])]
@@ -59,11 +58,10 @@ class NVE(LAMMPSFixMove):
             x[i,:] += self.dtv * v[i,:]
 
     def final_integrate(self):
-        nlocal = self.lmp.extract_global("nlocal", 0)
-        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
-        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
-        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
-        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        mass = self.lmp.numpy.extract_atom("mass")
+        atype = self.lmp.numpy.extract_atom("type")
+        v = self.lmp.numpy.extract_atom("v")
+        f = self.lmp.numpy.extract_atom("f")
 
         for i in range(v.shape[0]):
             dtfm = self.dtf / mass[int(atype[i])]
@@ -77,19 +75,19 @@ class NVE_Opt(LAMMPSFixMove):
         assert(self.group_name == "all")
 
     def init(self):
-        dt = self.lmp.extract_global("dt", 1)
-        ftm2v = self.lmp.extract_global("ftm2v", 1)
-        self.ntypes = self.lmp.extract_global("ntypes", 0)
+        dt = self.lmp.extract_global("dt")
+        ftm2v = self.lmp.extract_global("ftm2v")
+        self.ntypes = self.lmp.extract_global("ntypes")
         self.dtv = dt
         self.dtf = 0.5 * dt * ftm2v
-        self.mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
+        self.mass = self.lmp.numpy.extract_atom("mass")
 
     def initial_integrate(self, vflag):        
-        nlocal = self.lmp.extract_global("nlocal", 0)
-        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
-        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
-        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
-        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        nlocal = self.lmp.extract_global("nlocal")
+        atype = self.lmp.numpy.extract_atom("type")
+        x = self.lmp.numpy.extract_atom("x")
+        v = self.lmp.numpy.extract_atom("v")
+        f = self.lmp.numpy.extract_atom("f")
         dtf = self.dtf
         dtv = self.dtv
         mass = self.mass
@@ -102,13 +100,12 @@ class NVE_Opt(LAMMPSFixMove):
             x[:,d] += dtv * v[:,d]
 
     def final_integrate(self):
-        nlocal = self.lmp.extract_global("nlocal", 0)
-        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
-        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
-        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
-        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        nlocal = self.lmp.extract_global("nlocal")
+        mass = self.lmp.numpy.extract_atom("mass")
+        atype = self.lmp.numpy.extract_atom("type")
+        v = self.lmp.numpy.extract_atom("v")
+        f = self.lmp.numpy.extract_atom("f")
         dtf = self.dtf
-        dtv = self.dtv
         mass = self.mass
 
         dtfm = dtf / np.take(mass, atype)
diff --git a/python/examples/demo.py b/python/examples/demo.py
index 5d2b62fee3..0bd2c8f054 100755
--- a/python/examples/demo.py
+++ b/python/examples/demo.py
@@ -16,7 +16,7 @@ if len(argv) != 1:
   print("Syntax: demo.py")
   sys.exit()
 
-from lammps import lammps
+from lammps import lammps, LAMMPS_INT, LMP_STYLE_GLOBAL, LMP_VAR_EQUAL, LMP_VAR_ATOM
 lmp = lammps()
 
 # test out various library functions after running in.demo
@@ -25,18 +25,18 @@ lmp.file("in.demo")
 
 print("\nPython output:")
 
-natoms = lmp.extract_global("natoms",0)
-mass = lmp.extract_atom("mass",2)
-x = lmp.extract_atom("x",3)
+natoms = lmp.extract_global("natoms")
+mass = lmp.extract_atom("mass")
+x = lmp.extract_atom("x")
 print("Natoms, mass, x[0][0] coord =",natoms,mass[1],x[0][0])
 
-temp = lmp.extract_compute("thermo_temp",0,0)
+temp = lmp.extract_compute("thermo_temp", LMP_STYLE_GLOBAL, LAMMPS_INT)
 print("Temperature from compute =",temp)
 
-eng = lmp.extract_variable("eng",None,0)
+eng = lmp.extract_variable("eng",None, LMP_VAR_EQUAL)
 print("Energy from equal-style variable =",eng)
 
-vy = lmp.extract_variable("vy","all",1)
+vy = lmp.extract_variable("vy","all", LMP_VAR_ATOM)
 print("Velocity component from atom-style variable =",vy[1])
 
 vol = lmp.get_thermo("vol")
diff --git a/python/examples/mc.py b/python/examples/mc.py
index fb2bfabab9..fe7f6838c8 100755
--- a/python/examples/mc.py
+++ b/python/examples/mc.py
@@ -27,7 +27,7 @@ if len(argv) != 2:
 
 infile = sys.argv[1]
 
-from lammps import lammps
+from lammps import lammps, LAMMPS_INT, LMP_STYLE_GLOBAL, LMP_VAR_EQUAL
 lmp = lammps()
 
 # run infile one line at a time
@@ -42,14 +42,14 @@ lmp.command("variable e equal pe")
 
 lmp.command("run 0")
 
-natoms = lmp.extract_global("natoms",0)
-emin = lmp.extract_compute("thermo_pe",0,0) / natoms
+natoms = lmp.extract_global("natoms")
+emin = lmp.extract_compute("thermo_pe",LMP_STYLE_GLOBAL,LAMMPS_INT) / natoms
 lmp.command("variable emin equal $e")
 
 # disorder the system
 # estart = initial energy
 
-x = lmp.extract_atom("x",3)
+x = lmp.extract_atom("x")
 
 for i in range(natoms):
   x[i][0] += deltaperturb * (2*random.random()-1)
@@ -58,10 +58,10 @@ for i in range(natoms):
 lmp.command("variable elast equal $e")
 lmp.command("thermo_style custom step v_emin v_elast pe")
 lmp.command("run 0")
-x = lmp.extract_atom("x",3)
+x = lmp.extract_atom("x")
 lmp.command("variable elast equal $e")
   
-estart = lmp.extract_compute("thermo_pe",0,0) / natoms
+estart = lmp.extract_compute("thermo_pe", LMP_STYLE_GLOBAL, LAMMPS_INT) / natoms
 
 # loop over Monte Carlo moves
 # extract x after every run, in case reneighboring changed ptr in LAMMPS
@@ -78,8 +78,8 @@ for i in range(nloop):
   x[iatom][1] += deltamove * (2*random.random()-1)
 
   lmp.command("run 1 pre no post no")
-  x = lmp.extract_atom("x",3)
-  e = lmp.extract_compute("thermo_pe",0,0) / natoms
+  x = lmp.extract_atom("x")
+  e = lmp.extract_compute("thermo_pe", LMP_STYLE_GLOBAL, LAMMPS_INT) / natoms
 
   if e <= elast:
     elast = e
@@ -96,10 +96,10 @@ for i in range(nloop):
 # final energy and stats
 
 lmp.command("variable nbuild equal nbuild")
-nbuild = lmp.extract_variable("nbuild",None,0)
+nbuild = lmp.extract_variable("nbuild", None, LMP_VAR_EQUAL)
 
 lmp.command("run 0")
-estop = lmp.extract_compute("thermo_pe",0,0) / natoms
+estop = lmp.extract_compute("thermo_pe", LMP_STYLE_GLOBAL, LAMMPS_INT) / natoms
 
 print("MC stats:")
 print("  starting energy =",estart)
diff --git a/python/examples/simple.py b/python/examples/simple.py
index 632dde9ed3..9f85d8b4df 100755
--- a/python/examples/simple.py
+++ b/python/examples/simple.py
@@ -61,7 +61,7 @@ lmp.command("run 1");
 
 # extract force on single atom two different ways
 
-f = lmp.extract_atom("f",3)
+f = lmp.extract_atom("f")
 print("Force on 1 atom via extract_atom: ",f[0][0])
 
 fx = lmp.extract_variable("fx","all",1)
diff --git a/python/examples/split.py b/python/examples/split.py
index 77e605092e..bd2896c004 100755
--- a/python/examples/split.py
+++ b/python/examples/split.py
@@ -57,7 +57,7 @@ if color == 0:
   lmp.scatter_atoms("x",1,3,x)
   lmp.command("run 1");
 
-  f = lmp.extract_atom("f",3)
+  f = lmp.extract_atom("f")
   print("Force on 1 atom via extract_atom: ",f[0][0])
 
   fx = lmp.extract_variable("fx","all",1)
diff --git a/python/lammps.py b/python/lammps.py
index 6e910b49bf..832a7f3d4d 100644
--- a/python/lammps.py
+++ b/python/lammps.py
@@ -19,6 +19,7 @@ from __future__ import print_function
 # imports for simple LAMMPS python wrapper module "lammps"
 
 import sys,traceback,types
+import warnings
 from ctypes import *
 from os.path import dirname,abspath,join
 from inspect import getsourcefile
@@ -33,12 +34,13 @@ import sys
 
 # various symbolic constants to be used
 # in certain calls to select data formats
+LAMMPS_AUTODETECT = None
 LAMMPS_INT    = 0
-LAMMPS_INT2D  = 1
+LAMMPS_INT_2D  = 1
 LAMMPS_DOUBLE = 2
-LAMMPS_DOUBLE2D = 3
-LAMMPS_BIGINT = 4
-LAMMPS_TAGINT = 5
+LAMMPS_DOUBLE_2D = 3
+LAMMPS_INT64 = 4
+LAMMPS_INT64_2D = 5
 LAMMPS_STRING = 6
 
 # these must be kept in sync with the enums in library.h
@@ -313,6 +315,8 @@ class lammps(object):
     self.lib.lammps_get_last_error_message.restype = c_int
 
     self.lib.lammps_extract_global.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_extract_global_datatype.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_extract_global_datatype.restype = c_int
     self.lib.lammps_extract_compute.argtypes = [c_void_p, c_char_p, c_int, c_int]
 
     self.lib.lammps_get_thermo.argtypes = [c_void_p, c_char_p]
@@ -337,6 +341,8 @@ class lammps(object):
     self.lib.lammps_decode_image_flags.argtypes = [self.c_imageint, POINTER(c_int*3)]
 
     self.lib.lammps_extract_atom.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_extract_atom_datatype.argtypes = [c_void_p, c_char_p]
+    self.lib.lammps_extract_atom_datatype.restype = c_int
 
     self.lib.lammps_extract_fix.argtypes = [c_void_p, c_char_p, c_int, c_int, c_int, c_int]
 
@@ -473,7 +479,73 @@ class lammps(object):
             return np.int64
           return np.intc
 
+        def extract_atom(self, name, dtype=LAMMPS_AUTODETECT, nelem=LAMMPS_AUTODETECT, dim=LAMMPS_AUTODETECT):
+          """Retrieve per-atom properties from LAMMPS as NumPy arrays
+
+          This is a wrapper around the :cpp:func:`lammps_extract_atom`
+          function of the C-library interface. Its documentation includes a
+          list of the supported keywords and their data types.
+          Since Python needs to know the data type to be able to interpret
+          the result, by default, this function will try to auto-detect the data 
+          type by asking the library. You can also force a specific data type.
+          For that purpose the :py:mod:`lammps` module contains the constants
+          ``LAMMPS_INT``, ``LAMMPS_INT_2D``, ``LAMMPS_DOUBLE``,
+          ``LAMMPS_DOUBLE_2D``, ``LAMMPS_INT64``, ``LAMMPS_INT64_2D``, and
+          ``LAMMPS_STRING``.
+          This function returns ``None`` if either the keyword is not
+          recognized, or an invalid data type constant is used.
+
+          .. note::
+
+             While the returned arrays of per-atom data are dimensioned
+             for the range [0:nmax] - as is the underlying storage -
+             the data is usually only valid for the range of [0:nlocal],
+             unless the property of interest is also updated for ghost
+             atoms.  In some cases, this depends on a LAMMPS setting, see
+             for example :doc:`comm_modify vel yes <comm_modify>`.
+
+          :param name: name of the property
+          :type name:  string
+          :param dtype: type of the returned data (see :ref:`py_data_constants`)
+          :type dtype:  int, optional
+          :param nelem: number of elements in array
+          :type nelem:  int, optional
+          :param dim: dimension of each element
+          :type dim:  int, optional
+          :return: requested data as NumPy array with direct access to C data
+          :rtype: numpy.array
+          """
+          if dtype == LAMMPS_AUTODETECT:
+            dtype = self.lmp.extract_atom_datatype(name)
+
+          if nelem == LAMMPS_AUTODETECT:
+            if name == "mass":
+              nelem = self.lmp.extract_global("ntypes") + 1
+            else:
+              nelem = self.lmp.extract_global("nlocal")
+          if dim == LAMMPS_AUTODETECT:
+            if dtype in (LAMMPS_INT_2D, LAMMPS_DOUBLE_2D, LAMMPS_INT64_2D):
+              # TODO add other fields
+              if name in ("x", "v", "f", "angmom", "torque", "csforce", "vforce"):
+                dim = 3
+              else:
+                dim = 2
+            else:
+              dim = 1
+
+          raw_ptr = self.lmp.extract_atom(name, dtype)
+
+          if dtype in (LAMMPS_DOUBLE, LAMMPS_DOUBLE_2D):
+            return self.darray(raw_ptr, nelem, dim)
+          elif dtype in (LAMMPS_INT, LAMMPS_INT_2D):
+            return self.iarray(c_int32, raw_ptr, nelem, dim)
+          elif dtype in (LAMMPS_INT64, LAMMPS_INT64_2D):
+            return self.iarray(c_int64, raw_ptr, nelem, dim)
+          return raw_ptr
+
         def extract_atom_iarray(self, name, nelem, dim=1):
+          warnings.warn("deprecated, use extract_atom instead", DeprecationWarning)
+
           if name in ['id', 'molecule']:
             c_int_type = self.lmp.c_tagint
           elif name in ['image']:
@@ -484,15 +556,17 @@ class lammps(object):
           if dim == 1:
             raw_ptr = self.lmp.extract_atom(name, LAMMPS_INT)
           else:
-            raw_ptr = self.lmp.extract_atom(name, LAMMPS_INT2D)
+            raw_ptr = self.lmp.extract_atom(name, LAMMPS_INT_2D)
 
           return self.iarray(c_int_type, raw_ptr, nelem, dim)
 
         def extract_atom_darray(self, name, nelem, dim=1):
+          warnings.warn("deprecated, use extract_atom instead", DeprecationWarning)
+
           if dim == 1:
             raw_ptr = self.lmp.extract_atom(name, LAMMPS_DOUBLE)
           else:
-            raw_ptr = self.lmp.extract_atom(name, LAMMPS_DOUBLE2D)
+            raw_ptr = self.lmp.extract_atom(name, LAMMPS_DOUBLE_2D)
 
           return self.darray(raw_ptr, nelem, dim)
 
@@ -705,9 +779,9 @@ class lammps(object):
     underlying :cpp:func:`lammps_get_natoms` function returning a double.
 
     :return: number of atoms
-    :rtype: float
+    :rtype: int
     """
-    return self.lib.lammps_get_natoms(self.lmp)
+    return int(self.lib.lammps_get_natoms(self.lmp))
 
   # -------------------------------------------------------------------------
 
@@ -801,10 +875,35 @@ class lammps(object):
     else: return None
     return int(self.lib.lammps_extract_setting(self.lmp,name))
 
+  # -------------------------------------------------------------------------
+  # extract global info datatype
+
+  def extract_global_datatype(self, name):
+    """Retrieve global property datatype from LAMMPS
+
+    This is a wrapper around the :cpp:func:`lammps_extract_global_datatype`
+    function of the C-library interface. Its documentation includes a
+    list of the supported keywords.
+    This function returns ``None`` if the keyword is not
+    recognized. Otherwise it will return a positive integer value that
+    corresponds to one of the constants define in the :py:mod:`lammps` module:
+    ``LAMMPS_INT``, ``LAMMPS_INT_2D``, ``LAMMPS_DOUBLE``, ``LAMMPS_DOUBLE_2D``,
+    ``LAMMPS_INT64``, ``LAMMPS_INT64_2D``, and ``LAMMPS_STRING``. These values
+    are equivalent to the ones defined in :cpp:enum:`_LMP_DATATYPE_CONST`.
+
+    :param name: name of the property
+    :type name:  string
+    :return: datatype of global property
+    :rtype: int
+    """
+    if name: name = name.encode()
+    else: return None
+    return self.lib.lammps_extract_global_datatype(self.lmp, name)
+
   # -------------------------------------------------------------------------
   # extract global info
 
-  def extract_global(self, name, type):
+  def extract_global(self, name, dtype=LAMMPS_AUTODETECT):
     """Query LAMMPS about global settings of different types.
 
     This is a wrapper around the :cpp:func:`lammps_extract_global`
@@ -814,55 +913,87 @@ class lammps(object):
     of values.  The :cpp:func:`lammps_extract_global` documentation
     includes a list of the supported keywords and their data types.
     Since Python needs to know the data type to be able to interpret
-    the result, the type has to be provided as an argument.  For
-    that purpose the :py:mod:`lammps` module contains the constants
-    ``LAMMPS_INT``, ``LAMMPS_DOUBLE``, ``LAMMPS_BIGINT``,
-    ``LAMMPS_TAGINT``, and ``LAMMPS_STRING``.
-    This function returns ``None`` if either the keyword is not
-    recognized, or an invalid data type constant is used.
+    the result, by default, this function will try to auto-detect the data type
+    by asking the library. You can also force a specific data type.  For that
+    purpose the :py:mod:`lammps` module contains the constants ``LAMMPS_INT``,
+    ``LAMMPS_DOUBLE``, ``LAMMPS_INT64``, and ``LAMMPS_STRING``. These values
+    are equivalent to the ones defined in :cpp:enum:`_LMP_DATATYPE_CONST`.
+    This function returns ``None`` if either the keyword is not recognized,
+    or an invalid data type constant is used.
 
-    :param name: name of the setting
+    :param name: name of the property
     :type name:  string
-    :param type: type of the returned data
-    :type type:  int
-    :return: value of the setting
-    :rtype: integer or double or string or None
+    :param dtype: data type of the returned data (see :ref:`py_data_constants`)
+    :type dtype:  int, optional
+    :return: value of the property or None
+    :rtype: int, float, or NoneType
+    """
+    if dtype == LAMMPS_AUTODETECT:
+      dtype = self.extract_global_datatype(name)
+
+    if name: name = name.encode()
+    else: return None
+
+    if dtype == LAMMPS_INT:
+      self.lib.lammps_extract_global.restype = POINTER(c_int32)
+      target_type = int
+    elif dtype == LAMMPS_INT64:
+      self.lib.lammps_extract_global.restype = POINTER(c_int64)
+      target_type = int
+    elif dtype == LAMMPS_DOUBLE:
+      self.lib.lammps_extract_global.restype = POINTER(c_double)
+      target_type = float
+    elif dtype == LAMMPS_STRING:
+      self.lib.lammps_extract_global.restype = c_char_p
+      target_type = lambda x: str(x, 'ascii')
+
+    ptr = self.lib.lammps_extract_global(self.lmp, name)
+    if ptr:
+      return target_type(ptr[0])
+    return None
+
+
+  # -------------------------------------------------------------------------
+  # extract per-atom info datatype
+
+  def extract_atom_datatype(self, name):
+    """Retrieve per-atom property datatype from LAMMPS
+
+    This is a wrapper around the :cpp:func:`lammps_extract_atom_datatype`
+    function of the C-library interface. Its documentation includes a
+    list of the supported keywords.
+    This function returns ``None`` if the keyword is not
+    recognized. Otherwise it will return an integer value that
+    corresponds to one of the constants define in the :py:mod:`lammps` module:
+    ``LAMMPS_INT``, ``LAMMPS_INT_2D``, ``LAMMPS_DOUBLE``, ``LAMMPS_DOUBLE_2D``,
+    ``LAMMPS_INT64``, ``LAMMPS_INT64_2D``, and ``LAMMPS_STRING``. These values
+    are equivalent to the ones defined in :cpp:enum:`_LMP_DATATYPE_CONST`.
+
+    :param name: name of the property
+    :type name:  string
+    :return: data type of per-atom property (see :ref:`py_data_constants`)
+    :rtype: int
     """
     if name: name = name.encode()
     else: return None
-    if type == LAMMPS_INT:
-      self.lib.lammps_extract_global.restype = POINTER(c_int)
-    elif type == LAMMPS_DOUBLE:
-      self.lib.lammps_extract_global.restype = POINTER(c_double)
-    elif type == LAMMPS_BIGINT:
-      self.lib.lammps_extract_global.restype = POINTER(self.c_bigint)
-    elif type == LAMMPS_TAGINT:
-      self.lib.lammps_extract_global.restype = POINTER(self.c_tagint)
-    elif type == LAMMPS_STRING:
-      self.lib.lammps_extract_global.restype = c_char_p
-      ptr = self.lib.lammps_extract_global(self.lmp,name)
-      return str(ptr,'ascii')
-    else: return None
-    ptr = self.lib.lammps_extract_global(self.lmp,name)
-    if ptr: return ptr[0]
-    else: return None
+    return self.lib.lammps_extract_atom_datatype(self.lmp, name)
 
   # -------------------------------------------------------------------------
   # extract per-atom info
-  # NOTE: need to insure are converting to/from correct Python type
-  #   e.g. for Python list or NumPy or ctypes
 
-  def extract_atom(self,name,type):
+  def extract_atom(self, name, dtype=LAMMPS_AUTODETECT):
     """Retrieve per-atom properties from LAMMPS
 
     This is a wrapper around the :cpp:func:`lammps_extract_atom`
     function of the C-library interface. Its documentation includes a
     list of the supported keywords and their data types.
     Since Python needs to know the data type to be able to interpret
-    the result, the type has to be provided as an argument.  For
+    the result, by default, this function will try to auto-detect the data type
+    by asking the library. You can also force a specific data type.  For
     that purpose the :py:mod:`lammps` module contains the constants
-    ``LAMMPS_INT``, ``LAMMPS_INT2D``, ``LAMMPS_DOUBLE``,
-    and ``LAMMPS_DOUBLE2D``.
+    ``LAMMPS_INT``, ``LAMMPS_INT_2D``, ``LAMMPS_DOUBLE``, ``LAMMPS_DOUBLE_2D``,
+    ``LAMMPS_INT64``, ``LAMMPS_INT64_2D``, and ``LAMMPS_STRING``. These values
+    are equivalent to the ones defined in :cpp:enum:`_LMP_DATATYPE_CONST`.
     This function returns ``None`` if either the keyword is not
     recognized, or an invalid data type constant is used.
 
@@ -875,27 +1006,36 @@ class lammps(object):
        atoms.  In some cases, this depends on a LAMMPS setting, see
        for example :doc:`comm_modify vel yes <comm_modify>`.
 
-    :param name: name of the setting
+    :param name: name of the property
     :type name:  string
-    :param type: type of the returned data
-    :type type:  int
-    :return: requested data
-    :rtype: pointer to integer or double or None
+    :param dtype: data type of the returned data (see :ref:`py_data_constants`)
+    :type dtype:  int, optional
+    :return: requested data or ``None``
+    :rtype: ctypes.POINTER(ctypes.c_int32), ctypes.POINTER(ctypes.POINTER(ctypes.c_int32)),
+            ctypes.POINTER(ctypes.c_int64), ctypes.POINTER(ctypes.POINTER(ctypes.c_int64)),
+            ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
+            or NoneType
     """
-    ntypes = int(self.extract_setting('ntypes'))
-    nmax   = int(self.extract_setting('nmax'))
+    if dtype == LAMMPS_AUTODETECT:
+      dtype = self.extract_atom_datatype(name)
+
     if name: name = name.encode()
     else: return None
-    if type == LAMMPS_INT:
-      self.lib.lammps_extract_atom.restype = POINTER(c_int)
-    elif type == LAMMPS_INT2D:
-      self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_int))
-    elif type == LAMMPS_DOUBLE:
+
+    if dtype == LAMMPS_INT:
+      self.lib.lammps_extract_atom.restype = POINTER(c_int32)
+    elif dtype == LAMMPS_INT_2D:
+      self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_int32))
+    elif dtype == LAMMPS_DOUBLE:
       self.lib.lammps_extract_atom.restype = POINTER(c_double)
-    elif type == LAMMPS_DOUBLE2D:
+    elif dtype == LAMMPS_DOUBLE_2D:
       self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_double))
+    elif dtype == LAMMPS_INT64:
+      self.lib.lammps_extract_atom.restype = POINTER(c_int64)
+    elif dtype == LAMMPS_INT64_2D:
+      self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_int64))
     else: return None
-    ptr = self.lib.lammps_extract_atom(self.lmp,name)
+    ptr = self.lib.lammps_extract_atom(self.lmp, name)
     if ptr: return ptr
     else:   return None
 
@@ -1140,7 +1280,7 @@ class lammps(object):
 
   def gather_atoms(self,name,type,count):
     if name: name = name.encode()
-    natoms = self.lib.lammps_get_natoms(self.lmp)
+    natoms = self.get_natoms()
     if type == 0:
       data = ((count*natoms)*c_int)()
       self.lib.lammps_gather_atoms(self.lmp,name,type,count,data)
@@ -1154,7 +1294,7 @@ class lammps(object):
 
   def gather_atoms_concat(self,name,type,count):
     if name: name = name.encode()
-    natoms = self.lib.lammps_get_natoms(self.lmp)
+    natoms = self.get_natoms()
     if type == 0:
       data = ((count*natoms)*c_int)()
       self.lib.lammps_gather_atoms_concat(self.lmp,name,type,count,data)
@@ -1206,7 +1346,7 @@ class lammps(object):
   #   e.g. for Python list or NumPy or ctypes
   def gather(self,name,type,count):
     if name: name = name.encode()
-    natoms = self.lib.lammps_get_natoms(self.lmp)
+    natoms = self.get_natoms()
     if type == 0:
       data = ((count*natoms)*c_int)()
       self.lib.lammps_gather(self.lmp,name,type,count,data)
@@ -1218,7 +1358,7 @@ class lammps(object):
 
   def gather_concat(self,name,type,count):
     if name: name = name.encode()
-    natoms = self.lib.lammps_get_natoms(self.lmp)
+    natoms = self.get_natoms()
     if type == 0:
       data = ((count*natoms)*c_int)()
       self.lib.lammps_gather_concat(self.lmp,name,type,count,data)
@@ -1525,8 +1665,8 @@ class lammps(object):
   def available_styles(self, category):
     """Returns a list of styles available for a given category
 
-    This is a wrapper around the functions :cpp:func:`lammps_style_count`
-    and :cpp:func`lammps_style_name` of the library interface.
+    This is a wrapper around the functions :cpp:func:`lammps_style_count()`
+    and :cpp:func:`lammps_style_name()` of the library interface.
 
     :param category: name of category
     :type  category: string
@@ -1723,8 +1863,8 @@ class OutputCapture(object):
 # -------------------------------------------------------------------------
 
 class Variable(object):
-  def __init__(self, lammps_wrapper_instance, name, style, definition):
-    self.wrapper = lammps_wrapper_instance
+  def __init__(self, pylammps_instance, name, style, definition):
+    self._pylmp = pylammps_instance
     self.name = name
     self.style = style
     self.definition = definition.split()
@@ -1732,9 +1872,9 @@ class Variable(object):
   @property
   def value(self):
     if self.style == 'atom':
-      return list(self.wrapper.lmp.extract_variable(self.name, "all", 1))
+      return list(self._pylmp.lmp.extract_variable(self.name, "all", 1))
     else:
-      value = self.wrapper.lmp_print('"${%s}"' % self.name).strip()
+      value = self._pylmp.lmp_print('"${%s}"' % self.name).strip()
       try:
         return float(value)
       except ValueError:
@@ -1743,103 +1883,136 @@ class Variable(object):
 # -------------------------------------------------------------------------
 
 class AtomList(object):
-  def __init__(self, lammps_wrapper_instance):
-    self.lmp = lammps_wrapper_instance
-    self.natoms = self.lmp.system.natoms
-    self.dimensions = self.lmp.system.dimensions
+  """
+  A dynamic list of atoms that returns either an Atom or Atom2D instance for
+  each atom. Instances are only allocated when accessed.
+
+  :ivar natoms: total number of atoms
+  :ivar dimensions: number of dimensions in system
+  """
+  def __init__(self, pylammps_instance):
+    self._pylmp = pylammps_instance
+    self.natoms = self._pylmp.system.natoms
+    self.dimensions = self._pylmp.system.dimensions
+    self._loaded = {}
 
   def __getitem__(self, index):
-    if self.dimensions == 2:
-        return Atom2D(self.lmp, index + 1)
-    return Atom(self.lmp, index + 1)
+    """
+    Return Atom with given local index
+
+    :param index: Local index of atom
+    :type index: int
+    :rtype: Atom or Atom2D
+    """
+    if index not in self._loaded:
+        if self.dimensions == 2:
+            atom = Atom2D(self._pylmp, index + 1)
+        else:
+            atom = Atom(self._pylmp, index + 1)
+        self._loaded[index] = atom
+    return self._loaded[index]
+
+  def __len__(self):
+    return self.natoms
+
 
 # -------------------------------------------------------------------------
 
 class Atom(object):
-  def __init__(self, lammps_wrapper_instance, index):
-    self.lmp = lammps_wrapper_instance
+  """
+  A wrapper class then represents a single atom inside of LAMMPS
+
+  It provides access to properties of the atom and allows you to change some of them.
+  """
+  def __init__(self, pylammps_instance, index):
+    self._pylmp = pylammps_instance
     self.index = index
 
   @property
   def id(self):
-    return int(self.lmp.eval("id[%d]" % self.index))
+    return int(self._pylmp.eval("id[%d]" % self.index))
 
   @property
   def type(self):
-    return int(self.lmp.eval("type[%d]" % self.index))
+    return int(self._pylmp.eval("type[%d]" % self.index))
 
   @property
   def mol(self):
-    return self.lmp.eval("mol[%d]" % self.index)
+    return self._pylmp.eval("mol[%d]" % self.index)
 
   @property
   def mass(self):
-    return self.lmp.eval("mass[%d]" % self.index)
+    return self._pylmp.eval("mass[%d]" % self.index)
 
   @property
   def position(self):
-    return (self.lmp.eval("x[%d]" % self.index),
-            self.lmp.eval("y[%d]" % self.index),
-            self.lmp.eval("z[%d]" % self.index))
+    return (self._pylmp.eval("x[%d]" % self.index),
+            self._pylmp.eval("y[%d]" % self.index),
+            self._pylmp.eval("z[%d]" % self.index))
 
   @position.setter
   def position(self, value):
-     self.lmp.set("atom", self.index, "x", value[0])
-     self.lmp.set("atom", self.index, "y", value[1])
-     self.lmp.set("atom", self.index, "z", value[2])
+     self._pylmp.set("atom", self.index, "x", value[0])
+     self._pylmp.set("atom", self.index, "y", value[1])
+     self._pylmp.set("atom", self.index, "z", value[2])
 
   @property
   def velocity(self):
-    return (self.lmp.eval("vx[%d]" % self.index),
-            self.lmp.eval("vy[%d]" % self.index),
-            self.lmp.eval("vz[%d]" % self.index))
+    return (self._pylmp.eval("vx[%d]" % self.index),
+            self._pylmp.eval("vy[%d]" % self.index),
+            self._pylmp.eval("vz[%d]" % self.index))
 
   @velocity.setter
   def velocity(self, value):
-     self.lmp.set("atom", self.index, "vx", value[0])
-     self.lmp.set("atom", self.index, "vy", value[1])
-     self.lmp.set("atom", self.index, "vz", value[2])
+     self._pylmp.set("atom", self.index, "vx", value[0])
+     self._pylmp.set("atom", self.index, "vy", value[1])
+     self._pylmp.set("atom", self.index, "vz", value[2])
 
   @property
   def force(self):
-    return (self.lmp.eval("fx[%d]" % self.index),
-            self.lmp.eval("fy[%d]" % self.index),
-            self.lmp.eval("fz[%d]" % self.index))
+    return (self._pylmp.eval("fx[%d]" % self.index),
+            self._pylmp.eval("fy[%d]" % self.index),
+            self._pylmp.eval("fz[%d]" % self.index))
 
   @property
   def charge(self):
-    return self.lmp.eval("q[%d]" % self.index)
+    return self._pylmp.eval("q[%d]" % self.index)
 
 # -------------------------------------------------------------------------
 
 class Atom2D(Atom):
-  def __init__(self, lammps_wrapper_instance, index):
-    super(Atom2D, self).__init__(lammps_wrapper_instance, index)
+  """
+  A wrapper class then represents a single 2D atom inside of LAMMPS
+
+  It provides access to properties of the atom and allows you to change some of them.
+  """
+  def __init__(self, pylammps_instance, index):
+    super(Atom2D, self).__init__(pylammps_instance, index)
 
   @property
   def position(self):
-    return (self.lmp.eval("x[%d]" % self.index),
-            self.lmp.eval("y[%d]" % self.index))
+    return (self._pylmp.eval("x[%d]" % self.index),
+            self._pylmp.eval("y[%d]" % self.index))
 
   @position.setter
   def position(self, value):
-     self.lmp.set("atom", self.index, "x", value[0])
-     self.lmp.set("atom", self.index, "y", value[1])
+     self._pylmp.set("atom", self.index, "x", value[0])
+     self._pylmp.set("atom", self.index, "y", value[1])
 
   @property
   def velocity(self):
-    return (self.lmp.eval("vx[%d]" % self.index),
-            self.lmp.eval("vy[%d]" % self.index))
+    return (self._pylmp.eval("vx[%d]" % self.index),
+            self._pylmp.eval("vy[%d]" % self.index))
 
   @velocity.setter
   def velocity(self, value):
-     self.lmp.set("atom", self.index, "vx", value[0])
-     self.lmp.set("atom", self.index, "vy", value[1])
+     self._pylmp.set("atom", self.index, "vx", value[0])
+     self._pylmp.set("atom", self.index, "vy", value[1])
 
   @property
   def force(self):
-    return (self.lmp.eval("fx[%d]" % self.index),
-            self.lmp.eval("fy[%d]" % self.index))
+    return (self._pylmp.eval("fx[%d]" % self.index),
+            self._pylmp.eval("fy[%d]" % self.index))
 
 # -------------------------------------------------------------------------
 
@@ -1919,11 +2092,41 @@ def get_thermo_data(output):
 
 class PyLammps(object):
   """
-  More Python-like wrapper for LAMMPS (e.g., for IPython)
-  See examples/ipython for usage
+  This is a Python wrapper class around the lower-level
+  :py:class:`lammps` class, exposing a more Python-like,
+  object-oriented interface for prototyping system inside of IPython and
+  Jupyter notebooks.
+
+  It either creates its own instance of :py:class:`lammps` or can be
+  initialized with an existing instance. The arguments are the same of the
+  lower-level interface. The original interface can still be accessed via
+  :py:attr:`PyLammps.lmp`.
+
+  :param name: "machine" name of the shared LAMMPS library ("mpi" loads ``liblammps_mpi.so``, "" loads ``liblammps.so``)
+  :type  name: string
+  :param cmdargs: list of command line arguments to be passed to the :cpp:func:`lammps_open` function.  The executable name is automatically added.
+  :type  cmdargs: list
+  :param ptr: pointer to a LAMMPS C++ class instance when called from an embedded Python interpreter.  None means load symbols from shared library.
+  :type  ptr: pointer
+  :param comm: MPI communicator (as provided by `mpi4py <mpi4py_docs_>`_). ``None`` means use ``MPI_COMM_WORLD`` implicitly.
+  :type  comm: MPI_Comm
+
+  :ivar lmp:  instance of original LAMMPS Python interface
+  :vartype lmp: :py:class:`lammps`
+
+  :ivar runs:  list of completed runs, each storing the thermo output
+  :vartype run: list
   """
 
-  def __init__(self,name="",cmdargs=None,ptr=None,comm=None):
+  def __init__(self, name="", cmdargs=None, ptr=None, comm=None):
+    self.has_echo = False
+
+    if cmdargs:
+      if '-echo' in cmdargs:
+        idx = cmdargs.index('-echo')
+        # ensures that echo line is ignored during output capture
+        self.has_echo = idx+1 < len(cmdargs) and cmdargs[idx+1] in ('screen', 'both')
+
     if ptr:
       if isinstance(ptr,PyLammps):
         self.lmp = ptr.lmp
@@ -1942,26 +2145,65 @@ class PyLammps(object):
     self.lmp = None
 
   def close(self):
+    """Explicitly delete a LAMMPS instance
+
+    This is a wrapper around the :py:meth:`lammps.close` of the Python interface.
+    """
     if self.lmp: self.lmp.close()
     self.lmp = None
 
   def version(self):
+    """Return a numerical representation of the LAMMPS version in use.
+
+    This is a wrapper around the :py:meth:`lammps.version` function of the Python interface.
+
+    :return: version number
+    :rtype:  int
+    """
     return self.lmp.version()
 
-  def file(self,file):
+  def file(self, file):
+    """Read LAMMPS commands from a file.
+
+    This is a wrapper around the :py:meth:`lammps.file` function of the Python interface.
+
+    :param path: Name of the file/path with LAMMPS commands
+    :type path:  string
+    """
     self.lmp.file(file)
 
-  def write_script(self,filename):
-    """ Write LAMMPS script file containing all commands executed up until now """
-    with open(filename, "w") as f:
-      for cmd in self._cmd_history:
-        f.write("%s\n" % cmd)
+  def write_script(self, filepath):
+    """
+    Write LAMMPS script file containing all commands executed up until now
 
-  def command(self,cmd):
+    :param filepath: path to script file that should be written
+    :type filepath: string
+    """
+    with open(filepath, "w") as f:
+      for cmd in self._cmd_history:
+        print(cmd, file=f)
+
+  def command(self, cmd):
+    """
+    Execute LAMMPS command
+
+    All commands executed will be stored in a command history which can be
+    written to a file using :py:meth:`PyLammps.write_script()`
+
+    :param cmd: command string that should be executed
+    :type: cmd: string
+    """
     self.lmp.command(cmd)
     self._cmd_history.append(cmd)
 
   def run(self, *args, **kwargs):
+    """
+    Execute LAMMPS run command with given arguments
+
+    All thermo output during the run is captured and saved as new entry in
+    :py:attr:`PyLammps.runs`. The latest run can be retrieved by
+    :py:attr:`PyLammps.last_run`.
+    """
     output = self.__getattr__('run')(*args, **kwargs)
 
     if(self.has_mpi4py):
@@ -1972,48 +2214,102 @@ class PyLammps(object):
 
   @property
   def last_run(self):
+    """
+    Return data produced of last completed run command
+
+    :getter: Returns an object containing information about the last run command
+    :type: dict
+    """
     if len(self.runs) > 0:
         return self.runs[-1]
     return None
 
   @property
   def atoms(self):
+    """
+    All atoms of this LAMMPS instance
+
+    :getter: Returns a list of atoms currently in the system
+    :type: AtomList
+    """
     return AtomList(self)
 
   @property
   def system(self):
+    """
+    The system state of this LAMMPS instance
+
+    :getter: Returns an object with properties storing the current system state
+    :type: namedtuple
+    """
     output = self.info("system")
     d = self._parse_info_system(output)
     return namedtuple('System', d.keys())(*d.values())
 
   @property
   def communication(self):
+    """
+    The communication state of this LAMMPS instance
+
+    :getter: Returns an object with properties storing the current communication state
+    :type: namedtuple
+    """
     output = self.info("communication")
     d = self._parse_info_communication(output)
     return namedtuple('Communication', d.keys())(*d.values())
 
   @property
   def computes(self):
+    """
+    The list of active computes of this LAMMPS instance
+
+    :getter: Returns a list of computes that are currently active in this LAMMPS instance
+    :type: list
+    """
     output = self.info("computes")
     return self._parse_element_list(output)
 
   @property
   def dumps(self):
+    """
+    The list of active dumps of this LAMMPS instance
+
+    :getter: Returns a list of dumps that are currently active in this LAMMPS instance
+    :type: list
+    """
     output = self.info("dumps")
     return self._parse_element_list(output)
 
   @property
   def fixes(self):
+    """
+    The list of active fixes of this LAMMPS instance
+
+    :getter: Returns a list of fixes that are currently active in this LAMMPS instance
+    :type: list
+    """
     output = self.info("fixes")
     return self._parse_element_list(output)
 
   @property
   def groups(self):
+    """
+    The list of active atom groups of this LAMMPS instance
+
+    :getter: Returns a list of atom groups that are currently active in this LAMMPS instance
+    :type: list
+    """
     output = self.info("groups")
     return self._parse_groups(output)
 
   @property
   def variables(self):
+    """
+    Returns a dictionary of all variables defined in the current LAMMPS instance
+
+    :getter: Returns a dictionary of all variables that are defined in this LAMMPS instance
+    :type: dict
+    """
     output = self.info("variables")
     vars = {}
     for v in self._parse_element_list(output):
@@ -2021,6 +2317,15 @@ class PyLammps(object):
     return vars
 
   def eval(self, expr):
+    """
+    Evaluate expression
+
+    :param expr: the expression string that should be evaluated inside of LAMMPS
+    :type expr: string
+
+    :return: the value of the evaluated expression
+    :rtype: float if numeric, string otherwise
+    """
     value = self.lmp_print('"$(%s)"' % expr).strip()
     try:
       return float(value)
@@ -2156,11 +2461,23 @@ class PyLammps(object):
     'variable', 'velocity', 'write_restart']
 
   def __getattr__(self, name):
+    """
+    This method is where the Python 'magic' happens. If a method is not
+    defined by the class PyLammps, it assumes it is a LAMMPS command. It takes
+    all the arguments, concatinates them to a single string, and executes it using
+    :py:meth:`lammps.PyLammps.command()`.
+
+    :param verbose: Print output of command
+    :type verbose:  bool
+    :return: line or list of lines of output, None if no output
+    :rtype: list or string
+    """
     def handler(*args, **kwargs):
       cmd_args = [name] + [str(x) for x in args]
 
       with OutputCapture() as capture:
-        self.command(' '.join(cmd_args))
+        cmd = ' '.join(cmd_args)
+        self.command(cmd)
         output = capture.output
 
       if 'verbose' in kwargs and kwargs['verbose']:
@@ -2168,6 +2485,9 @@ class PyLammps(object):
 
       lines = output.splitlines()
 
+      if self.has_echo:
+        lines = lines[1:]
+
       if len(lines) > 1:
         return lines
       elif len(lines) == 1:
@@ -2179,14 +2499,56 @@ class PyLammps(object):
 
 class IPyLammps(PyLammps):
   """
-  IPython wrapper for LAMMPS which adds embedded graphics capabilities
+  IPython wrapper for LAMMPS which adds embedded graphics capabilities to PyLammmps interface
+
+  It either creates its own instance of :py:class:`lammps` or can be
+  initialized with an existing instance. The arguments are the same of the
+  lower-level interface. The original interface can still be accessed via
+  :py:attr:`PyLammps.lmp`.
+
+  :param name: "machine" name of the shared LAMMPS library ("mpi" loads ``liblammps_mpi.so``, "" loads ``liblammps.so``)
+  :type  name: string
+  :param cmdargs: list of command line arguments to be passed to the :cpp:func:`lammps_open` function.  The executable name is automatically added.
+  :type  cmdargs: list
+  :param ptr: pointer to a LAMMPS C++ class instance when called from an embedded Python interpreter.  None means load symbols from shared library.
+  :type  ptr: pointer
+  :param comm: MPI communicator (as provided by `mpi4py <mpi4py_docs_>`_). ``None`` means use ``MPI_COMM_WORLD`` implicitly.
+  :type  comm: MPI_Comm
   """
 
   def __init__(self,name="",cmdargs=None,ptr=None,comm=None):
     super(IPyLammps, self).__init__(name=name,cmdargs=cmdargs,ptr=ptr,comm=comm)
 
   def image(self, filename="snapshot.png", group="all", color="type", diameter="type",
-            size=None, view=None, center=None, up=None, zoom=1.0):
+            size=None, view=None, center=None, up=None, zoom=1.0, background_color="white"):
+    """ Generate image using write_dump command and display it
+
+    See :doc:`dump image <dump_image>` for more information.
+
+    :param filename: Name of the image file that should be generated. The extension determines whether it is PNG or JPEG
+    :type filename: string
+    :param group: the group of atoms write_image should use
+    :type group: string
+    :param color: name of property used to determine color
+    :type color: string
+    :param diameter: name of property used to determine atom diameter
+    :type diameter: string
+    :param size: dimensions of image
+    :type size: tuple (width, height)
+    :param view: view parameters
+    :type view: tuple (theta, phi)
+    :param center: center parameters
+    :type center: tuple (flag, center_x, center_y, center_z)
+    :param up: vector pointing to up direction
+    :type up: tuple (up_x, up_y, up_z)
+    :param zoom: zoom factor
+    :type zoom: float
+    :param background_color: background color of scene
+    :type background_color: string
+
+    :return: Image instance used to display image in notebook
+    :rtype: :py:class:`IPython.core.display.Image`
+    """
     cmd_args = [group, "image", filename, color, diameter]
 
     if size:
@@ -2215,12 +2577,22 @@ class IPyLammps(PyLammps):
     if zoom:
       cmd_args += ["zoom", zoom]
 
-    cmd_args.append("modify backcolor white")
+    cmd_args.append("modify backcolor " + background_color)
 
     self.write_dump(*cmd_args)
     from IPython.core.display import Image
-    return Image('snapshot.png')
+    return Image(filename)
 
   def video(self, filename):
+    """
+    Load video from file
+
+    Can be used to visualize videos from :doc:`dump movie <dump_image>`.
+
+    :param filename: Path to video file
+    :type filename: string
+    :return: HTML Video Tag used by notebook to embed a video
+    :rtype: :py:class:`IPython.display.HTML`
+    """
     from IPython.display import HTML
     return HTML("<video controls><source src=\"" + filename + "\"></video>")
diff --git a/src/atom.cpp b/src/atom.cpp
index 0322a55a96..fbf7a067a9 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -31,6 +31,8 @@
 #include "update.h"
 #include "variable.h"
 
+#include "library.h"
+
 #include <algorithm>
 #include <cstring>
 
@@ -2506,6 +2508,8 @@ length of the data area, and a short description.
    :cpp:func:`lammps_extract_atom`
 
 \endverbatim
+ *
+ * \sa extract_datatype
  *
  * \param  name  string with the keyword of the desired property.
                  Typically the name of the pointer variable returned
@@ -2515,6 +2519,8 @@ void *Atom::extract(const char *name)
 {
   // --------------------------------------------------------------------
   // 4th customization section: customize by adding new variable name
+  // please see the following function to set the type of the data
+  // so that programs can detect it dynamically at run time.
 
   /* NOTE: this array is only of length ntypes+1 */
   if (strcmp(name,"mass") == 0) return (void *) mass;
@@ -2583,6 +2589,89 @@ void *Atom::extract(const char *name)
   return nullptr;
 }
 
+
+/** Provide data type info about internal data of the Atom class
+ *
+\verbatim embed:rst
+
+.. versionadded:: 18Sep2020
+
+\endverbatim
+ *
+ * \sa extract
+ *
+ * \param  name  string with the keyword of the desired property.
+ * \return       data type constant for desired property or -1 */
+
+int Atom::extract_datatype(const char *name)
+{
+  // --------------------------------------------------------------------
+  // 5th customization section: customize by adding new variable name
+
+  if (strcmp(name,"mass") == 0) return LAMMPS_DOUBLE;
+
+  if (strcmp(name,"id") == 0) return LAMMPS_TAGINT;
+  if (strcmp(name,"type") == 0) return LAMMPS_INT;
+  if (strcmp(name,"mask") == 0) return LAMMPS_INT;
+  if (strcmp(name,"image") == 0) return LAMMPS_TAGINT;
+  if (strcmp(name,"x") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"v") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"f") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"molecule") == 0) return LAMMPS_TAGINT;
+  if (strcmp(name,"q") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"mu") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"omega") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"angmom") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"torque") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"radius") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"rmass") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"ellipsoid") == 0) return LAMMPS_INT;
+  if (strcmp(name,"line") == 0) return LAMMPS_INT;
+  if (strcmp(name,"tri") == 0) return LAMMPS_INT;
+  if (strcmp(name,"body") == 0) return LAMMPS_INT;
+
+  if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"x0") == 0) return LAMMPS_DOUBLE_2D;
+
+  if (strcmp(name,"spin") == 0) return LAMMPS_INT;
+  if (strcmp(name,"eradius") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"ervel") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"erforce") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"ervelforce") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"cs") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"csforce") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"vforce") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"etag") == 0) return LAMMPS_INT;
+
+  if (strcmp(name,"rho") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"drho") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"esph") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"desph") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"cv") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"vest") == 0) return LAMMPS_DOUBLE_2D;
+
+  // USER-MESONT package
+  if (strcmp(name,"length") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"buckling") == 0) return LAMMPS_INT;
+  if (strcmp(name,"bond_nt") == 0) return  LAMMPS_TAGINT_2D;
+
+  if (strcmp(name, "contact_radius") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name, "smd_data_9") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name, "smd_stress") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name, "eff_plastic_strain") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name, "eff_plastic_strain_rate") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name, "damage") == 0) return LAMMPS_DOUBLE;
+
+  if (strcmp(name,"dpdTheta") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"edpd_temp") == 0) return LAMMPS_DOUBLE;
+
+  // end of customization section
+  // --------------------------------------------------------------------
+
+  return -1;
+}
+
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
    call to avec tallies per-atom vectors
diff --git a/src/atom.h b/src/atom.h
index 10631d2435..bc69d3b27a 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -335,6 +335,7 @@ class Atom : protected Pointers {
   virtual void sync_modify(ExecutionSpace, unsigned int, unsigned int) {}
 
   void *extract(const char *);
+  int extract_datatype(const char *);
 
   inline int* get_map_array() {return map_array;};
   inline int get_map_size() {return map_tag_max+1;};
diff --git a/src/fmt/core.h b/src/fmt/core.h
index f9155000ec..6d87ab290a 100644
--- a/src/fmt/core.h
+++ b/src/fmt/core.h
@@ -18,7 +18,7 @@
 #include <vector>
 
 // The fmt library version in the form major * 10000 + minor * 100 + patch.
-#define FMT_VERSION 70002
+#define FMT_VERSION 70003
 
 #ifdef __clang__
 #  define FMT_CLANG_VERSION (__clang_major__ * 100 + __clang_minor__)
@@ -177,6 +177,12 @@
 #  endif
 #endif
 
+// LAMMPS customization
+// use 'v7_lmp' namespace instead of 'v7' so that our
+// bundled copy does not collide with linking other code
+// using system wide installations which may be using
+// a different version.
+
 #ifndef FMT_BEGIN_NAMESPACE
 #  if FMT_HAS_FEATURE(cxx_inline_namespaces) || FMT_GCC_VERSION >= 404 || \
       FMT_MSC_VER >= 1900
@@ -299,7 +305,7 @@ template <typename T> struct std_string_view {};
 
 #ifdef FMT_USE_INT128
 // Do nothing.
-#elif defined(__SIZEOF_INT128__) && !FMT_NVCC
+#elif defined(__SIZEOF_INT128__) && !FMT_NVCC && !(FMT_CLANG_VERSION && FMT_MSC_VER)
 #  define FMT_USE_INT128 1
 using int128_t = __int128_t;
 using uint128_t = __uint128_t;
@@ -489,6 +495,8 @@ constexpr basic_string_view<typename S::char_type> to_string_view(const S& s) {
   return s;
 }
 
+// LAMMPS customization using 'v7_lmp' instead of 'v7'
+
 namespace detail {
 void to_string_view(...);
 using fmt::v7_lmp::to_string_view;
@@ -1713,7 +1721,7 @@ template <typename Context> class basic_format_args {
   }
 
   template <typename Char> int get_id(basic_string_view<Char> name) const {
-    if (!has_named_args()) return {};
+    if (!has_named_args()) return -1;
     const auto& named_args =
         (is_packed() ? values_[-1] : args_[-1].value_).named_args;
     for (size_t i = 0; i < named_args.size; ++i) {
diff --git a/src/fmt/format.h b/src/fmt/format.h
index 5427042b4e..a4911b9fdb 100644
--- a/src/fmt/format.h
+++ b/src/fmt/format.h
@@ -69,6 +69,12 @@
 #  define FMT_NOINLINE
 #endif
 
+// LAMMPS customizations:
+// 1) Intel compilers on MacOS have __clang__ defined
+//    but fail to recognize [[clang::fallthrough]]
+// 2) Intel compilers on Linux identify as GCC compatible
+//    but fail to recognize [[gnu::fallthrough]]
+
 #if __cplusplus == 201103L || __cplusplus == 201402L
 #  if defined(__clang__) && !defined(__INTEL_COMPILER)
 #    define FMT_FALLTHROUGH [[clang::fallthrough]]
@@ -724,13 +730,18 @@ class FMT_API format_error : public std::runtime_error {
 
 namespace detail {
 
+template <typename T>
+using is_signed =
+    std::integral_constant<bool, std::numeric_limits<T>::is_signed ||
+                                     std::is_same<T, int128_t>::value>;
+
 // Returns true if value is negative, false otherwise.
 // Same as `value < 0` but doesn't produce warnings if T is an unsigned type.
-template <typename T, FMT_ENABLE_IF(std::numeric_limits<T>::is_signed)>
+template <typename T, FMT_ENABLE_IF(is_signed<T>::value)>
 FMT_CONSTEXPR bool is_negative(T value) {
   return value < 0;
 }
-template <typename T, FMT_ENABLE_IF(!std::numeric_limits<T>::is_signed)>
+template <typename T, FMT_ENABLE_IF(!is_signed<T>::value)>
 FMT_CONSTEXPR bool is_negative(T) {
   return false;
 }
@@ -745,9 +756,9 @@ FMT_CONSTEXPR bool is_supported_floating_point(T) {
 // Smallest of uint32_t, uint64_t, uint128_t that is large enough to
 // represent all values of T.
 template <typename T>
-using uint32_or_64_or_128_t = conditional_t<
-    num_bits<T>() <= 32, uint32_t,
-    conditional_t<num_bits<T>() <= 64, uint64_t, uint128_t>>;
+using uint32_or_64_or_128_t =
+    conditional_t<num_bits<T>() <= 32, uint32_t,
+                  conditional_t<num_bits<T>() <= 64, uint64_t, uint128_t>>;
 
 // Static data is placed in this class template for the header-only config.
 template <typename T = void> struct FMT_EXTERN_TEMPLATE_API basic_data {
@@ -1593,7 +1604,11 @@ template <typename OutputIt, typename Char, typename UInt> struct int_writer {
                               make_checked(p, s.size()));
     }
     if (prefix_size != 0) p[-1] = static_cast<Char>('-');
-    write(out, basic_string_view<Char>(buffer.data(), buffer.size()), specs);
+    using iterator = remove_reference_t<decltype(reserve(out, 0))>;
+    auto data = buffer.data();
+    out = write_padded<align::right>(out, specs, size, size, [=](iterator it) {
+      return copy_str<Char>(data, data + size, it);
+    });
   }
 
   void on_chr() { *out++ = static_cast<Char>(abs_value); }
diff --git a/src/library.cpp b/src/library.cpp
index 32c8728883..684cecc309 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -917,7 +917,7 @@ int lammps_extract_setting(void *handle, const char *keyword)
 This function returns a pointer to the location of some global property
 stored in one of the constituent classes of a LAMMPS instance.  The
 returned pointer is cast to ``void *`` and needs to be cast to a pointer
-of the type that the entity represents.  The pointers returned by this
+of the type that the entity represents. The pointers returned by this
 function are generally persistent; therefore it is not necessary to call
 the function again, unless a :doc:`clear` command is issued which wipes
 out and recreates the contents of the :cpp:class:`LAMMPS
@@ -1227,7 +1227,7 @@ of data type that the entity represents.
 A table with supported keywords is included in the documentation
 of the :cpp:func:`Atom::extract() <LAMMPS_NS::Atom::extract>` function.
 
-.. note::
+.. warning::
 
    The pointers returned by this function are generally not persistent
    since per-atom data may be re-distributed, re-allocated, and
@@ -1248,6 +1248,117 @@ void *lammps_extract_atom(void *handle, const char *name)
 
 /* ---------------------------------------------------------------------- */
 
+/** Get data type of internal global LAMMPS variables or arrays.
+ *
+\verbatim embed:rst
+
+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.
+
+.. versionadded:: 18Sep2020
+
+\endverbatim
+ *
+ * \param  handle   pointer to a previously created LAMMPS instance
+ * \param  name     string with the name of the extracted property
+ * \return          integer constant encoding the data type of the property
+ *                  or -1 if not found. */
+
+int lammps_extract_global_datatype(void *handle, const char *name)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+
+  if (strcmp(name,"units") == 0) return LAMMPS_STRING;
+  if (strcmp(name,"dt") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"ntimestep") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"boxlo") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxhi") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxxlo") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxxhi") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxylo") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxyhi") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxzlo") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"boxzhi") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"periodicity") == 0) return LAMMPS_INT;
+  if (strcmp(name,"triclinic") == 0) return LAMMPS_INT;
+  if (strcmp(name,"xy") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"xz") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"yz") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"natoms") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"nbonds") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"nangles") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"ndihedrals") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"nimpropers") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"nlocal") == 0) return LAMMPS_INT;
+  if (strcmp(name,"nghost") == 0) return LAMMPS_INT;
+  if (strcmp(name,"nmax") == 0) return LAMMPS_INT;
+  if (strcmp(name,"ntypes") == 0) return LAMMPS_INT;
+
+  if (strcmp(name,"q_flag") == 0) return LAMMPS_INT;
+
+  // update->atime can be referenced as a pointer
+  // thermo "timer" data cannot be, since it is computed on request
+  // lammps_get_thermo() can access all thermo keywords by value
+
+  if (strcmp(name,"atime") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"atimestep") == 0) return LAMMPS_BIGINT;
+
+  // global constants defined by units
+
+  if (strcmp(name,"boltz") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"hplanck") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"mvv2e") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"ftm2v") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"mv2d") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"nktv2p") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"qqr2e") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"qe2f") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"vxmu2f") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"xxt2kmu") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"dielectric") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"qqrd2e") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"e_mass") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"hhmrr2e") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"mvh2r") == 0) return LAMMPS_DOUBLE;
+
+  if (strcmp(name,"angstrom") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"femtosecond") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"qelectron") == 0) return LAMMPS_DOUBLE;
+
+  return -1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+/** Get data type of a LAMMPS per-atom property
+ *
+\verbatim embed:rst
+
+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.
+
+.. versionadded:: 18Sep2020
+
+\endverbatim
+ *
+ * \param  handle  pointer to a previously created LAMMPS instance
+ * \param  name    string with the name of the extracted property
+ * \return         integer constant encoding the data type of the property
+ *                 or -1 if not found.
+ * */
+
+int lammps_extract_atom_datatype(void *handle, const char *name)
+{
+  LAMMPS *lmp = (LAMMPS *) handle;
+  return lmp->atom->extract_datatype(name);
+}
+
+/* ---------------------------------------------------------------------- */
+
 /** Create N atoms from list of coordinates
  *
 \verbatim embed:rst
@@ -1472,11 +1583,13 @@ lists the available options.
      - ``int *``
      - Number of local data columns
 
-The pointers returned by this function are generally not persistent
-since the computed data may be re-distributed, re-allocated, and
-re-ordered at every invocation. It is advisable to re-invoke this
-function before the data is accessed, or make a copy if the data shall
-be used after other LAMMPS commands have been issued.
+.. warning::
+
+   The pointers returned by this function are generally not persistent
+   since the computed data may be re-distributed, re-allocated, and
+   re-ordered at every invocation. It is advisable to re-invoke this
+   function before the data is accessed, or make a copy if the data shall
+   be used after other LAMMPS commands have been issued.
 
 .. note::
 
@@ -1656,12 +1769,14 @@ The following table lists the available options.
      - ``int *``
      - Number of local data columns
 
-The pointers returned by this function for per-atom or local data are
-generally not persistent, since the computed data may be re-distributed,
-re-allocated, and re-ordered at every invocation of the fix.  It is thus
-advisable to re-invoke this function before the data is accessed, or
-make a copy, if the data shall be used after other LAMMPS commands have
-been issued.
+.. warning::
+
+   The pointers returned by this function for per-atom or local data are
+   generally not persistent, since the computed data may be re-distributed,
+   re-allocated, and re-ordered at every invocation of the fix.  It is thus
+   advisable to re-invoke this function before the data is accessed, or
+   make a copy, if the data shall be used after other LAMMPS commands have
+   been issued.
 
 .. note::
 
diff --git a/src/library.h b/src/library.h
index 4cf2dd6d19..2ddad86baa 100644
--- a/src/library.h
+++ b/src/library.h
@@ -40,6 +40,20 @@
 #include <inttypes.h>  /* for int64_t */
 #endif
 
+/** Data type constants for extracting data from atoms, computes and fixes
+ *
+ * Must be kept in sync with the equivalent constants in lammps.py */
+
+enum _LMP_DATATYPE_CONST {
+  LAMMPS_INT    = 0,       /*!< 32-bit integer (array) */
+  LAMMPS_INT_2D  = 1,      /*!< two-dimensional 32-bit integer array */
+  LAMMPS_DOUBLE = 2,       /*!< 64-bit double (array) */
+  LAMMPS_DOUBLE_2D = 3,    /*!< two-dimensional 64-bit double array */
+  LAMMPS_INT64 = 4,        /*!< 64-bit integer (array) */
+  LAMMPS_INT64_2D = 5,     /*!< two-dimensional 64-bit integer array */
+  LAMMPS_STRING = 6        /*!< C-String */
+};
+
 /** Style constants for extracting data from computes and fixes.
  *
  * Must be kept in sync with the equivalent constants in lammps.py */
@@ -113,6 +127,9 @@ int    lammps_extract_setting(void *handle, const char *keyword);
 void  *lammps_extract_global(void *handle, const char *name);
 void  *lammps_extract_atom(void *handle, const char *name);
 
+int lammps_extract_global_datatype(void *handle, const char *name);
+int lammps_extract_atom_datatype(void *handle, const char *name);
+
 #if !defined(LAMMPS_BIGBIG)
 int    lammps_create_atoms(void *handle, int n, int *id, int *type,
                            double *x, double *v, int *image, int bexpand);
diff --git a/src/lmptype.h b/src/lmptype.h
index e6b34c537f..e5dba94be0 100644
--- a/src/lmptype.h
+++ b/src/lmptype.h
@@ -102,6 +102,11 @@ typedef int64_t bigint;
 #define ATOTAGINT atoi
 #define ATOBIGINT ATOLL
 
+#define LAMMPS_TAGINT    LAMMPS_INT
+#define LAMMPS_TAGINT_2D LAMMPS_INT_2D
+#define LAMMPS_BIGINT    LAMMPS_INT64
+#define LAMMPS_BIGINT_2D LAMMPS_INT64_2D
+
 #define IMGMASK 1023
 #define IMGMAX 512
 #define IMGBITS 10
@@ -134,6 +139,11 @@ typedef int64_t bigint;
 #define ATOTAGINT ATOLL
 #define ATOBIGINT ATOLL
 
+#define LAMMPS_TAGINT    LAMMPS_INT64
+#define LAMMPS_TAGINT_2D LAMMPS_INT64_2D
+#define LAMMPS_BIGINT    LAMMPS_INT64
+#define LAMMPS_BIGINT_2D LAMMPS_INT64_2D
+
 #define IMGMASK 2097151
 #define IMGMAX 1048576
 #define IMGBITS 21
@@ -165,6 +175,11 @@ typedef int bigint;
 #define ATOTAGINT atoi
 #define ATOBIGINT atoi
 
+#define LAMMPS_TAGINT    LAMMPS_INT
+#define LAMMPS_TAGINT_2D LAMMPS_INT_2D
+#define LAMMPS_BIGINT    LAMMPS_INT
+#define LAMMPS_BIGINT_2D LAMMPS_INT_2D
+
 #define IMGMASK 1023
 #define IMGMAX 512
 #define IMGBITS 10
diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp
index eb0e1a7df6..11bdc28552 100644
--- a/unittest/c-library/test_library_properties.cpp
+++ b/unittest/c-library/test_library_properties.cpp
@@ -251,4 +251,90 @@ TEST_F(LibraryProperties, global)
     EXPECT_EQ((*b_ptr), 2);
     d_ptr = (double *)lammps_extract_global(lmp, "dt");
     EXPECT_DOUBLE_EQ((*d_ptr), 0.1);
+    int dtype = lammps_extract_global_datatype(lmp, "dt");
+    EXPECT_EQ(dtype, LAMMPS_DOUBLE);
 };
+
+class AtomProperties : public ::testing::Test {
+protected:
+    void *lmp;
+
+    AtomProperties(){};
+    ~AtomProperties() override{};
+
+    void SetUp() override
+    {
+        const char *args[] = {"LAMMPS_test", "-log",      "none",
+                              "-echo",       "screen",    "-nocite"};
+
+        char **argv = (char **)args;
+        int argc    = sizeof(args) / sizeof(char *);
+
+        ::testing::internal::CaptureStdout();
+        lmp                = lammps_open_no_mpi(argc, argv, NULL);
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        EXPECT_THAT(output, StartsWith("LAMMPS ("));
+        ::testing::internal::CaptureStdout();
+        lammps_command(lmp, "region box block 0 2 0 2 0 2");
+        lammps_command(lmp, "create_box 1 box");
+        lammps_command(lmp, "mass 1 3.0");
+        lammps_command(lmp, "create_atoms 1 single 1.0 1.0 1.5");
+        lammps_command(lmp, "create_atoms 1 single 0.2 0.1 0.1");
+        output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+    }
+    void TearDown() override
+    {
+        ::testing::internal::CaptureStdout();
+        lammps_close(lmp);
+        std::string output = ::testing::internal::GetCapturedStdout();
+        EXPECT_THAT(output, HasSubstr("Total wall time:"));
+        if (verbose) std::cout << output;
+        lmp = nullptr;
+    }
+};
+
+TEST_F(AtomProperties, invalid)
+{
+    ASSERT_EQ(lammps_extract_atom(lmp, "UNKNOWN"), nullptr);
+}
+
+TEST_F(AtomProperties, mass)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "mass"), LAMMPS_DOUBLE);
+    double * mass = (double *)lammps_extract_atom(lmp, "mass");
+    ASSERT_NE(mass, nullptr);
+    ASSERT_DOUBLE_EQ(mass[1], 3.0);
+}
+
+TEST_F(AtomProperties, id)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "id"), LAMMPS_TAGINT);
+    LAMMPS_NS::tagint * id = (LAMMPS_NS::tagint *)lammps_extract_atom(lmp, "id");
+    ASSERT_NE(id, nullptr);
+    ASSERT_EQ(id[0], 1);
+    ASSERT_EQ(id[1], 2);
+}
+
+TEST_F(AtomProperties, type)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "type"), LAMMPS_INT);
+    int * type = (int *)lammps_extract_atom(lmp, "type");
+    ASSERT_NE(type, nullptr);
+    ASSERT_EQ(type[0], 1);
+    ASSERT_EQ(type[1], 1);
+}
+
+TEST_F(AtomProperties, position)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "x"), LAMMPS_DOUBLE_2D);
+    double ** x = (double **)lammps_extract_atom(lmp, "x");
+    ASSERT_NE(x, nullptr);
+    EXPECT_DOUBLE_EQ(x[0][0], 1.0);
+    EXPECT_DOUBLE_EQ(x[0][1], 1.0);
+    EXPECT_DOUBLE_EQ(x[0][2], 1.5);
+    EXPECT_DOUBLE_EQ(x[1][0], 0.2);
+    EXPECT_DOUBLE_EQ(x[1][1], 0.1);
+    EXPECT_DOUBLE_EQ(x[1][2], 0.1);
+}
diff --git a/unittest/formats/CMakeLists.txt b/unittest/formats/CMakeLists.txt
index e8bfa64ef3..8b44bbc227 100644
--- a/unittest/formats/CMakeLists.txt
+++ b/unittest/formats/CMakeLists.txt
@@ -62,40 +62,40 @@ if (PKG_COMPRESS)
     set_tests_properties(DumpLocalGZ PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
     set_tests_properties(DumpLocalGZ PROPERTIES ENVIRONMENT "GZIP_BINARY=${GZIP_BINARY}")
 
-    if(Zstd_FOUND)
-        find_program(ZSTD_BINARY NAMES zstd)
+    find_package(PkgConfig REQUIRED)
+    pkg_check_modules(Zstd IMPORTED_TARGET libzstd>=1.4)
+    find_program(ZSTD_BINARY NAMES zstd)
 
-        if (ZSTD_BINARY)
-            add_executable(test_dump_atom_zstd test_dump_atom_zstd.cpp)
-            target_link_libraries(test_dump_atom_zstd PRIVATE lammps GTest::GMock GTest::GTest)
-            add_test(NAME DumpAtomZstd COMMAND test_dump_atom_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-            set_tests_properties(DumpAtomZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
-            set_tests_properties(DumpAtomZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
+    if(Zstd_FOUND AND ZSTD_BINARY)
+        add_executable(test_dump_atom_zstd test_dump_atom_zstd.cpp)
+        target_link_libraries(test_dump_atom_zstd PRIVATE lammps GTest::GMock GTest::GTest)
+        add_test(NAME DumpAtomZstd COMMAND test_dump_atom_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+        set_tests_properties(DumpAtomZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+        set_tests_properties(DumpAtomZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
 
-            add_executable(test_dump_custom_zstd test_dump_custom_zstd.cpp)
-            target_link_libraries(test_dump_custom_zstd PRIVATE lammps GTest::GMock GTest::GTest)
-            add_test(NAME DumpCustomZstd COMMAND test_dump_custom_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-            set_tests_properties(DumpCustomZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
-            set_tests_properties(DumpCustomZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
+        add_executable(test_dump_custom_zstd test_dump_custom_zstd.cpp)
+        target_link_libraries(test_dump_custom_zstd PRIVATE lammps GTest::GMock GTest::GTest)
+        add_test(NAME DumpCustomZstd COMMAND test_dump_custom_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+        set_tests_properties(DumpCustomZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+        set_tests_properties(DumpCustomZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
 
-            add_executable(test_dump_cfg_zstd test_dump_cfg_zstd.cpp)
-            target_link_libraries(test_dump_cfg_zstd PRIVATE lammps GTest::GMock GTest::GTest)
-            add_test(NAME DumpCfgZstd COMMAND test_dump_cfg_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-            set_tests_properties(DumpCfgZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
-            set_tests_properties(DumpCfgZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
+        add_executable(test_dump_cfg_zstd test_dump_cfg_zstd.cpp)
+        target_link_libraries(test_dump_cfg_zstd PRIVATE lammps GTest::GMock GTest::GTest)
+        add_test(NAME DumpCfgZstd COMMAND test_dump_cfg_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+        set_tests_properties(DumpCfgZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+        set_tests_properties(DumpCfgZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
 
-            add_executable(test_dump_xyz_zstd test_dump_xyz_zstd.cpp)
-            target_link_libraries(test_dump_xyz_zstd PRIVATE lammps GTest::GMock GTest::GTest)
-            add_test(NAME DumpXYZZstd COMMAND test_dump_xyz_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-            set_tests_properties(DumpXYZZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
-            set_tests_properties(DumpXYZZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
+        add_executable(test_dump_xyz_zstd test_dump_xyz_zstd.cpp)
+        target_link_libraries(test_dump_xyz_zstd PRIVATE lammps GTest::GMock GTest::GTest)
+        add_test(NAME DumpXYZZstd COMMAND test_dump_xyz_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+        set_tests_properties(DumpXYZZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+        set_tests_properties(DumpXYZZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
 
-            add_executable(test_dump_local_zstd test_dump_local_zstd.cpp)
-            target_link_libraries(test_dump_local_zstd PRIVATE lammps GTest::GMock GTest::GTest)
-            add_test(NAME DumpLocalZstd COMMAND test_dump_local_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-            set_tests_properties(DumpLocalZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
-            set_tests_properties(DumpLocalZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
-        endif()
+        add_executable(test_dump_local_zstd test_dump_local_zstd.cpp)
+        target_link_libraries(test_dump_local_zstd PRIVATE lammps GTest::GMock GTest::GTest)
+        add_test(NAME DumpLocalZstd COMMAND test_dump_local_zstd WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+        set_tests_properties(DumpLocalZstd PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+        set_tests_properties(DumpLocalZstd PROPERTIES ENVIRONMENT "ZSTD_BINARY=${ZSTD_BINARY}")
     endif()
 endif()
 
diff --git a/unittest/python/CMakeLists.txt b/unittest/python/CMakeLists.txt
index 574f7aab09..6082eb99b1 100644
--- a/unittest/python/CMakeLists.txt
+++ b/unittest/python/CMakeLists.txt
@@ -17,9 +17,9 @@ if (Python_EXECUTABLE)
   # prepare to augment the environment so that the LAMMPS python module and the shared library is found.
   set(PYTHON_TEST_ENVIRONMENT PYTHONPATH=${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH})
   if(APPLE)
-    list(APPEND PYTHON_TEST_ENVIRONMENT DYLD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$ENV{DYLD_LIBRARY_PATH})
+    list(APPEND PYTHON_TEST_ENVIRONMENT "DYLD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$ENV{DYLD_LIBRARY_PATH};LAMMPS_CMAKE_CACHE=${CMAKE_BINARY_DIR}/CMakeCache.txt")
   else()
-    list(APPEND PYTHON_TEST_ENVIRONMENT LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$ENV{LD_LIBRARY_PATH})
+    list(APPEND PYTHON_TEST_ENVIRONMENT "LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$ENV{LD_LIBRARY_PATH};LAMMPS_CMAKE_CACHE=${CMAKE_BINARY_DIR}/CMakeCache.txt")
   endif()
   if(LAMMPS_MACHINE)
     # convert from '_machine' to 'machine'
@@ -27,20 +27,43 @@ if (Python_EXECUTABLE)
     list(APPEND PYTHON_TEST_ENVIRONMENT LAMMPS_MACHINE_NAME=${LAMMPS_MACHINE_NAME})
   endif()
 
+  if(ENABLE_COVERAGE)
+    find_program(COVERAGE_BINARY coverage)
+    find_package_handle_standard_args(COVERAGE DEFAULT_MSG COVERAGE_BINARY)
+
+    if(COVERAGE_FOUND)
+        set(PYTHON_TEST_RUNNER ${Python_EXECUTABLE} -u ${COVERAGE_BINARY} run --parallel-mode --include=${LAMMPS_PYTHON_DIR}/lammps.py --omit=${LAMMPS_PYTHON_DIR}/install.py)
+    else()
+        set(PYTHON_TEST_RUNNER ${Python_EXECUTABLE} -u)
+    endif()
+  else()
+    set(PYTHON_TEST_RUNNER ${Python_EXECUTABLE})
+  endif()
+
   add_test(NAME PythonOpen
-           COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/python-open.py -v
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-open.py -v
            WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
   set_tests_properties(PythonOpen PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
 
   add_test(NAME PythonCommands
-           COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/python-commands.py -v
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-commands.py -v
            WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
   set_tests_properties(PythonCommands PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
 
   add_test(NAME PythonNumpy
-           COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/python-numpy.py -v
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-numpy.py -v
            WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
   set_tests_properties(PythonNumpy PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
+
+  add_test(NAME PythonCapabilities
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-capabilities.py -v
+           WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+  set_tests_properties(PythonCapabilities PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
+
+  add_test(NAME PythonPyLammps
+           COMMAND ${PYTHON_TEST_RUNNER} ${CMAKE_CURRENT_SOURCE_DIR}/python-pylammps.py -v
+           WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+  set_tests_properties(PythonPyLammps PROPERTIES ENVIRONMENT "${PYTHON_TEST_ENVIRONMENT}")
 else()
   message(STATUS "Skipping Tests for the LAMMPS Python Module: no suitable Python interpreter")
 endif()
diff --git a/unittest/python/python-capabilities.py b/unittest/python/python-capabilities.py
new file mode 100644
index 0000000000..27bdc9c561
--- /dev/null
+++ b/unittest/python/python-capabilities.py
@@ -0,0 +1,62 @@
+import sys,os,unittest
+from lammps import lammps
+
+class PythonCapabilities(unittest.TestCase):
+    def setUp(self):
+        machine = None
+        if 'LAMMPS_MACHINE_NAME' in os.environ:
+            machine=os.environ['LAMMPS_MACHINE_NAME']
+        self.lmp = lammps(name=machine,  cmdargs=['-nocite', '-log','none', '-echo','screen'])
+
+        if 'LAMMPS_CMAKE_CACHE' in os.environ:
+            self.cmake_cache = {}
+
+            with open(os.environ['LAMMPS_CMAKE_CACHE'], 'r') as f:
+                for line in f:
+                    line = line.strip()
+                    if not line or line.startswith('#') or line.startswith('//'): continue
+                    parts = line.split('=')
+                    key, value_type = parts[0].split(':')
+                    if len(parts) > 1:
+                        value = parts[1]
+                        if value_type == "BOOL":
+                            value = (value.upper() == "ON")
+                    else:
+                        value = None
+                    self.cmake_cache[key] = value
+
+    def tearDown(self):
+        del self.lmp
+
+    def test_version(self):
+        self.assertGreaterEqual(self.lmp.version(), 20200824)
+
+    def test_has_gzip_support(self):
+        self.assertEqual(self.lmp.has_gzip_support, self.cmake_cache['WITH_GZIP'])
+
+    def test_has_png_support(self):
+        self.assertEqual(self.lmp.has_png_support, self.cmake_cache['WITH_PNG'])
+
+    def test_has_jpeg_support(self):
+        self.assertEqual(self.lmp.has_jpeg_support, self.cmake_cache['WITH_JPEG'])
+
+    def test_has_ffmpeg_support(self):
+        self.assertEqual(self.lmp.has_ffmpeg_support, self.cmake_cache['WITH_FFMPEG'])
+
+    def test_installed_packages(self):
+        installed_packages = self.lmp.installed_packages
+        selected_packages = [key[4:] for key in self.cmake_cache.keys() if not key.startswith('PKG_CONFIG') and key.startswith('PKG_') and self.cmake_cache[key]]
+
+        for pkg in selected_packages:
+            self.assertIn(pkg, installed_packages)
+
+    def test_has_style(self):
+        self.assertTrue(self.lmp.has_style('pair', 'lj/cut'))
+        self.assertFalse(self.lmp.has_style('pair', 'lennard_jones'))
+
+    def test_available_styles(self):
+        pairs = self.lmp.available_styles('pair')
+        self.assertIn('lj/cut', pairs)
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py
index f6d48c8bba..0b853a207e 100644
--- a/unittest/python/python-commands.py
+++ b/unittest/python/python-commands.py
@@ -43,46 +43,46 @@ create_atoms 1 single &
     ##############################
     def testFile(self):
         """Test reading commands from a file"""
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,0)
         self.lmp.file(self.demo_file)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,1)
         self.lmp.file(self.cont_file)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,2)
 
     def testNoFile(self):
         """Test (not) reading commands from no file"""
         self.lmp.file(None)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,0)
 
     def testCommand(self):
         """Test executing individual commands"""
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,0)
         cmds = self.demo_input.splitlines()
         for cmd in cmds:
             self.lmp.command(cmd)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,1)
 
     def testCommandsList(self):
         """Test executing commands from list of strings"""
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,0)
         cmds = self.demo_input.splitlines()+self.cont_input.splitlines()
         self.lmp.commands_list(cmds)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,2)
 
     def testCommandsString(self):
         """Test executing block of commands from string"""
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,0)
         self.lmp.commands_string(self.demo_input+self.cont_input)
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,2)
 
 ##############################
diff --git a/unittest/python/python-numpy.py b/unittest/python/python-numpy.py
index 8fd41d2278..3c8ff9f512 100644
--- a/unittest/python/python-numpy.py
+++ b/unittest/python/python-numpy.py
@@ -1,7 +1,14 @@
 import sys,os,unittest
-from lammps import lammps, LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL, LMP_STYLE_ATOM, LMP_TYPE_VECTOR, LMP_TYPE_SCALAR, LMP_TYPE_ARRAY
+from lammps import lammps, LAMMPS_INT, LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL, LMP_STYLE_ATOM, LMP_TYPE_VECTOR, LMP_TYPE_SCALAR, LMP_TYPE_ARRAY
 from ctypes import c_void_p
 
+try:
+    import numpy
+    NUMPY_INSTALLED = True
+except ImportError:
+    NUMPY_INSTALLED = False
+
+@unittest.skipIf(not NUMPY_INSTALLED, "numpy is not available")
 class PythonNumpy(unittest.TestCase):
     def setUp(self):
         machine = None
@@ -25,7 +32,7 @@ class PythonNumpy(unittest.TestCase):
         self.lmp.command("create_atoms 1 single 1.0 1.0 1.0")
         self.lmp.command("create_atoms 1 single 1.0 1.0 1.5")
         self.lmp.command("compute coordsum all reduce sum x y z")
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,2)
         values = self.lmp.numpy.extract_compute("coordsum", LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR)
         self.assertEqual(len(values), 3)
@@ -43,7 +50,7 @@ class PythonNumpy(unittest.TestCase):
         self.lmp.command("create_atoms 1 single 1.0 1.0 1.0")
         self.lmp.command("create_atoms 1 single 1.0 1.0 1.5")
         self.lmp.command("compute ke all ke/atom")
-        natoms = int(self.lmp.get_natoms())
+        natoms = self.lmp.get_natoms()
         self.assertEqual(natoms,2)
         values = self.lmp.numpy.extract_compute("ke", LMP_STYLE_ATOM, LMP_TYPE_VECTOR)
         self.assertEqual(len(values), 2)
@@ -66,5 +73,67 @@ class PythonNumpy(unittest.TestCase):
         # TODO
         pass
 
+    def testExtractAtomDeprecated(self):
+        self.lmp.command("units lj")
+        self.lmp.command("atom_style atomic")
+        self.lmp.command("atom_modify map array")
+        self.lmp.command("region box block 0 2 0 2 0 2")
+        self.lmp.command("create_box 1 box")
+
+        x = [
+          1.0, 1.0, 1.0,
+          1.0, 1.0, 1.5
+        ]
+
+        types = [1, 1]
+
+        self.assertEqual(self.lmp.create_atoms(2, id=None, type=types, x=x), 2)
+        nlocal = self.lmp.extract_global("nlocal", LAMMPS_INT)
+        self.assertEqual(nlocal, 2)
+
+        ident = self.lmp.numpy.extract_atom_iarray("id", nlocal, dim=1)
+        self.assertEqual(len(ident), 2)
+
+        ntypes = self.lmp.extract_global("ntypes", LAMMPS_INT)
+        self.assertEqual(ntypes, 1)
+
+        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
+        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
+        self.assertEqual(len(x), 2)
+        self.assertTrue((x[0] == (1.0, 1.0, 1.0)).all())
+        self.assertTrue((x[1] == (1.0, 1.0, 1.5)).all())
+        self.assertEqual(len(v), 2)
+
+    def testExtractAtom(self):
+        self.lmp.command("units lj")
+        self.lmp.command("atom_style atomic")
+        self.lmp.command("atom_modify map array")
+        self.lmp.command("region box block 0 2 0 2 0 2")
+        self.lmp.command("create_box 1 box")
+
+        x = [
+          1.0, 1.0, 1.0,
+          1.0, 1.0, 1.5
+        ]
+
+        types = [1, 1]
+
+        self.assertEqual(self.lmp.create_atoms(2, id=None, type=types, x=x), 2)
+        nlocal = self.lmp.extract_global("nlocal")
+        self.assertEqual(nlocal, 2)
+
+        ident = self.lmp.numpy.extract_atom("id")
+        self.assertEqual(len(ident), 2)
+
+        ntypes = self.lmp.extract_global("ntypes")
+        self.assertEqual(ntypes, 1)
+
+        x = self.lmp.numpy.extract_atom("x")
+        v = self.lmp.numpy.extract_atom("v")
+        self.assertEqual(len(x), 2)
+        self.assertTrue((x[0] == (1.0, 1.0, 1.0)).all())
+        self.assertTrue((x[1] == (1.0, 1.0, 1.5)).all())
+        self.assertEqual(len(v), 2)
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/unittest/python/python-pylammps.py b/unittest/python/python-pylammps.py
new file mode 100644
index 0000000000..f703d206fa
--- /dev/null
+++ b/unittest/python/python-pylammps.py
@@ -0,0 +1,88 @@
+import sys,os,unittest
+from lammps import PyLammps
+
+class PythonPyLammps(unittest.TestCase):
+    def setUp(self):
+        machine = None
+        if 'LAMMPS_MACHINE_NAME' in os.environ:
+            machine=os.environ['LAMMPS_MACHINE_NAME']
+        self.pylmp = PyLammps(name=machine,  cmdargs=['-nocite', '-log','none', '-echo', 'screen'])
+        self.pylmp.units("lj")
+        self.pylmp.atom_style("atomic")
+        self.pylmp.atom_modify("map array")
+
+        if 'LAMMPS_CMAKE_CACHE' in os.environ:
+            self.cmake_cache = {}
+
+            with open(os.environ['LAMMPS_CMAKE_CACHE'], 'r') as f:
+                for line in f:
+                    line = line.strip()
+                    if not line or line.startswith('#') or line.startswith('//'): continue
+                    parts = line.split('=')
+                    key, value_type = parts[0].split(':')
+                    if len(parts) > 1:
+                        value = parts[1]
+                        if value_type == "BOOL":
+                            value = (value.upper() == "ON")
+                    else:
+                        value = None
+                    self.cmake_cache[key] = value
+
+    def tearDown(self):
+        self.pylmp.close()
+        del self.pylmp
+
+    def test_version(self):
+        self.assertGreaterEqual(self.pylmp.version(), 20200824)
+
+    def test_create_atoms(self):
+        self.pylmp.region("box block", 0, 2, 0, 2, 0, 2)
+        self.pylmp.create_box(1, "box")
+
+        x = [
+          1.0, 1.0, 1.0,
+          1.0, 1.0, 1.5
+        ]
+
+        types = [1, 1]
+
+        self.assertEqual(self.pylmp.lmp.create_atoms(2, id=None, type=types, x=x), 2)
+        self.assertEqual(self.pylmp.system.natoms, 2)
+        self.assertEqual(len(self.pylmp.atoms), 2)
+        self.assertEqual(self.pylmp.atoms[0].position, tuple(x[0:3]))
+        self.assertEqual(self.pylmp.atoms[1].position, tuple(x[3:6]))
+        self.assertEqual(self.pylmp.last_run, None)
+
+
+    def test_write_script(self):
+        outfile = 'in.test_write_script'
+        self.pylmp.write_script(outfile)
+        self.assertTrue(os.path.exists(outfile))
+        os.remove(outfile)
+
+    def test_runs(self):
+        self.pylmp.lattice("fcc", 0.8442),
+        self.pylmp.region("box block", 0, 4, 0, 4, 0, 4)
+        self.pylmp.create_box(1, "box")
+        self.pylmp.create_atoms(1, "box")
+        self.pylmp.mass(1, 1.0)
+        self.pylmp.velocity("all create", 1.44, 87287, "loop geom")
+        self.pylmp.pair_style("lj/cut", 2.5)
+        self.pylmp.pair_coeff(1, 1, 1.0, 1.0, 2.5)
+        self.pylmp.neighbor(0.3, "bin")
+        self.pylmp.neigh_modify("delay 0 every 20 check no")
+        self.pylmp.fix("1 all nve")
+        self.pylmp.variable("fx atom fx")
+        self.pylmp.run(10)
+
+        self.assertEqual(len(self.pylmp.runs), 1)
+        self.assertEqual(self.pylmp.last_run, self.pylmp.runs[0])
+        self.assertEqual(len(self.pylmp.last_run.thermo.Step), 2)
+        self.assertEqual(len(self.pylmp.last_run.thermo.Temp), 2)
+        self.assertEqual(len(self.pylmp.last_run.thermo.E_pair), 2)
+        self.assertEqual(len(self.pylmp.last_run.thermo.E_mol), 2)
+        self.assertEqual(len(self.pylmp.last_run.thermo.TotEng), 2)
+        self.assertEqual(len(self.pylmp.last_run.thermo.Press), 2)
+
+if __name__ == "__main__":
+    unittest.main()