diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst
index 992df6ba63..17b4715dc7 100644
--- a/doc/src/Developer_utils.rst
+++ b/doc/src/Developer_utils.rst
@@ -101,6 +101,9 @@ and parsing files or arguments.
 .. doxygenfunction:: split_words
    :project: progguide
 
+.. doxygenfunction:: split_lines
+   :project: progguide
+
 .. doxygenfunction:: strmatch
    :project: progguide
 
diff --git a/lib/gpu/lal_neighbor.cpp b/lib/gpu/lal_neighbor.cpp
index aabba49575..a0d2eaa8c3 100644
--- a/lib/gpu/lal_neighbor.cpp
+++ b/lib/gpu/lal_neighbor.cpp
@@ -740,6 +740,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
 
   // If binning on GPU, do this now
   if (_gpu_nbor==1) {
+    mn = _max_nbors;
     const numtyp i_cell_size=static_cast<numtyp>(1.0/_cell_size);
     const int neigh_block=_block_cell_id;
     const int GX=(int)ceil((float)nall/neigh_block);
diff --git a/python/lammps/__init__.py b/python/lammps/__init__.py
index 48839273c5..e6ffd779a9 100644
--- a/python/lammps/__init__.py
+++ b/python/lammps/__init__.py
@@ -15,7 +15,7 @@ from .pylammps import *
 
 # convert module string version to numeric version
 def get_version_number():
-    from datetime import datetime
+    import time
     from sys import version_info
     vstring = None
     if version_info.major == 3 and version_info.minor >= 8:
@@ -32,7 +32,7 @@ def get_version_number():
     if not vstring:
         return 0
 
-    d = datetime.strptime(vstring, "%d%b%Y")
-    return d.year*10000 + d.month*100 + d.day
+    t = time.strptime(vstring, "%d%b%Y")
+    return t.tm_year*10000 + t.tm_mon*100 + t.tm_mday
 
 __version__ = get_version_number()
diff --git a/python/lammps/core.py b/python/lammps/core.py
index be026d5e10..3118cb3d99 100644
--- a/python/lammps/core.py
+++ b/python/lammps/core.py
@@ -729,12 +729,11 @@ class lammps(object):
   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`
-    function of the C-library interface.  Unlike the C function
-    this method returns the value and not a pointer and thus can
-    only return the first value for keywords representing a list
-    of values.  The :cpp:func:`lammps_extract_global` documentation
-    includes a list of the supported keywords and their data types.
+    This is a wrapper around the :cpp:func:`lammps_extract_global` function
+    of the C-library interface.  Since there are no pointers in Python, this
+    method will - unlike the C function - return the value or a list 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, 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
@@ -746,12 +745,23 @@ class lammps(object):
     :type name:  string
     :param dtype: data type of the returned data (see :ref:`py_datatype_constants`)
     :type dtype:  int, optional
-    :return: value of the property or None
-    :rtype: int, float, or NoneType
+    :return: value of the property or list of values or None
+    :rtype: int, float, list, or NoneType
     """
+
     if dtype == LAMMPS_AUTODETECT:
       dtype = self.extract_global_datatype(name)
 
+    # set length of vector for items that are not a scalar
+    vec_dict = { 'boxlo':3, 'boxhi':3, 'sublo':3, 'subhi':3,
+                 'sublo_lambda':3, 'subhi_lambda':3, 'periodicity':3 }
+    if name in vec_dict:
+      veclen = vec_dict[name]
+    elif name == 'respa_dt':
+      veclen = self.extract_global('respa_levels',LAMMPS_INT)
+    else:
+      veclen = 1
+
     if name: name = name.encode()
     else: return None
 
@@ -766,14 +776,19 @@ class lammps(object):
       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])
+      if dtype == LAMMPS_STRING:
+        return ptr.decode('utf-8')
+      if veclen > 1:
+        result = []
+        for i in range(0,veclen):
+          result.append(target_type(ptr[i]))
+        return result
+      else: return target_type(ptr[0])
     return None
 
-
   # -------------------------------------------------------------------------
   # extract per-atom info datatype
 
diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp
index c428163688..f12c1774d2 100644
--- a/src/KIM/kim_interactions.cpp
+++ b/src/KIM/kim_interactions.cpp
@@ -12,11 +12,11 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing authors: Axel Kohlmeyer (Temple U),
-                         Ryan S. Elliott (UMN)
-                         Ellad B. Tadmor (UMN)
-                         Ronald Miller   (Carleton U)
-                         Yaser Afshar (UMN)
+   Contributing authors: Axel Kohlmeyer  (Temple U),
+                         Ryan S. Elliott (UMN),
+                         Ellad B. Tadmor (UMN),
+                         Ronald Miller   (Carleton U),
+                         Yaser Afshar    (UMN)
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
@@ -65,6 +65,7 @@
 #include "error.h"
 #include "fix_store_kim.h"
 #include "input.h"
+#include "variable.h"
 #include "modify.h"
 #include "update.h"
 
@@ -129,6 +130,7 @@ void KimInteractions::do_setup(int narg, char **arg)
                     "=======\n");
 
   if (simulatorModel) {
+    auto first_visit = input->variable->find("kim_update");
     if (!fixed_types) {
       std::string atom_type_sym_list =
         fmt::format("{}", fmt::join(arg, arg + narg, " "));
@@ -198,15 +200,17 @@ void KimInteractions::do_setup(int narg, char **arg)
 
       const std::string sim_field_str(sim_field);
       if (sim_field_str == "model-defn") {
+        if (first_visit < 0) input->one("variable kim_update equal 0");
+        else  input->one("variable kim_update equal 1");
         if (domain->periodicity[0] &&
             domain->periodicity[1] &&
             domain->periodicity[2])
           input->one("variable kim_periodic equal 1");
-        else if (domain->periodicity[0] &&
-                 domain->periodicity[1] &&
+        else if (!domain->periodicity[0] &&
+                 !domain->periodicity[1] &&
                  !domain->periodicity[2])
-          input->one("variable kim_periodic equal 2");
-        else input->one("variable kim_periodic equal 0");
+          input->one("variable kim_periodic equal 0");
+        else input->one("variable kim_periodic equal 2");
 
         // KIM Simulator Model has a Model definition
         no_model_definition = false;
@@ -241,6 +245,7 @@ void KimInteractions::do_setup(int narg, char **arg)
 
   } else {
 
+
     // not a simulator model. issue pair_style and pair_coeff commands.
 
     if (fixed_types)
diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp
index d1f602a1ea..4c43ca3744 100644
--- a/src/PYTHON/python_impl.cpp
+++ b/src/PYTHON/python_impl.cpp
@@ -55,6 +55,16 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp)
 
   nfunc = 0;
   pfuncs = nullptr;
+
+  // check for PYTHONUNBUFFERED environment variable
+  const char * PYTHONUNBUFFERED = getenv("PYTHONUNBUFFERED");
+
+  if (PYTHONUNBUFFERED != nullptr && strcmp(PYTHONUNBUFFERED, "1") == 0) {
+    // Python Global configuration variable
+    // Force the stdout and stderr streams to be unbuffered.
+    Py_UnbufferedStdioFlag = 1;
+  }
+
   // one-time initialization of Python interpreter
   // pyMain stores pointer to main module
   external_interpreter = Py_IsInitialized();
@@ -496,6 +506,7 @@ int PythonImpl::create_entry(char *name)
                  "cannot be used unless output is a string");
     pfuncs[ifunc].length_longstr = length_longstr;
     pfuncs[ifunc].longstr = new char[length_longstr+1];
+    pfuncs[ifunc].longstr[length_longstr] = '\0';
   }
 
   if (strstr(ostr,"v_") != ostr) error->all(FLERR,"Invalid python command");
diff --git a/src/USER-CGDNA/pair_oxdna2_coaxstk.cpp b/src/USER-CGDNA/pair_oxdna2_coaxstk.cpp
index f0b3386478..47d16c4806 100644
--- a/src/USER-CGDNA/pair_oxdna2_coaxstk.cpp
+++ b/src/USER-CGDNA/pair_oxdna2_coaxstk.cpp
@@ -551,7 +551,7 @@ void PairOxdna2Coaxstk::coeff(int narg, char **arg)
 {
   int count;
 
-  if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk");
+  if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/coaxstk");
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
@@ -673,7 +673,7 @@ void PairOxdna2Coaxstk::coeff(int narg, char **arg)
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk");
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/coaxstk");
 
 }
 
diff --git a/src/USER-CGDNA/pair_oxdna2_dh.cpp b/src/USER-CGDNA/pair_oxdna2_dh.cpp
index d877e0bf9a..6b0df2cde6 100644
--- a/src/USER-CGDNA/pair_oxdna2_dh.cpp
+++ b/src/USER-CGDNA/pair_oxdna2_dh.cpp
@@ -273,7 +273,7 @@ void PairOxdna2Dh::coeff(int narg, char **arg)
 {
   int count;
 
-  if (narg != 5) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
+  if (narg != 5) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/dh");
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
@@ -356,7 +356,7 @@ void PairOxdna2Dh::coeff(int narg, char **arg)
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/dh");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/USER-CGDNA/pair_oxrna2_stk.cpp b/src/USER-CGDNA/pair_oxrna2_stk.cpp
index 33480aa380..3b02a57501 100644
--- a/src/USER-CGDNA/pair_oxrna2_stk.cpp
+++ b/src/USER-CGDNA/pair_oxrna2_stk.cpp
@@ -1022,7 +1022,7 @@ void PairOxrna2Stk::coeff(int narg, char **arg)
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk");
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/stk");
 
 }
 
diff --git a/src/USER-CGDNA/pair_oxrna2_xstk.cpp b/src/USER-CGDNA/pair_oxrna2_xstk.cpp
index c3b97bfa9c..71c6ca0356 100644
--- a/src/USER-CGDNA/pair_oxrna2_xstk.cpp
+++ b/src/USER-CGDNA/pair_oxrna2_xstk.cpp
@@ -576,7 +576,7 @@ void PairOxrna2Xstk::coeff(int narg, char **arg)
 {
   int count;
 
-  if (narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk");
+  if (narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/xstk");
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
@@ -707,7 +707,7 @@ void PairOxrna2Xstk::coeff(int narg, char **arg)
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk");
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/xstk");
 
 }
 
diff --git a/src/library.cpp b/src/library.cpp
index 8de36a299f..a950d31e4e 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -31,12 +31,14 @@
 #include "group.h"
 #include "info.h"
 #include "input.h"
+#include "integrate.h"
 #include "memory.h"
 #include "modify.h"
 #include "molecule.h"
 #include "neigh_list.h"
 #include "neighbor.h"
 #include "region.h"
+#include "respa.h"
 #include "output.h"
 #include "thermo.h"
 #include "timer.h"
@@ -984,12 +986,14 @@ to then decide how to cast the (void*) pointer and access the data.
  * \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)
+int lammps_extract_global_datatype(void * /*handle*/, const char *name)
 {
   if (strcmp(name,"dt") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"ntimestep") == 0) return LAMMPS_BIGINT;
   if (strcmp(name,"atime") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"atimestep") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"respa_levels") == 0) return LAMMPS_INT;
+  if (strcmp(name,"respa_dt") == 0) return LAMMPS_DOUBLE;
 
   if (strcmp(name,"boxlo") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"boxhi") == 0) return LAMMPS_DOUBLE;
@@ -1116,6 +1120,14 @@ report the "native" data type.  The following tables are provided:
      - bigint
      - 1
      - the number of the timestep when "atime" was last updated.
+   * - respa_levels
+     - int
+     - 1
+     - number of r-RESPA levels. See :doc:`run_style`.
+   * - respa_dt
+     - double
+     - number of r-RESPA levels
+     - length of the time steps with r-RESPA. See :doc:`run_style`.
 
 .. _extract_box_settings:
 
@@ -1366,12 +1378,22 @@ void *lammps_extract_global(void *handle, const char *name)
   if (strcmp(name,"atime") == 0) return (void *) &lmp->update->atime;
   if (strcmp(name,"atimestep") == 0) return (void *) &lmp->update->atimestep;
 
+  if (utils::strmatch(lmp->update->integrate_style,"^respa")) {
+    Respa *respa = (Respa *)lmp->update->integrate;
+    if (strcmp(name,"respa_levels") == 0) return (void *) &respa->nlevels;
+    if (strcmp(name,"respa_dt") == 0) return (void *) respa->step;
+  }
   if (strcmp(name,"boxlo") == 0) return (void *) lmp->domain->boxlo;
   if (strcmp(name,"boxhi") == 0) return (void *) lmp->domain->boxhi;
   if (strcmp(name,"sublo") == 0) return (void *) lmp->domain->sublo;
   if (strcmp(name,"subhi") == 0) return (void *) lmp->domain->subhi;
-  if (strcmp(name,"sublo_lambda") == 0) return (void *) lmp->domain->sublo_lamda;
-  if (strcmp(name,"subhi_lambda") == 0) return (void *) lmp->domain->subhi_lamda;
+  // these are only valid for a triclinic cell
+  if (lmp->domain->triclinic) {
+    if (strcmp(name,"sublo_lambda") == 0)
+      return (void *) lmp->domain->sublo_lamda;
+    if (strcmp(name,"subhi_lambda") == 0)
+      return (void *) lmp->domain->subhi_lamda;
+  }
   if (strcmp(name,"boxxlo") == 0) return (void *) &lmp->domain->boxlo[0];
   if (strcmp(name,"boxxhi") == 0) return (void *) &lmp->domain->boxhi[0];
   if (strcmp(name,"boxylo") == 0) return (void *) &lmp->domain->boxlo[1];
diff --git a/src/utils.cpp b/src/utils.cpp
index b5576b9f27..e733d7eaae 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -851,6 +851,14 @@ std::vector<std::string> utils::split_words(const std::string &text)
   return list;
 }
 
+/* ----------------------------------------------------------------------
+   Convert multi-line string into lines
+------------------------------------------------------------------------- */
+std::vector<std::string> utils::split_lines(const std::string &text)
+{
+  return Tokenizer(text, "\n").as_vector();
+}
+
 /* ----------------------------------------------------------------------
    Return whether string is a valid integer number
 ------------------------------------------------------------------------- */
diff --git a/src/utils.h b/src/utils.h
index eab81f1343..ec4dd6ae85 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -321,6 +321,12 @@ namespace LAMMPS_NS {
 
     std::vector<std::string> split_words(const std::string &text);
 
+    /** Take multi-line text and split into lines
+     *
+     * \param text string that should be split
+     * \return STL vector with the lines */
+    std::vector<std::string> split_lines(const std::string &text);
+
     /** Check if string can be converted to valid integer
      *
      * \param str string that should be checked
diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt
index 02300b5ea2..60e83e7e16 100644
--- a/unittest/force-styles/CMakeLists.txt
+++ b/unittest/force-styles/CMakeLists.txt
@@ -124,7 +124,7 @@ file(GLOB FIX_TIMESTEP_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/fix-tim
 foreach(TEST ${FIX_TIMESTEP_TESTS})
   string(REGEX REPLACE "^.*fix-timestep-(.*)\.yaml" "FixTimestep:\\1" TNAME ${TEST})
   add_test(NAME ${TNAME} COMMAND test_fix_timestep ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
+  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH}")
 endforeach()
 
 # dihedral style tester
diff --git a/unittest/force-styles/tests/fix-timestep-python_move_nve.yaml b/unittest/force-styles/tests/fix-timestep-python_move_nve.yaml
new file mode 100644
index 0000000000..9c78dbe416
--- /dev/null
+++ b/unittest/force-styles/tests/fix-timestep-python_move_nve.yaml
@@ -0,0 +1,73 @@
+---
+lammps_version: 10 Mar 2021
+date_generated: Wed Mar 24 18:57:26 202
+epsilon: 9e-12
+prerequisites: ! |
+  atom full
+  fix python/move
+pre_commands: ! ""
+post_commands: ! |
+  fix test all python/move py_nve.NVE
+input_file: in.fourmol
+natoms: 29
+run_pos: ! |2
+    1 -2.7045559775384037e-01  2.4912159905679729e+00 -1.6695851791541885e-01
+    2  3.1004029573899528e-01  2.9612354631094391e+00 -8.5466363037021464e-01
+    3 -7.0398551400789477e-01  1.2305509955830618e+00 -6.2777526944456274e-01
+    4 -1.5818159336499285e+00  1.4837407818929933e+00 -1.2538710836062004e+00
+    5 -9.0719763672789266e-01  9.2652103885675297e-01  3.9954210488374786e-01
+    6  2.4831720524855985e-01  2.8313021497871271e-01 -1.2314233331711453e+00
+    7  3.4143527641386412e-01 -2.2646551041391422e-02 -2.5292291414903052e+00
+    8  1.1743552229100009e+00 -4.8863228565853950e-01 -6.3783432910825522e-01
+    9  1.3800524229500313e+00 -2.5274721030406683e-01  2.8353985887095157e-01
+   10  2.0510765220543883e+00 -1.4604063740302866e+00 -9.8323745081712954e-01
+   11  1.7878031944442556e+00 -1.9921863272948861e+00 -1.8890602447625777e+00
+   12  3.0063007039340053e+00 -4.9013350496963298e-01 -1.6231898107386229e+00
+   13  4.0515402959192999e+00 -8.9202011606653986e-01 -1.6400005529924957e+00
+   14  2.6066963345543819e+00 -4.1789253965514150e-01 -2.6634003608794394e+00
+   15  2.9695287185712913e+00  5.5422613165234036e-01 -1.2342022021790127e+00
+   16  2.6747029695228521e+00 -2.4124119054564295e+00 -2.3435746150616148e-02
+   17  2.2153577785283796e+00 -2.0897985186907717e+00  1.1963150794479436e+00
+   18  2.1369701704115704e+00  3.0158507413630606e+00 -3.5179348337215015e+00
+   19  1.5355837136087378e+00  2.6255292355375675e+00 -4.2353987779879052e+00
+   20  2.7727573005678776e+00  3.6923910449610169e+00 -3.9330842459133493e+00
+   21  4.9040128073204299e+00 -4.0752348172957946e+00 -3.6210314709891711e+00
+   22  4.3582355554440841e+00 -4.2126119427287048e+00 -4.4612844196314052e+00
+   23  5.7439382849307599e+00 -3.5821957939275029e+00 -3.8766361295935821e+00
+   24  2.0689243582422630e+00  3.1513346907271012e+00  3.1550389754828800e+00
+   25  1.3045351331492134e+00  3.2665125705842848e+00  2.5111855257433504e+00
+   26  2.5809237402711274e+00  4.0117602605482832e+00  3.2212060529089896e+00
+   27 -1.9611343130357228e+00 -4.3563411931359752e+00  2.1098293115523705e+00
+   28 -2.7473562684513411e+00 -4.0200819932379330e+00  1.5830052163433954e+00
+   29 -1.3126000191359855e+00 -3.5962518039482929e+00  2.2746342468737835e+00
+run_vel: ! |2
+    1  8.1705744183262364e-03  1.6516406176274288e-02  4.7902264318912978e-03
+    2  5.4501493445687759e-03  5.1791699408496334e-03 -1.4372931530376651e-03
+    3 -8.2298292722385643e-03 -1.2926551614621381e-02 -4.0984181178163881e-03
+    4 -3.7699042590093588e-03 -6.5722892098813799e-03 -1.1184640360133230e-03
+    5 -1.1021961004346586e-02 -9.8906780939335987e-03 -2.8410737829284395e-03
+    6 -3.9676663166400034e-02  4.6817061464710256e-02  3.7148491979476124e-02
+    7  9.1033953013898580e-04 -1.0128524411938794e-02 -5.1568251805019748e-02
+    8  7.9064712058855707e-03 -3.3507254552631767e-03  3.4557098492564629e-02
+    9  1.5644176117320923e-03  3.7365546102722164e-03  1.5047408822037646e-02
+   10  2.9201446820573174e-02 -2.9249578745486147e-02 -1.5018077424322538e-02
+   11 -4.7835961513517560e-03 -3.7481385134185206e-03 -2.3464104142290089e-03
+   12  2.2696451841920521e-03 -3.4774154398129479e-04 -3.0640770327796806e-03
+   13  2.7531740451953168e-03  5.8171061612840667e-03 -7.9467454022160518e-04
+   14  3.5246182371994252e-03 -5.7939995585585468e-03 -3.9478431172751344e-03
+   15 -1.8547943640122894e-03 -5.8554729942777743e-03  6.2938485140538649e-03
+   16  1.8681499973445245e-02 -1.3262466204585335e-02 -4.5638651457003243e-02
+   17 -1.2896269981100382e-02  9.7527665265956451e-03  3.7296535360836762e-02
+   18 -8.0065794848261610e-04 -8.6270473212554395e-04 -1.4483040697508738e-03
+   19  1.2452390836182623e-03 -2.5061097118772701e-03  7.2998631009712975e-03
+   20  3.5930060229597042e-03  3.6938860309252966e-03  3.2322732687893028e-03
+   21 -1.4689220370766550e-03 -2.7352129761527741e-04  7.0581624215243391e-04
+   22 -7.0694199254630339e-03 -4.2577148924878554e-03  2.8079117614251796e-04
+   23  6.0446963117374913e-03 -1.4000131614795382e-03  2.5819754847014255e-03
+   24  3.1926367902287940e-04 -9.9445664749276438e-04  1.4999996959365452e-04
+   25  1.3789754514814662e-04 -4.4335894884532569e-03 -8.1808136725080281e-04
+   26  2.0485904035217549e-03  2.7813358633835984e-03  4.3245727149206692e-03
+   27  4.5604120293369857e-04 -1.0305523026921137e-03  2.1188058381358511e-04
+   28 -6.2544520861855116e-03  1.4127711176146942e-03 -1.8429821884794269e-03
+   29  6.4110631534401762e-04  3.1273432719593867e-03  3.7253671105656715e-03
+...
diff --git a/unittest/force-styles/tests/py_nve.py b/unittest/force-styles/tests/py_nve.py
new file mode 100644
index 0000000000..57592ea074
--- /dev/null
+++ b/unittest/force-styles/tests/py_nve.py
@@ -0,0 +1,104 @@
+from __future__ import print_function
+from lammps import lammps
+
+class LAMMPSFix(object):
+    def __init__(self, ptr, group_name="all"):
+        self.lmp = lammps(ptr=ptr)
+        self.group_name = group_name
+
+class LAMMPSFixMove(LAMMPSFix):
+    def __init__(self, ptr, group_name="all"):
+        super(LAMMPSFixMove, self).__init__(ptr, group_name)
+
+    def init(self):
+        pass
+
+    def initial_integrate(self, vflag):
+        pass
+
+    def final_integrate(self):
+        pass
+
+    def initial_integrate_respa(self, vflag, ilevel, iloop):
+        pass
+
+    def final_integrate_respa(self, ilevel, iloop):
+        pass
+
+    def reset_dt(self):
+        pass
+
+
+class NVE(LAMMPSFixMove):
+    """ Python implementation of fix/nve """
+    def __init__(self, ptr, group_name="all"):
+        super(NVE, self).__init__(ptr)
+        assert(self.group_name == "all")
+        self._step_respa = None
+
+    def init(self):
+        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
+
+    @property
+    def step_respa(self):
+        if not self._step_respa:
+            self._step_respa = self.lmp.extract_global("respa_dt")
+        return self._step_respa
+
+    def initial_integrate(self, vflag):
+        nlocal = self.lmp.extract_global("nlocal")
+        mass = self.lmp.extract_atom("mass")
+        atype = self.lmp.extract_atom("type")
+        x = self.lmp.extract_atom("x")
+        v = self.lmp.extract_atom("v")
+        f = self.lmp.extract_atom("f")
+
+        for i in range(nlocal):
+            dtfm = self.dtf / mass[int(atype[i])]
+            v[i][0] += dtfm * f[i][0]
+            v[i][1] += dtfm * f[i][1]
+            v[i][2] += dtfm * f[i][2]
+            x[i][0] += self.dtv * v[i][0]
+            x[i][1] += self.dtv * v[i][1]
+            x[i][2] += self.dtv * v[i][2]
+
+    def final_integrate(self):
+        nlocal = self.lmp.extract_global("nlocal")
+        mass = self.lmp.extract_atom("mass")
+        atype = self.lmp.extract_atom("type")
+        v = self.lmp.extract_atom("v")
+        f = self.lmp.extract_atom("f")
+
+        for i in range(nlocal):
+            dtfm = self.dtf / mass[int(atype[i])]
+            v[i][0] += dtfm * f[i][0]
+            v[i][1] += dtfm * f[i][1]
+            v[i][2] += dtfm * f[i][2]
+
+    def initial_integrate_respa(self, vflag, ilevel, iloop):
+        ftm2v = self.lmp.extract_global("ftm2v")
+        self.dtv = self.step_respa[ilevel]
+        self.dtf = 0.5 * self.step_respa[ilevel] * ftm2v
+
+        # innermost level - NVE update of v and x
+        # all other levels - NVE update of v
+
+        if ilevel == 0:
+            self.initial_integrate(vflag)
+        else:
+            self.final_integrate()
+
+    def final_integrate_respa(self, ilevel, iloop):
+        ftm2v = self.lmp.extract_global("ftm2v")
+        self.dtf = 0.5 * self.step_respa[ilevel] * ftm2v
+        self.final_integrate()
+
+    def reset_dt(self):
+        dt = self.lmp.extract_global("dt")
+        ftm2v = self.lmp.extract_global("ftm2v")
+        self.dtv = dt;
+        self.dtf = 0.5 * dt * ftm2v;
diff --git a/unittest/python/CMakeLists.txt b/unittest/python/CMakeLists.txt
index d508602c93..d1db17c941 100644
--- a/unittest/python/CMakeLists.txt
+++ b/unittest/python/CMakeLists.txt
@@ -8,7 +8,7 @@ add_executable(test_python_package test_python_package.cpp)
 target_link_libraries(test_python_package PRIVATE lammps GTest::GMock GTest::GTest)
 target_compile_definitions(test_python_package PRIVATE -DTEST_INPUT_FOLDER=${TEST_INPUT_FOLDER})
 add_test(NAME PythonPackage COMMAND test_python_package WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-set_tests_properties(PythonPackage PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH}")
+set_tests_properties(PythonPackage PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH};PYTHONUNBUFFERED=1")
 
 # we must have shared libraries enabled for testing the python module
 if(NOT BUILD_SHARED_LIBS)
diff --git a/unittest/python/func.py b/unittest/python/func.py
index 27704660a6..cf8db41670 100644
--- a/unittest/python/func.py
+++ b/unittest/python/func.py
@@ -4,6 +4,16 @@ from __future__ import print_function
 def square(val):
     return val*val
 
+def bool_to_val(txt):
+    if txt.upper() in ["TRUE", "YES"]:
+        return 1.0
+    return 0.0
+
+def val_to_bool(val):
+    if val != 0:
+        return "True"
+    return "False"
+
 def printnum():
     print("2.25")
 
@@ -11,8 +21,10 @@ def printtxt():
     print("sometext")
 
 def getidxvar(lmpptr):
-    from lammps import lammps, LMP_VAR_EQUAL
+    from lammps import lammps
     lmp = lammps(ptr=lmpptr)
-
-    val = lmp.extract_variable("idx",None,LMP_VAR_EQUAL)
+    val = lmp.extract_variable("idx")
     print(val)
+
+def longstr():
+    return "Lorem ipsum dolor sit amet, consectetur adipiscing elit. Praesent metus."
diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py
index c82e2fdb26..3661feb8a0 100644
--- a/unittest/python/python-commands.py
+++ b/unittest/python/python-commands.py
@@ -259,15 +259,43 @@ create_atoms 1 single &
             result = self.lmp.get_thermo(key)
             self.assertEqual(value, result, key)
 
-    def test_extract_global_double(self):
+    def test_extract_global(self):
         self.lmp.command("region box block -1 1 -2 2 -3 3")
         self.lmp.command("create_box 1 box")
+        self.assertEqual(self.lmp.extract_global("units"), "lj")
+        self.assertEqual(self.lmp.extract_global("ntimestep"), 0)
+        self.assertEqual(self.lmp.extract_global("dt"), 0.005)
+
         self.assertEqual(self.lmp.extract_global("boxxlo"), -1.0)
         self.assertEqual(self.lmp.extract_global("boxxhi"), 1.0)
         self.assertEqual(self.lmp.extract_global("boxylo"), -2.0)
         self.assertEqual(self.lmp.extract_global("boxyhi"), 2.0)
         self.assertEqual(self.lmp.extract_global("boxzlo"), -3.0)
         self.assertEqual(self.lmp.extract_global("boxzhi"), 3.0)
+        self.assertEqual(self.lmp.extract_global("boxlo"), [-1.0, -2.0, -3.0])
+        self.assertEqual(self.lmp.extract_global("boxhi"), [1.0, 2.0, 3.0])
+        self.assertEqual(self.lmp.extract_global("sublo"), [-1.0, -2.0, -3.0])
+        self.assertEqual(self.lmp.extract_global("subhi"), [1.0, 2.0, 3.0])
+        self.assertEqual(self.lmp.extract_global("periodicity"), [1,1,1])
+        self.assertEqual(self.lmp.extract_global("triclinic"), 0)
+        self.assertEqual(self.lmp.extract_global("sublo_lambda"), None)
+        self.assertEqual(self.lmp.extract_global("subhi_lambda"), None)
+        self.assertEqual(self.lmp.extract_global("respa_levels"), None)
+        self.assertEqual(self.lmp.extract_global("respa_dt"), None)
+
+        # set and initialize r-RESPA
+        self.lmp.command("run_style respa 3 5 2 pair 2 kspace 3")
+        self.lmp.command("mass * 1.0")
+        self.lmp.command("run 1 post no")
+        self.assertEqual(self.lmp.extract_global("ntimestep"), 1)
+        self.assertEqual(self.lmp.extract_global("respa_levels"), 3)
+        self.assertEqual(self.lmp.extract_global("respa_dt"), [0.0005, 0.0025, 0.005])
+
+        # checks only for triclinic boxes
+        self.lmp.command("change_box all triclinic")
+        self.assertEqual(self.lmp.extract_global("triclinic"), 1)
+        self.assertEqual(self.lmp.extract_global("sublo_lambda"), [0.0, 0.0, 0.0])
+        self.assertEqual(self.lmp.extract_global("subhi_lambda"), [1.0, 1.0, 1.0])
 
 ##############################
 if __name__ == "__main__":
diff --git a/unittest/python/run.py b/unittest/python/run.py
new file mode 100644
index 0000000000..7cdb205f50
--- /dev/null
+++ b/unittest/python/run.py
@@ -0,0 +1,2 @@
+from __future__ import print_function
+print("Lorem ipsum dolor sit amet, consectetur adipiscing elit. Praesent metus.")
diff --git a/unittest/python/test_python_package.cpp b/unittest/python/test_python_package.cpp
index cba77ee2b0..c240d23875 100644
--- a/unittest/python/test_python_package.cpp
+++ b/unittest/python/test_python_package.cpp
@@ -14,19 +14,25 @@
 #include "atom.h"
 #include "info.h"
 #include "input.h"
+#include "variable.h"
+#include "library.h"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 
 #include <cmath>
 #include <cstring>
 #include <vector>
+#include <functional>
+
+#include "../testing/core.h"
+#include "../testing/systems/melt.h"
 
 // location of '*.py' files required by tests
 #define STRINGIFY(val) XSTR(val)
 #define XSTR(val) #val
 std::string INPUT_FOLDER = STRINGIFY(TEST_INPUT_FOLDER);
 
-// whether to print verbose output (i.e. not capturing LAMMPS screen output).
+const char * LOREM_IPSUM = "Lorem ipsum dolor sit amet, consectetur adipiscing elit. Praesent metus.";
 bool verbose = false;
 
 using LAMMPS_NS::utils::split_words;
@@ -34,95 +40,287 @@ using LAMMPS_NS::utils::split_words;
 namespace LAMMPS_NS {
 using ::testing::MatchesRegex;
 using ::testing::StrEq;
+using ::testing::Eq;
+using ::testing::HasSubstr;
 
-class PythonPackageTest : public ::testing::Test {
+class PythonPackageTest : public LAMMPSTest {
 protected:
-    LAMMPS *lmp;
-    Info *info;
-
-    void SetUp() override
+    void InitSystem() override
     {
-        const char *args[] = {"PythonPackageTest", "-log", "none", "-echo", "screen", "-nocite"};
-        char **argv        = (char **)args;
-        int argc           = sizeof(args) / sizeof(char *);
-        if (!verbose) ::testing::internal::CaptureStdout();
-        lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
-        if (!verbose) ::testing::internal::GetCapturedStdout();
-        ASSERT_NE(lmp, nullptr);
-        info = new Info(lmp);
-        if (!verbose) ::testing::internal::CaptureStdout();
-        lmp->input->one("units real");
-        lmp->input->one("dimension 3");
-        lmp->input->one("region box block -4 4 -4 4 -4 4");
-        lmp->input->one("create_box 1 box");
-        lmp->input->one("create_atoms 1 single  0.0  0.0 0.0    units box");
-        lmp->input->one("create_atoms 1 single  1.9 -1.9 1.9999 units box");
-        lmp->input->one("pair_style zero 2.0");
-        lmp->input->one("pair_coeff * *");
-        lmp->input->one("mass * 1.0");
-        lmp->input->one("variable input_dir index " + INPUT_FOLDER);
-        if (!verbose) ::testing::internal::GetCapturedStdout();
-    }
+        if (!info->has_package("PYTHON")) GTEST_SKIP();
 
-    void TearDown() override
-    {
-        if (!verbose) ::testing::internal::CaptureStdout();
-        delete info;
-        delete lmp;
-        if (!verbose) ::testing::internal::GetCapturedStdout();
+        HIDE_OUTPUT([&] {
+            command("units real");
+            command("dimension 3");
+            command("region box block -4 4 -4 4 -4 4");
+            command("create_box 1 box");
+            command("create_atoms 1 single  0.0  0.0 0.0    units box");
+            command("create_atoms 1 single  1.9 -1.9 1.9999 units box");
+            command("pair_style zero 2.0");
+            command("pair_coeff * *");
+            command("mass * 1.0");
+            command("variable input_dir index " + INPUT_FOLDER);
+        });
     }
 };
 
-TEST_F(PythonPackageTest, python_invoke)
+class FixPythonInvokeTest : public MeltTest {
+protected:
+    void InitSystem() override
+    {
+        if (!info->has_package("PYTHON")) GTEST_SKIP();
+
+        MeltTest::InitSystem();
+    }
+};
+
+TEST_F(PythonPackageTest, InvokeFunctionFromFile)
 {
-    if (!info->has_style("command", "python")) GTEST_SKIP();
     // execute python function from file
-    if (!verbose) ::testing::internal::CaptureStdout();
-    lmp->input->one("python printnum file ${input_dir}/func.py");
-    if (!verbose) ::testing::internal::GetCapturedStdout();
-    ::testing::internal::CaptureStdout();
-    lmp->input->one("python printnum invoke");
-    std::string output = ::testing::internal::GetCapturedStdout();
-    if (verbose) std::cout << output;
-    ASSERT_THAT(output, MatchesRegex("python.*2.25.*"));
+    HIDE_OUTPUT([&] {
+        command("python printnum file ${input_dir}/func.py");
+    });
 
+    auto output = CAPTURE_OUTPUT([&]() {
+        command("python printnum invoke");
+    });
+    ASSERT_THAT(output, HasSubstr("2.25\n"));
+}
+
+TEST_F(PythonPackageTest, InvokeFunctionPassInt)
+{
+    // execute python function, passing integer as argument
+    HIDE_OUTPUT([&] {
+        command("variable sq python square");
+        command("python square input 1 2 format ii return v_sq file ${input_dir}/func.py");
+        command("python square invoke");
+    });
+
+    ASSERT_EQ(get_variable_value("sq"), 4.0);
+}
+
+TEST_F(PythonPackageTest, InvokeFunctionPassFloat)
+{
+    // execute python function, passing float as argument
+    HIDE_OUTPUT([&] {
+        command("variable sq python square");
+        command("python square input 1 2.5 format ff return v_sq file ${input_dir}/func.py");
+    });
+
+    ASSERT_EQ(get_variable_value("sq"), 6.25);
+}
+
+TEST_F(PythonPackageTest, InvokeFunctionPassString)
+{
+    // execute python function, passing string as argument
+    HIDE_OUTPUT([&] {
+        command("variable val python bool_to_val");
+        command("python bool_to_val input 1 \"true\" format sf return v_val file ${input_dir}/func.py");
+    });
+
+    ASSERT_EQ(get_variable_value("val"), 1.0);
+}
+
+TEST_F(PythonPackageTest, InvokeFunctionPassStringVariable)
+{
+    // execute python function, passing string variable as argument
+    HIDE_OUTPUT([&] {
+        command("variable val python bool_to_val");
+        command("python bool_to_val input 1 v_str format sf return v_val file ${input_dir}/func.py");
+    });
+
+    HIDE_OUTPUT([&] {
+        command("variable str string \"true\"");
+    });
+
+    ASSERT_EQ(get_variable_value("val"), 1.0);
+
+    HIDE_OUTPUT([&] {
+        command("variable str string \"false\"");
+    });
+
+    ASSERT_EQ(get_variable_value("val"), 0.0);
+}
+
+TEST_F(PythonPackageTest, InvokeStringFunction)
+{
+    // execute python function, passing string variable as argument
+    HIDE_OUTPUT([&] {
+        command("variable str python val_to_bool");
+        command("python val_to_bool input 1 v_val format is return v_str file ${input_dir}/func.py");
+    });
+
+    HIDE_OUTPUT([&] {
+        command("variable val equal 0");
+    });
+
+    ASSERT_THAT(get_variable_string("str"), StrEq("False"));
+
+    HIDE_OUTPUT([&] {
+        command("variable val equal 1");
+    });
+
+    ASSERT_THAT(get_variable_string("str"), StrEq("True"));
+}
+
+TEST_F(PythonPackageTest, InvokeLongStringFunction)
+{
+    // execute python function, passing string variable as argument
+    HIDE_OUTPUT([&] {
+        command("variable str python longstr");
+        command("python longstr format s length 72 return v_str file ${input_dir}/func.py");
+    });
+
+    ASSERT_THAT(get_variable_string("str"), StrEq(LOREM_IPSUM));
+}
+
+TEST_F(PythonPackageTest, InvokeOtherFunctionFromFile)
+{
     // execute another python function from same file
-    if (!verbose) ::testing::internal::CaptureStdout();
-    lmp->input->one("python printtxt exists");
-    if (!verbose) ::testing::internal::GetCapturedStdout();
-    ::testing::internal::CaptureStdout();
-    lmp->input->one("python printtxt invoke");
-    output = ::testing::internal::GetCapturedStdout();
-    if (verbose) std::cout << output;
-    ASSERT_THAT(output, MatchesRegex("python.*sometext.*"));
+    HIDE_OUTPUT([&] {
+        command("python printnum file ${input_dir}/func.py");
+        command("python printtxt exists");
+    });
 
+    auto output = CAPTURE_OUTPUT([&] {
+        command("python printtxt invoke");
+    });
+    ASSERT_THAT(output, HasSubstr("sometext\n"));
+}
+
+TEST_F(PythonPackageTest, InvokeFunctionThatUsesLAMMPSModule)
+{
     // execute python function that uses the LAMMPS python module
-    if (!verbose) ::testing::internal::CaptureStdout();
-    lmp->input->one("variable idx equal 2.25");
-    lmp->input->one("python getidxvar input 1 SELF format p exists");
-    if (!verbose) ::testing::internal::GetCapturedStdout();
-    ::testing::internal::CaptureStdout();
-    lmp->input->one("python getidxvar invoke");
-    output = ::testing::internal::GetCapturedStdout();
-    if (verbose) std::cout << output;
-    ASSERT_THAT(output, MatchesRegex("python.*2.25.*"));
+    HIDE_OUTPUT([&] {
+        command("python printnum file ${input_dir}/func.py");
+        command("variable idx equal 2.25");
+        command("python getidxvar input 1 SELF format p exists");
+    });
+    auto output = CAPTURE_OUTPUT([&] {
+        command("python getidxvar invoke");
+    });
+    ASSERT_THAT(output, HasSubstr("2.25\n"));
 }
 
 TEST_F(PythonPackageTest, python_variable)
 {
-    if (!info->has_style("command", "python")) GTEST_SKIP();
-    if (!verbose) ::testing::internal::CaptureStdout();
-    lmp->input->one("variable sq python square");
-    lmp->input->one("variable val index 1.5");
-    lmp->input->one("python square input 1 v_val return v_sq format ff file ${input_dir}/func.py");
-    if (!verbose) ::testing::internal::GetCapturedStdout();
-    ::testing::internal::CaptureStdout();
-    lmp->input->one("print \"${sq}\"");
-    std::string output = ::testing::internal::GetCapturedStdout();
-    if (verbose) std::cout << output;
+    // define variable that evaluates a python function
+    HIDE_OUTPUT([&] {
+        command("variable sq python square");
+        command("variable val index 1.5");
+        command("python square input 1 v_val return v_sq format ff file ${input_dir}/func.py");
+    });
+    std::string output = CAPTURE_OUTPUT([&] {
+        command("print \"${sq}\"");
+    });
     ASSERT_THAT(output, MatchesRegex("print.*2.25.*"));
 }
 
+TEST_F(PythonPackageTest, InlineFunction)
+{
+    // define variable that evaluates a python function
+    HIDE_OUTPUT([&] {
+        command("variable fact python factorial");
+        command("python factorial input 1 v_n return v_fact format ii here \"\"\"\n"
+                "def factorial(n):\n"
+                "  if n == 0 or n == 1: return 1\n"
+                "  return n*factorial(n-1)\n"
+                "\"\"\"");
+    });
+
+    HIDE_OUTPUT([&] {
+        command("variable n equal 1");
+    });
+
+    ASSERT_EQ(get_variable_value("fact"), 1.0);
+
+    HIDE_OUTPUT([&] {
+        command("variable n equal 2");
+    });
+
+    ASSERT_EQ(get_variable_value("fact"), 2.0);
+
+    HIDE_OUTPUT([&] {
+        command("variable n equal 3");
+    });
+
+    ASSERT_EQ(get_variable_value("fact"), 6.0);
+}
+
+TEST_F(PythonPackageTest, RunSource)
+{
+    // execute python script from file
+    auto output = CAPTURE_OUTPUT([&] {
+        command("python xyz source ${input_dir}/run.py");
+    });
+
+    ASSERT_THAT(output, HasSubstr(LOREM_IPSUM));
+}
+
+TEST_F(PythonPackageTest, RunSourceInline)
+{
+    // execute inline python script
+    auto output = CAPTURE_OUTPUT([&] {
+        command("python xyz source \"\"\"\n"
+                "from __future__ import print_function\n"
+                "print(2+2)\n"
+                "\"\"\""
+        );
+    });
+
+    ASSERT_THAT(output, HasSubstr("4"));
+}
+
+TEST_F(FixPythonInvokeTest, end_of_step)
+{
+    HIDE_OUTPUT([&] {
+        command("python end_of_step_callback here \"\"\"\n"
+                "from __future__ import print_function\n"
+                "def end_of_step_callback(ptr):\n"
+                "    print(\"PYTHON_END_OF_STEP\")\n"
+                "\"\"\"");
+        command("fix eos all python/invoke 10 end_of_step end_of_step_callback");
+    });
+
+    auto output = CAPTURE_OUTPUT([&] {
+        command("run 50");
+    });
+
+    auto lines = utils::split_lines(output);
+    int count = 0;
+
+    for(auto & line : lines) {
+        if (line == "PYTHON_END_OF_STEP") ++count;
+    }
+
+    ASSERT_EQ(count, 5);
+}
+
+TEST_F(FixPythonInvokeTest, post_force)
+{
+    HIDE_OUTPUT([&] {
+        command("python post_force_callback here \"\"\"\n"
+                "from __future__ import print_function\n"
+                "def post_force_callback(ptr, vflag):\n"
+                "    print(\"PYTHON_POST_FORCE\")\n"
+                "\"\"\"");
+        command("fix pf all python/invoke 10 post_force post_force_callback");
+    });
+
+    auto output = CAPTURE_OUTPUT([&] {
+        command("run 50");
+    });
+
+    auto lines = utils::split_lines(output);
+    int count = 0;
+
+    for(auto & line : lines) {
+        if (line == "PYTHON_POST_FORCE") ++count;
+    }
+
+    ASSERT_EQ(count, 5);
+}
+
 } // namespace LAMMPS_NS
 
 int main(int argc, char **argv)
@@ -132,7 +330,7 @@ int main(int argc, char **argv)
 
     // handle arguments passed via environment variable
     if (const char *var = getenv("TEST_ARGS")) {
-        std::vector<std::string> env = split_words(var);
+        auto env = split_words(var);
         for (auto arg : env) {
             if (arg == "-v") {
                 verbose = true;
diff --git a/unittest/testing/systems/melt.h b/unittest/testing/systems/melt.h
index 4189a6b771..5ac92f562b 100644
--- a/unittest/testing/systems/melt.h
+++ b/unittest/testing/systems/melt.h
@@ -19,23 +19,25 @@ class MeltTest : public LAMMPSTest {
 protected:
     virtual void InitSystem() override
     {
-        command("units           lj");
-        command("atom_style      atomic");
-        command("atom_modify     map yes");
+        HIDE_OUTPUT([&] {
+            command("units           lj");
+            command("atom_style      atomic");
+            command("atom_modify     map yes");
 
-        command("lattice         fcc 0.8442");
-        command("region          box block 0 2 0 2 0 2");
-        command("create_box      1 box");
-        command("create_atoms    1 box");
-        command("mass            1 1.0");
+            command("lattice         fcc 0.8442");
+            command("region          box block 0 2 0 2 0 2");
+            command("create_box      1 box");
+            command("create_atoms    1 box");
+            command("mass            1 1.0");
 
-        command("velocity        all create 3.0 87287");
+            command("velocity        all create 3.0 87287");
 
-        command("pair_style      lj/cut 2.5");
-        command("pair_coeff      1 1 1.0 1.0 2.5");
+            command("pair_style      lj/cut 2.5");
+            command("pair_coeff      1 1 1.0 1.0 2.5");
 
-        command("neighbor        0.3 bin");
-        command("neigh_modify    every 20 delay 0 check no");
+            command("neighbor        0.3 bin");
+            command("neigh_modify    every 20 delay 0 check no");
+        });
     }
 };
 
diff --git a/unittest/utils/test_utils.cpp b/unittest/utils/test_utils.cpp
index 56a88e330a..81dda11615 100644
--- a/unittest/utils/test_utils.cpp
+++ b/unittest/utils/test_utils.cpp
@@ -113,7 +113,7 @@ TEST(Utils, count_words_with_extra_spaces)
 
 TEST(Utils, split_words_simple)
 {
-    std::vector<std::string> list = utils::split_words("one two three");
+    auto list = utils::split_words("one two three");
     ASSERT_EQ(list.size(), 3);
     ASSERT_THAT(list[0], StrEq("one"));
     ASSERT_THAT(list[1], StrEq("two"));
@@ -122,7 +122,7 @@ TEST(Utils, split_words_simple)
 
 TEST(Utils, split_words_quoted)
 {
-    std::vector<std::string> list = utils::split_words("one 'two' \"three\"");
+    auto list = utils::split_words("one 'two' \"three\"");
     ASSERT_EQ(list.size(), 3);
     ASSERT_THAT(list[0], StrEq("one"));
     ASSERT_THAT(list[1], StrEq("two"));
@@ -131,7 +131,7 @@ TEST(Utils, split_words_quoted)
 
 TEST(Utils, split_words_escaped)
 {
-    std::vector<std::string> list = utils::split_words("1\\' '\"two\"' 3\\\"");
+    auto list = utils::split_words("1\\' '\"two\"' 3\\\"");
     ASSERT_EQ(list.size(), 3);
     ASSERT_THAT(list[0], StrEq("1\\'"));
     ASSERT_THAT(list[1], StrEq("\"two\""));
@@ -140,13 +140,22 @@ TEST(Utils, split_words_escaped)
 
 TEST(Utils, split_words_quote_in_quoted)
 {
-    std::vector<std::string> list = utils::split_words("one 't\\'wo' \"th\\\"ree\"");
+    auto list = utils::split_words("one 't\\'wo' \"th\\\"ree\"");
     ASSERT_EQ(list.size(), 3);
     ASSERT_THAT(list[0], StrEq("one"));
     ASSERT_THAT(list[1], StrEq("t\\'wo"));
     ASSERT_THAT(list[2], StrEq("th\\\"ree"));
 }
 
+TEST(Utils, split_lines)
+{
+    auto list = utils::split_lines(" line 1\nline 2 \n line 3 \n");
+    ASSERT_EQ(list.size(), 3);
+    ASSERT_THAT(list[0], StrEq(" line 1"));
+    ASSERT_THAT(list[1], StrEq("line 2 "));
+    ASSERT_THAT(list[2], StrEq(" line 3 "));
+}
+
 TEST(Utils, valid_integer1)
 {
     ASSERT_TRUE(utils::is_integer("10"));