diff --git a/.github/workflows/push_pull.yml b/.github/workflows/push_pull.yml
index 3d8f0376700..26dfebcad9f 100644
--- a/.github/workflows/push_pull.yml
+++ b/.github/workflows/push_pull.yml
@@ -16,7 +16,7 @@ jobs:
       - name: Setup Python environment
         uses: actions/setup-python@v2
         with:
-          python-version: '3.7'
+          python-version: '3.8'
       - name: Check without sanitizer
         uses: ./.github/actions/build_and_check
         with:
@@ -33,7 +33,7 @@ jobs:
       - name: Setup Python environment
         uses: actions/setup-python@v2
         with:
-          python-version: '3.7'
+          python-version: '3.8'
       - name: Check with sanitizer
         uses: ./.github/actions/build_and_check
         with:
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index c8a54090ae1..2fec11330c3 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -21,7 +21,9 @@ stages:
   timeout: 40m
   interruptible: false
   tags:
-    - linux
+    - docker
+    - espresso
+    - no-cuda
 
 variables:
   GIT_SUBMODULE_STRATEGY: recursive
@@ -44,7 +46,8 @@ style:
     - sh maintainer/CI/fix_style.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
   variables:
     GIT_SUBMODULE_STRATEGY: none
   artifacts:
@@ -65,7 +68,8 @@ style_doxygen:
     - sh ../maintainer/CI/dox_warnings.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
 
 ### Builds without CUDA
 
@@ -85,7 +89,8 @@ default:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
 
 maxset:
   <<: *global_job_definition
@@ -105,7 +110,9 @@ maxset:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
+    - numa
 
 no_rotation:
   <<: *global_job_definition
@@ -122,7 +129,9 @@ no_rotation:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
+    - numa
 
 ubuntu:wo-dependencies:
   <<: *global_job_definition
@@ -138,7 +147,8 @@ ubuntu:wo-dependencies:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
 
 ### Builds with different distributions
 
@@ -155,7 +165,8 @@ debian:10:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
 
 fedora:34:
   <<: *global_job_definition
@@ -170,7 +181,8 @@ fedora:34:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
 
 ### Builds with CUDA
 
@@ -195,8 +207,9 @@ clang-sanitizer:
   timeout: 2h
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 fast_math:
   <<: *global_job_definition
@@ -215,7 +228,7 @@ fast_math:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
   when: manual
 
@@ -238,8 +251,9 @@ cuda11-maxset:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 cuda10-maxset:
   <<: *global_job_definition
@@ -263,8 +277,9 @@ cuda10-maxset:
     expire_in: 1 week
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 tutorials-samples-maxset:
   <<: *global_job_definition
@@ -287,7 +302,7 @@ tutorials-samples-maxset:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
 
 tutorials-samples-default:
@@ -310,7 +325,7 @@ tutorials-samples-default:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
   only:
     - schedules
@@ -336,7 +351,7 @@ tutorials-samples-empty:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
   only:
     - schedules
@@ -362,7 +377,8 @@ tutorials-samples-no-gpu:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
   only:
     - schedules
 
@@ -397,7 +413,7 @@ installation:
     - make -j2 check_tutorials
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
   when: manual
 
@@ -417,8 +433,9 @@ empty:
     - bash maintainer/CI/build_cmake.sh
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 check_sphinx:
   <<: *global_job_definition
@@ -440,8 +457,9 @@ check_sphinx:
     expire_in: 1 week
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 run_tutorials:
   <<: *global_job_definition
@@ -464,8 +482,9 @@ run_tutorials:
     expire_in: 1 week
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
   only:
     - schedules
 
@@ -486,7 +505,9 @@ run_doxygen:
     expire_in: 1 week
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
+    - numa
 
 maxset_no_gpu:
   <<: *global_job_definition
@@ -500,7 +521,9 @@ maxset_no_gpu:
     - make -t && make check
   tags:
     - docker
-    - linux
+    - espresso
+    - no-cuda
+    - numa
 
 maxset_3_cores:
   <<: *global_job_definition
@@ -514,8 +537,9 @@ maxset_3_cores:
     - make -t && make check_unit_tests && make check_python_parallel_odd
   tags:
     - docker
-    - linux
+    - espresso
     - cuda
+    - numa
 
 status_success:
   <<: *notification_job_definition
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9c9a0286b9..7a47427941e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -170,7 +170,7 @@ if(WITH_CUDA)
   endif()
 endif(WITH_CUDA)
 
-find_package(PythonInterp 3.7 REQUIRED)
+find_package(PythonInterp 3.8 REQUIRED)
 
 if(WITH_PYTHON)
   find_package(Cython 0.29.14 REQUIRED)
diff --git a/doc/sphinx/inter_bonded.rst b/doc/sphinx/inter_bonded.rst
index 465fb40f046..cf9b921bcc6 100644
--- a/doc/sphinx/inter_bonded.rst
+++ b/doc/sphinx/inter_bonded.rst
@@ -399,6 +399,9 @@ angle between the planes defined by the particle triples :math:`p_1`,
 Together with appropriate Lennard-Jones interactions, this potential can
 mimic a large number of atomic torsion potentials.
 
+Note that there is a singularity in the forces, but not in the energy, when
+:math:`\phi = 0` and :math:`\phi = \pi`.
+
 
 Tabulated dihedral potential
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -414,6 +417,8 @@ A tabulated dihedral interaction can be instantiated via
 The energy and force tables must be sampled from :math:`0` to :math:`2\pi`.
 For details of the interpolation, see :ref:`Tabulated interaction`.
 
+Note that there is a singularity in the forces, but not in the energy, when
+:math:`\phi = 0` and :math:`\phi = \pi`.
 
 .. _Immersed Boundary Method interactions:
 
diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst
index 038520c1038..62dc61fe2a7 100644
--- a/doc/sphinx/system_setup.rst
+++ b/doc/sphinx/system_setup.rst
@@ -158,7 +158,7 @@ Lees--Edwards boundary conditions
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Lees--Edwards boundary conditions (LEbc) are special periodic boundary
-conditions to simulation systems under shear stress :cite:`lees72a`.
+conditions to simulate systems under shear stress :cite:`lees72a`.
 Periodic images of particles across the shear boundary appear with a
 time-dependent position offset. When a particle crosses the shear boundary,
 it appears to the opposite side of the simulation box with a position offset
@@ -227,13 +227,13 @@ The properties of the cell system can be accessed by
 
 * :py:attr:`~espressomd.cellsystem.CellSystem.node_grid`
 
-  (int[3]) 3D node grid for real space domain decomposition (optional, if
+  3D node grid for real space domain decomposition (optional, if
   unset an optimal set is chosen automatically). The domain decomposition
   can be visualized with :file:`samples/visualization_cellsystem.py`.
 
 * :py:attr:`~espressomd.cellsystem.CellSystem.skin`
 
-  (float) Skin for the Verlet list. This value has to be set, otherwise the simulation will not start.
+  Skin for the Verlet list. This value has to be set, otherwise the simulation will not start.
 
 Details about the cell system can be obtained by :meth:`espressomd.system.System.cell_system.get_state() <espressomd.cellsystem.CellSystem.get_state>`:
 
diff --git a/maintainer/CI/build_cmake.sh b/maintainer/CI/build_cmake.sh
index 91b44886c75..29a553b8f54 100755
--- a/maintainer/CI/build_cmake.sh
+++ b/maintainer/CI/build_cmake.sh
@@ -72,6 +72,12 @@ set_default_value() {
     fi
 }
 
+# the number of available processors depends on the CI runner
+if grep -q "i7-3820" /proc/cpuinfo; then
+   ci_procs=2
+else
+   ci_procs=4
+fi
 
 # handle environment variables
 set_default_value srcdir "$(pwd)"
@@ -82,7 +88,7 @@ set_default_value with_ubsan false
 set_default_value with_asan false
 set_default_value with_static_analysis false
 set_default_value myconfig "default"
-set_default_value build_procs 2
+set_default_value build_procs ${ci_procs}
 set_default_value check_procs ${build_procs}
 set_default_value check_odd_only false
 set_default_value check_gpu_only false
diff --git a/maintainer/CI/deploy_tutorials.py b/maintainer/CI/deploy_tutorials.py
index 39494a2f7cf..adae4b36dc1 100755
--- a/maintainer/CI/deploy_tutorials.py
+++ b/maintainer/CI/deploy_tutorials.py
@@ -20,20 +20,18 @@
 
 """List all tutorial files to deploy (PDF, HTML and figures)"""
 
-import os
-import glob
+import pathlib
 import lxml.html
 
-deploy_list = glob.glob('*/*.pdf')
-for filepath in glob.glob('*/*.html'):
+deploy_list = list(pathlib.Path().glob('*/*.pdf'))
+for filepath in pathlib.Path().glob('*/*.html'):
     deploy_list.append(filepath)
     # extract all figures
-    dirname = os.path.dirname(filepath)
     with open(filepath, encoding='utf-8') as f:
         html = lxml.html.parse(f)
     figures = filter(lambda src: not src.startswith('data:image'),
                      html.xpath('//img/@src'))
-    deploy_list += list(map(lambda src: os.path.join(dirname, src), figures))
+    deploy_list += list(map(lambda src: filepath.parent / src, figures))
 
 with open('deploy_list.txt', 'w') as f:
-    f.write('\n'.join(deploy_list))
+    f.write('\n'.join(map(str, deploy_list)))
diff --git a/maintainer/CI/dox_warnings.py b/maintainer/CI/dox_warnings.py
index 897b37f8512..035bf1d4f9d 100755
--- a/maintainer/CI/dox_warnings.py
+++ b/maintainer/CI/dox_warnings.py
@@ -18,7 +18,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 import re
-import os
+import pathlib
 
 # collect list of Doxygen warnings
 with open('doc/doxygen/warnings.log') as f:
@@ -38,7 +38,7 @@
     for line in content.strip().split('\n'):
         m = re.search(r'^(.+):(\d+):[\s\*]*([@\\]t?param).*\s(\S+)\s*$', line)
         filepath, lineno, paramtype, varname = m.groups()
-        ext = os.path.splitext(filepath)[1]
+        ext = pathlib.Path(filepath).suffix
         if ext.lower() not in source_code_ext:
             continue
         warning = (f'argument \'{varname}\' of {paramtype} has no description,'
diff --git a/maintainer/CI/jupyter_warnings.py b/maintainer/CI/jupyter_warnings.py
index 19bda8e477e..ce0ec832711 100755
--- a/maintainer/CI/jupyter_warnings.py
+++ b/maintainer/CI/jupyter_warnings.py
@@ -22,9 +22,8 @@
 pointing to sections of the notebooks or to the online Sphinx documentation.
 """
 
-import os
 import sys
-import glob
+import pathlib
 
 import lxml.etree
 import nbformat
@@ -61,7 +60,7 @@ def detect_invalid_urls(nb, sphinx_root='.'):
     root = lxml.etree.fromstring(html_string, parser=html_parser)
     # process all links
     espressomd_website_root = 'https://espressomd.github.io/doc/'
-    sphinx_html_root = os.path.join(sphinx_root, 'doc', 'sphinx', 'html')
+    sphinx_html_root = pathlib.Path(sphinx_root) / 'doc' / 'sphinx' / 'html'
     broken_links = []
     for link in root.xpath('//a'):
         url = link.attrib.get('href', '')
@@ -75,15 +74,15 @@ def detect_invalid_urls(nb, sphinx_root='.'):
                 url = url.split('?', 1)[0]
             # check file exists
             basename = url.split(espressomd_website_root, 1)[1]
-            filepath = os.path.join(sphinx_html_root, basename)
-            if not os.path.isfile(filepath):
+            filepath = sphinx_html_root / basename
+            if not filepath.is_file():
                 broken_links.append(f'{url} does not exist')
                 continue
             # check anchor exists
             if anchor is not None:
                 if filepath not in sphinx_docs:
                     sphinx_docs[filepath] = lxml.etree.parse(
-                        filepath, parser=html_parser)
+                        str(filepath), parser=html_parser)
                 doc = sphinx_docs[filepath]
                 nodes = doc.xpath(f'//*[@id="{anchor}"]')
                 if not nodes:
@@ -99,14 +98,13 @@ def detect_invalid_urls(nb, sphinx_root='.'):
 
 if __name__ == '__main__':
     error_code = 0
-    for nb_filepath in sorted(glob.glob('doc/tutorials/*/*.ipynb')):
+    for nb_filepath in sorted(pathlib.Path().glob('doc/tutorials/*/*.ipynb')):
         with open(nb_filepath, encoding='utf-8') as f:
             nb = nbformat.read(f, as_version=4)
         issues = detect_invalid_urls(nb)
         if issues:
             error_code = 1
-            basename = os.path.basename(nb_filepath)
-            print(f'In notebook {basename}:', file=sys.stderr)
+            print(f'In notebook {nb_filepath.name}:', file=sys.stderr)
             for issue in issues:
                 print(f'* {issue}', file=sys.stderr)
     if not error_code:
diff --git a/maintainer/benchmarks/benchmarks.py b/maintainer/benchmarks/benchmarks.py
index 2e37a6a46d4..dfb3fc81f33 100644
--- a/maintainer/benchmarks/benchmarks.py
+++ b/maintainer/benchmarks/benchmarks.py
@@ -16,9 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-import os
 import sys
 import time
+import pathlib
 import numpy as np
 
 
@@ -127,12 +127,12 @@ def write_report(filepath, n_proc, timings, n_steps, label=''):
         Label to distinguish e.g. MD from MC or LB steps.
 
     '''
-    script = os.path.basename(sys.argv[0])
+    script = pathlib.Path(sys.argv[0]).name
     cmd = " ".join(x for x in sys.argv[1:] if not x.startswith("--output"))
     avg, ci = get_average_time(timings)
     header = '"script","arguments","cores","mean","ci","nsteps","duration","label"\n'
     report = f'"{script}","{cmd}",{n_proc},{avg:.3e},{ci:.3e},{n_steps},{np.sum(timings):.1f},"{label}"\n'
-    if os.path.isfile(filepath):
+    if pathlib.Path(filepath).is_file():
         header = ''
     with open(filepath, "a") as f:
         f.write(header + report)
diff --git a/src/core/bonded_interactions/bonded_tab.hpp b/src/core/bonded_interactions/bonded_tab.hpp
index 4863697bd1d..99498d5c4bb 100644
--- a/src/core/bonded_interactions/bonded_tab.hpp
+++ b/src/core/bonded_interactions/bonded_tab.hpp
@@ -203,7 +203,8 @@ inline double TabulatedAngleBond::energy(Utils::Vector3d const &r_mid,
 }
 
 /** Compute the four-body dihedral interaction force.
- *  This function is not tested yet.
+ *  The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$
+ *  (see @cite swope92a page 592).
  *
  *  @param[in]  r1        Position of the first particle.
  *  @param[in]  r2        Position of the second particle.
@@ -220,18 +221,18 @@ TabulatedDihedralBond::forces(Utils::Vector3d const &r1,
   /* vectors for dihedral angle calculation */
   Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34;
   double l_v12Xv23, l_v23Xv34;
-  /* dihedral angle, cosine of the dihedral angle, cosine of the bond angles */
+  /* dihedral angle, cosine of the dihedral angle */
   double phi, cos_phi;
 
   /* dihedral angle */
-  calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
-                      v23Xv34, &l_v23Xv34, &cos_phi, &phi);
+  auto const angle_is_undefined =
+      calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23,
+                          v23Xv34, l_v23Xv34, cos_phi, phi);
   /* dihedral angle not defined - force zero */
-  if (phi == -1.0) {
+  if (angle_is_undefined) {
     return {};
   }
 
-  /* calculate force components (directions) */
   auto const f1 = (v23Xv34 - cos_phi * v12Xv23) / l_v12Xv23;
   auto const f4 = (v12Xv23 - cos_phi * v23Xv34) / l_v23Xv34;
 
@@ -252,7 +253,7 @@ TabulatedDihedralBond::forces(Utils::Vector3d const &r1,
 }
 
 /** Compute the four-body dihedral interaction energy.
- *  This function is not tested yet.
+ *  The energy doesn't have any singularity if the angle phi is well-defined.
  *
  *  @param[in]  r1        Position of the first particle.
  *  @param[in]  r2        Position of the second particle.
@@ -267,8 +268,14 @@ inline boost::optional<double> TabulatedDihedralBond::energy(
   double l_v12Xv23, l_v23Xv34;
   /* dihedral angle, cosine of the dihedral angle */
   double phi, cos_phi;
-  calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
-                      v23Xv34, &l_v23Xv34, &cos_phi, &phi);
+  auto const angle_is_undefined =
+      calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23,
+                          v23Xv34, l_v23Xv34, cos_phi, phi);
+  /* dihedral angle not defined - energy zero */
+  if (angle_is_undefined) {
+    return {};
+  }
+
   return pot->energy(phi);
 }
 
diff --git a/src/core/bonded_interactions/dihedral.hpp b/src/core/bonded_interactions/dihedral.hpp
index f02ea9c26d2..0991f8304ee 100644
--- a/src/core/bonded_interactions/dihedral.hpp
+++ b/src/core/bonded_interactions/dihedral.hpp
@@ -83,7 +83,7 @@ struct DihedralBond {
  * specified by the particle triples (p1,p2,p3) and (p2,p3,p4).
  * Vectors a, b and c are the bond vectors between consecutive particles.
  * If the a,b or b,c are parallel the dihedral angle is not defined in which
- * case the routine returns phi=-1. Calling functions should check for that
+ * case the function returns true. Calling functions should check for that.
  *
  * @param[in]  r1 , r2 , r3 , r4 Positions of the particles forming the dihedral
  * @param[out] a Vector from @p p1 to @p p2
@@ -94,14 +94,15 @@ struct DihedralBond {
  * @param[out] bXc Vector product of b and c
  * @param[out] l_bXc |bXc|
  * @param[out] cosphi Cosine of the dihedral angle
- * @param[out] phi Dihedral angle
+ * @param[out] phi Dihedral angle in the range [0, pi]
+ * @return Whether the angle is undefined.
  */
-inline void
+inline bool
 calc_dihedral_angle(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
                     Utils::Vector3d const &r3, Utils::Vector3d const &r4,
                     Utils::Vector3d &a, Utils::Vector3d &b, Utils::Vector3d &c,
-                    Utils::Vector3d &aXb, double *l_aXb, Utils::Vector3d &bXc,
-                    double *l_bXc, double *cosphi, double *phi) {
+                    Utils::Vector3d &aXb, double &l_aXb, Utils::Vector3d &bXc,
+                    double &l_bXc, double &cosphi, double &phi) {
   a = box_geo.get_mi_vector(r2, r1);
   b = box_geo.get_mi_vector(r3, r2);
   c = box_geo.get_mi_vector(r4, r3);
@@ -111,31 +112,34 @@ calc_dihedral_angle(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
   bXc = vector_product(b, c);
 
   /* calculate the unit vectors */
-  *l_aXb = aXb.norm();
-  *l_bXc = bXc.norm();
+  l_aXb = aXb.norm();
+  l_bXc = bXc.norm();
 
   /* catch case of undefined dihedral angle */
-  if (*l_aXb <= TINY_LENGTH_VALUE || *l_bXc <= TINY_LENGTH_VALUE) {
-    *phi = -1.0;
-    *cosphi = 0;
-    return;
+  if (l_aXb <= TINY_LENGTH_VALUE || l_bXc <= TINY_LENGTH_VALUE) {
+    phi = -1.0;
+    cosphi = 0.0;
+    return true;
   }
 
-  aXb /= *l_aXb;
-  bXc /= *l_bXc;
+  aXb /= l_aXb;
+  bXc /= l_bXc;
 
-  *cosphi = aXb * bXc;
+  cosphi = aXb * bXc;
 
-  if (fabs(fabs(*cosphi) - 1) < TINY_SIN_VALUE)
-    *cosphi = std::round(*cosphi);
+  if (fabs(fabs(cosphi) - 1) < TINY_SIN_VALUE)
+    cosphi = std::round(cosphi);
 
   /* Calculate dihedral angle */
-  *phi = acos(*cosphi);
+  phi = acos(cosphi);
   if ((aXb * c) < 0.0)
-    *phi = (2.0 * Utils::pi()) - *phi;
+    phi = (2.0 * Utils::pi()) - phi;
+  return false;
 }
 
 /** Compute the four-body dihedral interaction force.
+ *  The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$
+ *  (see @cite swope92a page 592).
  *
  *  @param[in]  r1        Position of the first particle.
  *  @param[in]  r2        Position of the second particle.
@@ -152,20 +156,19 @@ DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
   Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34;
   double l_v12Xv23, l_v23Xv34;
   /* dihedral angle, cosine of the dihedral angle */
-  double phi, cosphi, sinmphi_sinphi;
-  /* force factors */
-  double fac;
+  double phi, cos_phi, sin_mphi_over_sin_phi;
 
   /* dihedral angle */
-  calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
-                      v23Xv34, &l_v23Xv34, &cosphi, &phi);
+  auto const angle_is_undefined =
+      calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23,
+                          v23Xv34, l_v23Xv34, cos_phi, phi);
   /* dihedral angle not defined - force zero */
-  if (phi == -1.0) {
+  if (angle_is_undefined) {
     return {};
   }
 
-  auto const f1 = (v23Xv34 - cosphi * v12Xv23) / l_v12Xv23;
-  auto const f4 = (v12Xv23 - cosphi * v23Xv34) / l_v23Xv34;
+  auto const f1 = (v23Xv34 - cos_phi * v12Xv23) / l_v12Xv23;
+  auto const f4 = (v12Xv23 - cos_phi * v23Xv34) / l_v23Xv34;
 
   auto const v23Xf1 = vector_product(v23, f1);
   auto const v23Xf4 = vector_product(v23, f4);
@@ -173,19 +176,17 @@ DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
   auto const v12Xf1 = vector_product(v12, f1);
 
   /* calculate force magnitude */
-  fac = -bend * mult;
+  auto fac = -bend * mult;
 
   if (fabs(sin(phi)) < TINY_SIN_VALUE) {
-    /*(comes from taking the first term of the MacLaurin expansion of
-      sin(n*phi - phi0) and sin(phi) and then making the division).
-      The original code had a 2PI term in the cosine (cos(2PI - nPhi))
-      but I removed it because it wasn't doing anything. AnaVV*/
-    sinmphi_sinphi = mult * cos(mult * phi - phase) / cosphi;
+    /* comes from taking the first term of the MacLaurin expansion of
+     * sin(n * phi - phi0) and sin(phi) and then making the division */
+    sin_mphi_over_sin_phi = mult * cos(mult * phi - phase) / cos_phi;
   } else {
-    sinmphi_sinphi = sin(mult * phi - phase) / sin(phi);
+    sin_mphi_over_sin_phi = sin(mult * phi - phase) / sin(phi);
   }
 
-  fac *= sinmphi_sinphi;
+  fac *= sin_mphi_over_sin_phi;
 
   /* store dihedral forces */
   auto const force1 = fac * v23Xf1;
@@ -196,6 +197,7 @@ DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
 }
 
 /** Compute the four-body dihedral interaction energy.
+ *  The energy doesn't have any singularity if the angle phi is well-defined.
  *
  *  @param[in]  r1        Position of the first particle.
  *  @param[in]  r2        Position of the second particle.
@@ -210,12 +212,13 @@ DihedralBond::energy(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
   Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34;
   double l_v12Xv23, l_v23Xv34;
   /* dihedral angle, cosine of the dihedral angle */
-  double phi, cosphi;
+  double phi, cos_phi;
 
-  calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
-                      v23Xv34, &l_v23Xv34, &cosphi, &phi);
-  /* dihedral angle not defined - force zero */
-  if (phi == -1.0) {
+  auto const angle_is_undefined =
+      calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23,
+                          v23Xv34, l_v23Xv34, cos_phi, phi);
+  /* dihedral angle not defined - energy zero */
+  if (angle_is_undefined) {
     return {};
   }
 
diff --git a/src/core/lees_edwards/protocols.hpp b/src/core/lees_edwards/protocols.hpp
index 65a275497b4..29737168b40 100644
--- a/src/core/lees_edwards/protocols.hpp
+++ b/src/core/lees_edwards/protocols.hpp
@@ -37,7 +37,8 @@ struct Off {
 
 /** Lees-Edwards protocol for linear shearing */
 struct LinearShear {
-  LinearShear() : m_initial_pos_offset{0}, m_shear_velocity{0}, m_time_0{0} {}
+  LinearShear()
+      : m_initial_pos_offset{0.}, m_shear_velocity{0.}, m_time_0{0.} {}
   LinearShear(double initial_offset, double shear_velocity, double time_0)
       : m_initial_pos_offset{initial_offset},
         m_shear_velocity{shear_velocity}, m_time_0{time_0} {}
@@ -52,15 +53,20 @@ struct LinearShear {
 
 /** Lees-Edwards protocol for oscillatory shearing */
 struct OscillatoryShear {
-  OscillatoryShear() : m_amplitude{0}, m_omega{0}, m_time_0{0} {}
-  OscillatoryShear(double amplitude, double omega, double time_0)
-      : m_amplitude{amplitude}, m_omega{omega}, m_time_0{time_0} {}
+  OscillatoryShear()
+      : m_initial_pos_offset{0.}, m_amplitude{0.}, m_omega{0.}, m_time_0{0.} {}
+  OscillatoryShear(double initial_offset, double amplitude, double omega,
+                   double time_0)
+      : m_initial_pos_offset{initial_offset},
+        m_amplitude{amplitude}, m_omega{omega}, m_time_0{time_0} {}
   double pos_offset(double time) const {
-    return m_amplitude * std::sin(m_omega * (time - m_time_0));
+    return m_initial_pos_offset +
+           m_amplitude * std::sin(m_omega * (time - m_time_0));
   }
   double shear_velocity(double time) const {
     return m_omega * m_amplitude * std::cos(m_omega * (time - m_time_0));
   }
+  double m_initial_pos_offset;
   double m_amplitude;
   double m_omega;
   double m_time_0;
diff --git a/src/core/unit_tests/lees_edwards_test.cpp b/src/core/unit_tests/lees_edwards_test.cpp
index 0f9003f3eaf..8739059a363 100644
--- a/src/core/unit_tests/lees_edwards_test.cpp
+++ b/src/core/unit_tests/lees_edwards_test.cpp
@@ -132,10 +132,12 @@ BOOST_AUTO_TEST_CASE(protocol_lin) {
 
 BOOST_AUTO_TEST_CASE(protocol_osc) {
   auto const t0 = 1.2;
+  auto const x0 = 0.1;
   auto const a = 3.1;
   auto const o = 2.1;
-  auto osc = OscillatoryShear(a, o, t0);
-  BOOST_CHECK_CLOSE(get_pos_offset(3.3, osc), a * sin(o * (3.3 - t0)), tol);
+  auto osc = OscillatoryShear(x0, a, o, t0);
+  BOOST_CHECK_CLOSE(get_pos_offset(3.3, osc), x0 + a * sin(o * (3.3 - t0)),
+                    tol);
   BOOST_CHECK_CLOSE(get_shear_velocity(3.3, osc), a * o * cos(o * (3.3 - t0)),
                     tol);
 }
diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx
index 48cd2253111..d14b9b8cf8c 100644
--- a/src/python/espressomd/analyze.pyx
+++ b/src/python/espressomd/analyze.pyx
@@ -374,10 +374,7 @@ class Analysis:
 
         """
         self.check_topology(chain_start, number_of_chains, chain_length)
-        re = analyze.calc_re(
-            chain_start,
-            number_of_chains,
-            chain_length)
+        re = analyze.calc_re(chain_start, number_of_chains, chain_length)
         return np.array([re[0], re[1], re[2], re[3]])
 
     def calc_rg(self, chain_start=None, number_of_chains=None,
@@ -410,10 +407,7 @@ class Analysis:
 
         """
         self.check_topology(chain_start, number_of_chains, chain_length)
-        rg = analyze.calc_rg(
-            chain_start,
-            number_of_chains,
-            chain_length)
+        rg = analyze.calc_rg(chain_start, number_of_chains, chain_length)
         return np.array([rg[0], rg[1], rg[2], rg[3]])
 
     def calc_rh(self, chain_start=None, number_of_chains=None,
@@ -445,10 +439,7 @@ class Analysis:
         """
 
         self.check_topology(chain_start, number_of_chains, chain_length)
-        rh = analyze.calc_rh(
-            chain_start,
-            number_of_chains,
-            chain_length)
+        rh = analyze.calc_rh(chain_start, number_of_chains, chain_length)
         return np.array([rh[0], rh[1]])
 
     def check_topology(self, chain_start=None, number_of_chains=None,
diff --git a/src/python/espressomd/lees_edwards.py b/src/python/espressomd/lees_edwards.py
index 8b64df06059..b4e18a929e5 100644
--- a/src/python/espressomd/lees_edwards.py
+++ b/src/python/espressomd/lees_edwards.py
@@ -22,9 +22,23 @@
 @script_interface_register
 class LeesEdwards(ScriptInterfaceHelper):
 
-    """Interface to the Lees-Edwards boundary conditions.
+    """
+    Interface to the :ref:`Lees-Edwards boundary conditions`.
 
-       See documentation.
+    Attributes
+    ----------
+    protocol : :obj:`object`
+        Lees--Edwards protocol.
+    shear_velocity: :obj:`float`
+        Current shear velocity.
+    pos_offset : :obj:`float`
+        Current position offset
+    shear_direction : :obj:`int`
+        Shear direction: 0 for the *x*-axis, 1 for the *y*-axis,
+        2 for the *z*-axis.
+    shear_plane_normal : :obj:`int`
+        Shear plane normal: 0 for the *x*-axis, 1 for the *y*-axis,
+        2 for the *z*-axis.
 
     """
 
@@ -34,25 +48,21 @@ class LeesEdwards(ScriptInterfaceHelper):
 @script_interface_register
 class Off(ScriptInterfaceHelper):
 
-    """Lees-Edwards protocol resulting in un-shifted boundaries."""
+    """Lees--Edwards protocol resulting in un-shifted boundaries."""
     _so_name = "LeesEdwards::Off"
 
 
 @script_interface_register
 class LinearShear(ScriptInterfaceHelper):
 
-    """Lees-Edwards protocol for linear shear.
+    """Lees--Edwards protocol for linear shear.
 
     Parameters
     ----------
-    shear_direction
-       Cartesian coordinate of the shear direction (0=x, 1=y, 2=z)
-    shear_plane_normal
-       Cartesian coordinate of the shear plane normal
-    initial_pos_offset
-       Positional offset at the Lees-Edwards boundary at t=0
-    shear_velocity
-       Shear velocity (velocity jump) across the Lees-Edwards boundary
+    initial_pos_offset : :obj:`float`
+       Positional offset at the Lees--Edwards boundary at t=0.
+    shear_velocity : :obj:`float`
+       Shear velocity (velocity jump) across the Lees--Edwards boundary.
 
     """
     _so_name = "LeesEdwards::LinearShear"
@@ -61,21 +71,18 @@ class LinearShear(ScriptInterfaceHelper):
 @script_interface_register
 class OscillatoryShear(ScriptInterfaceHelper):
 
-    """Lees-Edwards protocol for oscillatory shear.
+    """Lees--Edwards protocol for oscillatory shear.
 
     Parameters
     ----------
-    shear_direction
-       Cartesian coordinate of the shear direction (0=x, 1=y, 2=z)
-    shear_plane_normal
-       Cartesian coordinate of the shear plane normal
-    amplitude
-       Maximum amplitude of the positional offset at the Lees-Edwards boundary
-    frequency
-       Frequency of the shear
-    time_0
-       Time offset of the oscillation
-
+    initial_pos_offset : :obj:`float`
+       Positional offset at the Lees--Edwards boundary at t=0.
+    amplitude : :obj:`float`
+       Maximum amplitude of the positional offset at the Lees--Edwards boundary.
+    omega : :obj:`float`
+       Radian frequency of the oscillation.
+    time_0 : :obj:`float`
+       Time offset of the oscillation.
 
     """
     _so_name = "LeesEdwards::OscillatoryShear"
diff --git a/src/script_interface/lees_edwards/OscillatoryShear.hpp b/src/script_interface/lees_edwards/OscillatoryShear.hpp
index 9129d3ad4b3..7bf9923b193 100644
--- a/src/script_interface/lees_edwards/OscillatoryShear.hpp
+++ b/src/script_interface/lees_edwards/OscillatoryShear.hpp
@@ -37,7 +37,10 @@ class OscillatoryShear : public Protocol {
       : m_protocol{std::make_shared<::LeesEdwards::ActiveProtocol>(
             ::LeesEdwards::OscillatoryShear())} {
     add_parameters(
-        {{"amplitude",
+        {{"initial_pos_offset",
+          boost::get<::LeesEdwards::OscillatoryShear>(*m_protocol)
+              .m_initial_pos_offset},
+         {"amplitude",
           boost::get<::LeesEdwards::OscillatoryShear>(*m_protocol).m_amplitude},
          {"omega",
           boost::get<::LeesEdwards::OscillatoryShear>(*m_protocol).m_omega},
diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py
index 655a6e7de48..4f2d0025692 100644
--- a/testsuite/python/analyze_chains.py
+++ b/testsuite/python/analyze_chains.py
@@ -133,6 +133,16 @@ def test_radii(self):
         self.system.box_l = self.system.box_l / 2.
         all_partcls.pos = old_pos
 
+    def test_exceptions(self):
+        err_msg = """particle with id 10 does not exist
+cannot perform analysis on the range chain_start=0, number_of_chains=2, chain_length=10
+please provide a contiguous range of particle ids"""
+        analysis = self.system.analysis
+        for method in (analysis.calc_re, analysis.calc_rg, analysis.calc_rh):
+            with self.assertRaisesRegex(ValueError, err_msg):
+                method(chain_start=0, number_of_chains=self.num_poly,
+                       chain_length=2 * self.num_mono)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/python/cluster_analysis.py b/testsuite/python/cluster_analysis.py
index 89945dfa1e9..aad912b4cfe 100644
--- a/testsuite/python/cluster_analysis.py
+++ b/testsuite/python/cluster_analysis.py
@@ -149,28 +149,26 @@ def test_single_cluster_analysis(self):
         # Unknown method should return None
         self.assertIsNone(c.call_method("unknown"))
 
-        # Fractal dimension calc require gsl
-        if not espressomd.has_features("GSL"):
-            print("Skipping fractal dimension tests due to missing GSL dependency")
-            return
-
-        # The fractal dimension of a line should be 1
-        dr = 0.
-        self.assertAlmostEqual(c.fractal_dimension(dr=dr)[0], 1, delta=0.05)
-
-        # The fractal dimension of a disk should be close to 2
-        self.system.part.clear()
-        center = np.array((0.1, .02, 0.15))
-        for _ in range(3000):
-            r_inv, phi = np.random.random(2) * np.array((0.2, 2 * np.pi))
-            r = 1 / r_inv
-            self.system.part.add(
-                pos=center + r * np.array((np.sin(phi), np.cos(phi), 0)))
-        self.cs.clear()
-        self.cs.run_for_all_pairs()
-        cid = self.cs.cluster_ids()[0]
-        df = self.cs.clusters[cid].fractal_dimension(dr=0.001)
-        self.assertAlmostEqual(df[0], 2, delta=0.08)
+        # Fractal dimension calculation requires gsl
+        if espressomd.has_features("GSL"):
+            # The fractal dimension of a line should be 1
+            dr = 0.
+            self.assertAlmostEqual(
+                c.fractal_dimension(dr=dr)[0], 1, delta=0.05)
+
+            # The fractal dimension of a disk should be close to 2
+            self.system.part.clear()
+            center = np.array((0.1, .02, 0.15))
+            for _ in range(3000):
+                r_inv, phi = np.random.random(2) * np.array((0.2, 2 * np.pi))
+                r = 1 / r_inv
+                self.system.part.add(
+                    pos=center + r * np.array((np.sin(phi), np.cos(phi), 0)))
+            self.cs.clear()
+            self.cs.run_for_all_pairs()
+            cid = self.cs.cluster_ids()[0]
+            df = self.cs.clusters[cid].fractal_dimension(dr=0.001)
+            self.assertAlmostEqual(df[0], 2, delta=0.08)
 
     def test_analysis_for_bonded_particles(self):
         self.set_two_clusters()
diff --git a/testsuite/python/coulomb_cloud_wall.py b/testsuite/python/coulomb_cloud_wall.py
index c661660f837..26d3f048c8a 100644
--- a/testsuite/python/coulomb_cloud_wall.py
+++ b/testsuite/python/coulomb_cloud_wall.py
@@ -36,8 +36,8 @@ class CoulombCloudWall(ut.TestCase):
     """
 
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
-    data = np.genfromtxt(tests_common.abspath(
-        "data/coulomb_cloud_wall_system.data"))
+    data = np.genfromtxt(tests_common.data_path(
+        "coulomb_cloud_wall_system.data"))
 
     tolerance = 1E-3
 
diff --git a/testsuite/python/coulomb_cloud_wall_duplicated.py b/testsuite/python/coulomb_cloud_wall_duplicated.py
index 79289a0f360..492f170dd91 100644
--- a/testsuite/python/coulomb_cloud_wall_duplicated.py
+++ b/testsuite/python/coulomb_cloud_wall_duplicated.py
@@ -35,7 +35,7 @@ class CoulombCloudWall(ut.TestCase):
     system.time_step = 0.01
     system.cell_system.skin = 0.4
     data = np.genfromtxt(
-        tests_common.abspath("data/coulomb_cloud_wall_duplicated_system.data"))
+        tests_common.data_path("coulomb_cloud_wall_duplicated_system.data"))
 
     tolerance = 1E-3
 
diff --git a/testsuite/python/coulomb_mixed_periodicity.py b/testsuite/python/coulomb_mixed_periodicity.py
index 706f6774690..730e42a0993 100644
--- a/testsuite/python/coulomb_mixed_periodicity.py
+++ b/testsuite/python/coulomb_mixed_periodicity.py
@@ -31,8 +31,8 @@ class CoulombMixedPeriodicity(ut.TestCase):
     """Test mixed periodicity electrostatics"""
 
     system = espressomd.System(box_l=[10, 10, 10])
-    data = np.genfromtxt(tests_common.abspath(
-        "data/coulomb_mixed_periodicity_system.data"))
+    data = np.genfromtxt(tests_common.data_path(
+        "coulomb_mixed_periodicity_system.data"))
 
     tolerance_force = 5E-4
     tolerance_energy = 1.8E-3
diff --git a/testsuite/python/coulomb_tuning.py b/testsuite/python/coulomb_tuning.py
index 737cb5c4d47..dbe70141ded 100644
--- a/testsuite/python/coulomb_tuning.py
+++ b/testsuite/python/coulomb_tuning.py
@@ -39,7 +39,7 @@ def setUp(self):
         self.system.time_step = 0.01
         self.system.cell_system.skin = 0.4
 
-        data = np.load(tests_common.abspath("data/coulomb_tuning_system.npz"))
+        data = np.load(tests_common.data_path("coulomb_tuning_system.npz"))
         self.forces = data['forces']
         self.system.part.add(pos=data['pos'], q=data['charges'])
 
diff --git a/testsuite/python/decorators.py b/testsuite/python/decorators.py
index 139dcfa6b9f..1907441cad0 100644
--- a/testsuite/python/decorators.py
+++ b/testsuite/python/decorators.py
@@ -54,6 +54,9 @@ def test_missing_modules(self):
         decorator = utx.skipIfMissingModules(['numpy___', 'scipy___'])
         args = self.get_skip_reason(decorator)
         self.assertEqual(args, (err_msg + 'modules numpy___, scipy___',))
+        decorator = utx.skipIfMissingModules(['espressomd'])
+        args = self.get_skip_reason(decorator)
+        self.assertIsNone(args)
 
     def test_missing_gpu(self):
         espressomd.gpu_available = lambda: False
diff --git a/testsuite/python/dipolar_direct_summation.py b/testsuite/python/dipolar_direct_summation.py
index acb688ef84d..f341261e504 100644
--- a/testsuite/python/dipolar_direct_summation.py
+++ b/testsuite/python/dipolar_direct_summation.py
@@ -19,15 +19,15 @@
 import espressomd
 import espressomd.magnetostatics
 import espressomd.magnetostatic_extensions
-import os
+import pathlib
 import numpy as np
 import unittest as ut
 import unittest_decorators as utx
 import tests_common
-OPEN_BOUNDARIES_REF_ENERGY = tests_common.abspath(
-    "data/dipolar_open_boundaries_energy.npy")
-OPEN_BOUNDARIES_REF_ARRAYS = tests_common.abspath(
-    "data/dipolar_open_boundaries_arrays.npy")
+OPEN_BOUNDARIES_REF_ENERGY = tests_common.data_path(
+    "dipolar_open_boundaries_energy.npy")
+OPEN_BOUNDARIES_REF_ARRAYS = tests_common.data_path(
+    "dipolar_open_boundaries_arrays.npy")
 
 
 @utx.skipIfMissingFeatures(["DIPOLES"])
@@ -114,12 +114,11 @@ def test_gen_reference_data(self):
         filepaths = ('dipolar_direct_summation_energy.npy',
                      'dipolar_direct_summation_arrays.npy')
         for filepath in filepaths:
-            if os.path.isfile(filepath):
-                os.remove(filepath)
+            pathlib.Path(filepath).unlink(missing_ok=True)
 
         self.gen_reference_data(filepaths[0], filepaths[1])
         for filepath in filepaths:
-            self.assertTrue(os.path.isfile(filepath))
+            self.assertTrue(pathlib.Path(filepath).is_file())
 
     def gen_reference_data(self, filepath_energy=OPEN_BOUNDARIES_REF_ENERGY,
                            filepath_arrays=OPEN_BOUNDARIES_REF_ARRAYS):
diff --git a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
index bdcca976a83..c70460dce7d 100644
--- a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
+++ b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
@@ -62,14 +62,14 @@ def test_mdlc(self):
         particle_radius = 0.5
         box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius
         s.box_l = 3 * [box_l]
-        ref_E_path = tests_common.abspath(
-            "data/mdlc_reference_data_energy.dat")
+        ref_E_path = tests_common.data_path(
+            "mdlc_reference_data_energy.dat")
         ref_E = float(np.genfromtxt(ref_E_path)) * DIPOLAR_PREFACTOR
         gap_size = 2.0
 
         # Particles
         data = np.genfromtxt(
-            tests_common.abspath("data/mdlc_reference_data_forces_torques.dat"))
+            tests_common.data_path("mdlc_reference_data_forces_torques.dat"))
         partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
         partcls.rotation = 3 * [True]
 
@@ -126,7 +126,7 @@ def test_p3m(self):
 
         # Particles
         data = np.genfromtxt(
-            tests_common.abspath("data/p3m_magnetostatics_system.data"))
+            tests_common.data_path("p3m_magnetostatics_system.data"))
         partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
         partcls.rotation = 3 * [True]
 
@@ -135,7 +135,7 @@ def test_p3m(self):
         s.actors.add(p3m)
         s.integrator.run(0)
         expected = np.genfromtxt(
-            tests_common.abspath("data/p3m_magnetostatics_expected.data"))[:, 1:]
+            tests_common.data_path("p3m_magnetostatics_expected.data"))[:, 1:]
         err_f = self.vector_error(
             partcls.f, expected[:, 0:3] * DIPOLAR_PREFACTOR)
         err_t = self.vector_error(
@@ -166,7 +166,7 @@ def test_scafacos_dipoles(self):
 
         # Particles
         data = np.genfromtxt(
-            tests_common.abspath("data/p3m_magnetostatics_system.data"))
+            tests_common.data_path("p3m_magnetostatics_system.data"))
         partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
         partcls.rotation = 3 * [True]
 
@@ -186,7 +186,7 @@ def test_scafacos_dipoles(self):
         s.actors.add(scafacos)
         s.integrator.run(0)
         expected = np.genfromtxt(
-            tests_common.abspath("data/p3m_magnetostatics_expected.data"))[:, 1:]
+            tests_common.data_path("p3m_magnetostatics_expected.data"))[:, 1:]
         err_f = self.vector_error(
             partcls.f, expected[:, 0:3] * DIPOLAR_PREFACTOR)
         err_t = self.vector_error(
diff --git a/testsuite/python/ek_eof_one_species.py b/testsuite/python/ek_eof_one_species.py
index 420b5d70ea6..b347339ec70 100644
--- a/testsuite/python/ek_eof_one_species.py
+++ b/testsuite/python/ek_eof_one_species.py
@@ -19,18 +19,16 @@
 import unittest_decorators as utx
 import unittest_generator as utg
 import pathlib
+import tempfile
+import contextlib
 
 import sys
 import math
 import numpy as np
-try:
-    import vtk
-    from vtk.util import numpy_support as VN
-    skipIfMissingPythonPackage = utx.no_skip
-except ImportError:
-    skipIfMissingPythonPackage = ut.skip(
-        "Python module vtk not available, skipping test!")
 
+with contextlib.suppress(ImportError):
+    import vtk
+    import vtk.util.numpy_support
 
 import espressomd
 import espressomd.electrokinetics
@@ -213,13 +211,14 @@ def bisection():
 
 @utx.skipIfMissingGPU()
 @utx.skipIfMissingFeatures(["ELECTROKINETICS", "EK_BOUNDARIES"])
+@utx.skipIfMissingModules("vtk")
 class ek_eof_one_species(ut.TestCase):
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
     xi = bisection()
 
     def parse_vtk(self, filepath, name, shape):
         reader = vtk.vtkStructuredPointsReader()
-        reader.SetFileName(filepath)
+        reader.SetFileName(str(filepath))
         reader.ReadAllVectorsOn()
         reader.ReadAllScalarsOn()
         reader.Update()
@@ -227,7 +226,8 @@ def parse_vtk(self, filepath, name, shape):
         data = reader.GetOutput()
         points = data.GetPointData()
 
-        return VN.vtk_to_numpy(points.GetArray(name)).reshape(shape, order='F')
+        return vtk.util.numpy_support.vtk_to_numpy(
+            points.GetArray(name)).reshape(shape, order='F')
 
     @classmethod
     def setUpClass(cls):
@@ -425,7 +425,7 @@ def test(self):
         self.assertLess(total_pressure_difference_xz, 1.0e-04,
                         "Pressure accuracy xz component not achieved")
 
-    @skipIfMissingPythonPackage
+    @utx.skipIfMissingModules("vtk")
     def test_vtk(self):
         ek = self.ek
         counterions = self.counterions
@@ -433,41 +433,44 @@ def test_vtk(self):
             map(int, np.round(self.system.box_l / params_base['agrid'])))
 
         # write VTK files
-        vtk_root = f"vtk_out/ek_eof_{axis}"
-        pathlib.Path(vtk_root).mkdir(parents=True, exist_ok=True)
-        path_vtk_boundary = f"{vtk_root}/boundary.vtk"
-        path_vtk_velocity = f"{vtk_root}/velocity.vtk"
-        path_vtk_potential = f"{vtk_root}/potential.vtk"
-        path_vtk_lbdensity = f"{vtk_root}/density.vtk"
-        path_vtk_lbforce = f"{vtk_root}/lbforce.vtk"
-        path_vtk_density = f"{vtk_root}/lbdensity.vtk"
-        path_vtk_flux = f"{vtk_root}/flux.vtk"
-        path_vtk_flux_link = f"{vtk_root}/flux_link.vtk"
-        if espressomd.has_features('EK_DEBUG'):
-            path_vtk_flux_fluc = f"{vtk_root}/flux_fluc.vtk"
-        ek.write_vtk_boundary(path_vtk_boundary)
-        ek.write_vtk_velocity(path_vtk_velocity)
-        ek.write_vtk_potential(path_vtk_potential)
-        ek.write_vtk_density(path_vtk_lbdensity)
-        ek.write_vtk_lbforce(path_vtk_lbforce)
-        counterions.write_vtk_density(path_vtk_density)
-        counterions.write_vtk_flux(path_vtk_flux)
-        if espressomd.has_features('EK_DEBUG'):
-            counterions.write_vtk_flux_fluc(path_vtk_flux_fluc)
-        counterions.write_vtk_flux_link(path_vtk_flux_link)
-
-        # load VTK files to check they are correctly formatted
-        get_vtk = self.parse_vtk
-        vtk_boundary = get_vtk(path_vtk_boundary, "boundary", grid_dims)
-        vtk_velocity = get_vtk(path_vtk_velocity, "velocity", grid_dims + [3])
-        vtk_potential = get_vtk(path_vtk_potential, "potential", grid_dims)
-        vtk_lbdensity = get_vtk(path_vtk_lbdensity, "density_lb", grid_dims)
-        get_vtk(path_vtk_lbforce, "lbforce", grid_dims + [3])
-        vtk_density = get_vtk(path_vtk_density, "density_1", grid_dims)
-        vtk_flux = get_vtk(path_vtk_flux, "flux_1", grid_dims + [3])
-        if espressomd.has_features('EK_DEBUG'):
-            get_vtk(path_vtk_flux_fluc, "flux_fluc_1", grid_dims + [4])
-        get_vtk(path_vtk_flux_link, "flux_link_1", grid_dims + [13])
+        with tempfile.TemporaryDirectory() as tmp_directory:
+            path_vtk_root = pathlib.Path(tmp_directory)
+            path_vtk_root.mkdir(parents=True, exist_ok=True)
+            path_vtk_boundary = path_vtk_root / "boundary.vtk"
+            path_vtk_velocity = path_vtk_root / "velocity.vtk"
+            path_vtk_potential = path_vtk_root / "potential.vtk"
+            path_vtk_lbdensity = path_vtk_root / "density.vtk"
+            path_vtk_lbforce = path_vtk_root / "lbforce.vtk"
+            path_vtk_density = path_vtk_root / "lbdensity.vtk"
+            path_vtk_flux = path_vtk_root / "flux.vtk"
+            path_vtk_flux_link = path_vtk_root / "flux_link.vtk"
+            if espressomd.has_features('EK_DEBUG'):
+                path_vtk_flux_fluc = path_vtk_root / "flux_fluc.vtk"
+            ek.write_vtk_boundary(str(path_vtk_boundary))
+            ek.write_vtk_velocity(str(path_vtk_velocity))
+            ek.write_vtk_potential(str(path_vtk_potential))
+            ek.write_vtk_density(str(path_vtk_lbdensity))
+            ek.write_vtk_lbforce(str(path_vtk_lbforce))
+            counterions.write_vtk_density(str(path_vtk_density))
+            counterions.write_vtk_flux(str(path_vtk_flux))
+            if espressomd.has_features('EK_DEBUG'):
+                counterions.write_vtk_flux_fluc(str(path_vtk_flux_fluc))
+            counterions.write_vtk_flux_link(str(path_vtk_flux_link))
+
+            # load VTK files to check they are correctly formatted
+            get_vtk = self.parse_vtk
+            vtk_boundary = get_vtk(path_vtk_boundary, "boundary", grid_dims)
+            vtk_velocity = get_vtk(
+                path_vtk_velocity, "velocity", grid_dims + [3])
+            vtk_potential = get_vtk(path_vtk_potential, "potential", grid_dims)
+            vtk_lbdensity = get_vtk(
+                path_vtk_lbdensity, "density_lb", grid_dims)
+            get_vtk(path_vtk_lbforce, "lbforce", grid_dims + [3])
+            vtk_density = get_vtk(path_vtk_density, "density_1", grid_dims)
+            vtk_flux = get_vtk(path_vtk_flux, "flux_1", grid_dims + [3])
+            if espressomd.has_features('EK_DEBUG'):
+                get_vtk(path_vtk_flux_fluc, "flux_fluc_1", grid_dims + [4])
+            get_vtk(path_vtk_flux_link, "flux_link_1", grid_dims + [13])
 
         # check VTK files against the EK grid
         species_density = np.zeros(grid_dims)
diff --git a/testsuite/python/h5md.py b/testsuite/python/h5md.py
index 47dbdf51023..a4919a22396 100644
--- a/testsuite/python/h5md.py
+++ b/testsuite/python/h5md.py
@@ -20,7 +20,7 @@
 """
 Testmodule for the H5MD interface.
 """
-import os
+import pathlib
 import sys
 import unittest as ut
 import unittest_decorators as utx
@@ -31,19 +31,17 @@
 import espressomd.lees_edwards
 import espressomd.version
 import tempfile
-try:
+import contextlib
+
+with contextlib.suppress(ImportError):
     import h5py  # h5py has to be imported *after* espressomd (MPI)
-    skipIfMissingPythonPackage = utx.no_skip
-except ImportError:
-    skipIfMissingPythonPackage = ut.skip(
-        "Python module h5py not available, skipping test!")
 
 
 N_PART = 26
 
 
 @utx.skipIfMissingFeatures(['H5MD'])
-@skipIfMissingPythonPackage
+@utx.skipIfMissingModules("h5py")
 class H5mdTests(ut.TestCase):
     """
     Test the core implementation of writing hdf5 files.
@@ -83,11 +81,12 @@ class H5mdTests(ut.TestCase):
     @classmethod
     def setUpClass(cls):
         cls.temp_dir = tempfile.TemporaryDirectory()
-        cls.temp_file = os.path.join(cls.temp_dir.name, 'test.h5')
+        cls.temp_path = pathlib.Path(cls.temp_dir.name)
+        cls.temp_file = cls.temp_path / 'test.h5'
         h5_units = espressomd.io.writer.h5md.UnitSystem(
             time='ps', mass='u', length='m', charge='e')
         h5 = espressomd.io.writer.h5md.H5md(
-            file_path=cls.temp_file, unit_system=h5_units)
+            file_path=str(cls.temp_file), unit_system=h5_units)
         h5.write()
         h5.write()
         h5.flush()
@@ -116,18 +115,18 @@ def tearDownClass(cls):
         cls.temp_dir.cleanup()
 
     def test_opening(self):
-        h5 = espressomd.io.writer.h5md.H5md(file_path=self.temp_file)
+        h5 = espressomd.io.writer.h5md.H5md(file_path=str(self.temp_file))
         h5.close()
 
     def test_appending(self):
         # write one frame to the file
-        temp_file = os.path.join(self.temp_dir.name, 'appending.h5')
-        h5 = espressomd.io.writer.h5md.H5md(file_path=temp_file)
+        temp_file = self.temp_path / 'appending.h5'
+        h5 = espressomd.io.writer.h5md.H5md(file_path=str(temp_file))
         h5.write()
         h5.flush()
         h5.close()
         # append one frame to the file
-        h5 = espressomd.io.writer.h5md.H5md(file_path=temp_file)
+        h5 = espressomd.io.writer.h5md.H5md(file_path=str(temp_file))
         h5.write()
         h5.flush()
         h5.close()
@@ -147,46 +146,45 @@ def predicate(cur, key):
     def test_exceptions(self):
         h5md = espressomd.io.writer.h5md
         h5_units = h5md.UnitSystem(time='ps', mass='u', length='m', charge='e')
-        temp_file = os.path.join(self.temp_dir.name, 'exceptions.h5')
+        temp_file = self.temp_path / 'exceptions.h5'
         # write a non-compliant file
-        with open(temp_file, 'wb'):
-            pass
+        temp_file.write_bytes(b'')
         with self.assertRaisesRegex(RuntimeError, 'not a valid HDF5 file'):
-            h5md.H5md(file_path=temp_file, unit_system=h5_units)
+            h5md.H5md(file_path=str(temp_file), unit_system=h5_units)
         # cannot append to a closed file with a leftover backup file
-        main_file = os.path.join(self.temp_dir.name, 'main.h5')
-        h5 = h5md.H5md(file_path=main_file)
+        main_file = self.temp_path / 'main.h5'
+        h5 = h5md.H5md(file_path=str(main_file))
         h5.write()
         h5.flush()
         h5.close()
-        with open(main_file + '.bak', 'wb') as backup:
-            with open(main_file, 'rb') as original:
-                backup.write(original.read())
+        main_file.with_suffix(temp_file.suffix + '.bak').write_bytes(b'')
         with self.assertRaisesRegex(RuntimeError, 'A backup of the .h5 file exists'):
-            h5md.H5md(file_path=main_file)
-        # cannot create a new file with a leftover backup file
-        os.remove(main_file)
+            h5md.H5md(file_path=str(main_file))
+        # cannot create a new file when a leftover backup file exists
+        main_file.unlink()
         with self.assertRaisesRegex(RuntimeError, 'A backup of the .h5 file exists'):
-            h5md.H5md(file_path=main_file)
+            h5md.H5md(file_path=str(main_file))
         # open a file with different specifications
-        temp_file = os.path.join(self.temp_dir.name, 'wrong_spec.h5')
-        h5 = espressomd.io.writer.h5md.H5md(file_path=temp_file, fields=[])
+        temp_file = self.temp_path / 'wrong_spec.h5'
+        h5 = espressomd.io.writer.h5md.H5md(
+            file_path=str(temp_file), fields=[])
         h5.write()
         h5.flush()
         h5.close()
         with self.assertRaisesRegex(RuntimeError, "The given .h5 file does not match the specifications in 'fields'"):
-            h5md.H5md(file_path=temp_file, fields='all')
+            h5md.H5md(file_path=str(temp_file), fields='all')
         # open a file with invalid specifications
         with self.assertRaisesRegex(ValueError, "Unknown field 'lb'"):
-            h5md.H5md(file_path=temp_file, fields='lb')
+            h5md.H5md(file_path=str(temp_file), fields='lb')
         # check read-only parameters
         for key in self.h5_obj.get_params():
             with self.assertRaisesRegex(RuntimeError, f"Parameter '{key}' is read-only"):
                 setattr(self.h5_obj, key, None)
 
     def test_empty(self):
-        temp_file = os.path.join(self.temp_dir.name, 'empty.h5')
-        h5 = espressomd.io.writer.h5md.H5md(file_path=temp_file, fields=[])
+        temp_file = self.temp_path / 'empty.h5'
+        h5 = espressomd.io.writer.h5md.H5md(
+            file_path=str(temp_file), fields=[])
         h5.write()
         h5.flush()
         h5.close()
@@ -249,7 +247,7 @@ def test_img(self):
 
     @utx.skipIfMissingFeatures("MASS")
     def test_mass(self):
-        """Test if masses have been written correct."""
+        """Test if masses have been written properly."""
         np.testing.assert_allclose(self.py_mass, 2.3)
 
     @utx.skipIfMissingFeatures(['ELECTROSTATICS'])
@@ -291,9 +289,9 @@ def test_script(self):
         data = self.py_file['parameters/files'].attrs['script'].decode('utf-8')
         self.assertEqual(data, ref)
         # case #2: running an interactive Python session
-        temp_file = os.path.join(self.temp_dir.name, 'no_script.h5')
+        temp_file = self.temp_path / 'no_script.h5'
         sys.argv[0] = ''
-        h5 = espressomd.io.writer.h5md.H5md(file_path=temp_file)
+        h5 = espressomd.io.writer.h5md.H5md(file_path=str(temp_file))
         sys.argv[0] = __file__
         h5.write()
         h5.flush()
@@ -317,9 +315,9 @@ def get_unit(path):
         self.assertEqual(get_unit('particles/atoms/velocity/value'), b'm ps-1')
 
     def test_getters(self):
-        self.assertEqual(self.h5_params['file_path'], self.temp_file)
-        self.assertEqual(os.path.abspath(self.h5_params['script_path']),
-                         os.path.abspath(__file__))
+        self.assertEqual(self.h5_params['file_path'], str(self.temp_file))
+        self.assertEqual(pathlib.Path(self.h5_params['script_path']).resolve(),
+                         pathlib.Path(__file__).resolve())
         self.assertEqual(self.h5_params['fields'], ['all'])
         self.assertEqual(self.h5_params['time_unit'], 'ps')
         if espressomd.has_features(['ELECTROSTATICS']):
diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py
index 95a8828022f..3e9eaa44427 100644
--- a/testsuite/python/integrator_npt.py
+++ b/testsuite/python/integrator_npt.py
@@ -154,10 +154,8 @@ def test_dp3m_exception(self):
             self.run_with_p3m(
                 dp3m, cubic_box=False, direction=(False, True, True))
         self.tearDown()
-        try:
-            self.run_with_p3m(dp3m)
-        except Exception as err:
-            self.fail(f'integrator raised ValueError("{err}")')
+        # should not raise an exception
+        self.run_with_p3m(dp3m)
 
     @utx.skipIfMissingFeatures(["P3M"])
     def test_p3m_exception(self):
@@ -170,10 +168,8 @@ def test_p3m_exception(self):
             self.run_with_p3m(
                 p3m, cubic_box=False, direction=(False, True, True))
         self.tearDown()
-        try:
-            self.run_with_p3m(p3m)
-        except Exception as err:
-            self.fail(f'integrator raised ValueError("{err}")')
+        # should not raise an exception
+        self.run_with_p3m(p3m)
 
     @utx.skipIfMissingGPU()
     @utx.skipIfMissingFeatures(["P3M"])
diff --git a/testsuite/python/integrator_npt_stats.py b/testsuite/python/integrator_npt_stats.py
index f9d15d428e2..54e2cde202f 100644
--- a/testsuite/python/integrator_npt_stats.py
+++ b/testsuite/python/integrator_npt_stats.py
@@ -44,7 +44,7 @@ def test_compressibility(self):
         system = self.system
         system.box_l = [5.86326165] * 3
 
-        data = np.genfromtxt(tests_common.abspath("data/npt_lj_system.data"))
+        data = np.genfromtxt(tests_common.data_path("npt_lj_system.data"))
         p_ext = 2.0
 
         system.part.add(pos=data[:, :3], v=data[:, 3:])
diff --git a/testsuite/python/interactions_dihedral.py b/testsuite/python/interactions_dihedral.py
index 91c1cbb5c45..0c68450c0e2 100644
--- a/testsuite/python/interactions_dihedral.py
+++ b/testsuite/python/interactions_dihedral.py
@@ -22,28 +22,19 @@
 
 import espressomd
 
-# Dihedral interaction needs more rigorous tests.
-# The geometry checked here is rather simple and special.
-# I also found that as the dihedral angle approaches to 0, the simulation
-# values deviate from the analytic values by roughly 10%.
-
 
 def rotate_vector(v, k, phi):
     """Rotates vector v around unit vector k by angle phi.
     Uses Rodrigues' rotation formula."""
-    vrot = v * np.cos(phi) + np.cross(k, v) * \
-        np.sin(phi) + k * np.dot(k, v) * (1.0 - np.cos(phi))
+    vrot = np.array(v) * np.cos(phi) + np.cross(k, v) * \
+        np.sin(phi) + np.array(k) * np.dot(k, v) * (1.0 - np.cos(phi))
     return vrot
 
 
-def dihedral_potential(k, phi, n, phase):
-    if phi == -1:
-        return 0
-    else:
-        return k * (1 - np.cos(n * phi - phase))
-
-
-def dihedral_force(k, n, phase, p1, p2, p3, p4):
+def dihedral_potential_and_forces(k, n, phase, p1, p2, p3, p4):
+    """
+    Calculate the potential and forces for a dihedral angle.
+    """
     v12 = p2 - p1
     v23 = p3 - p2
     v34 = p4 - p3
@@ -52,158 +43,126 @@ def dihedral_force(k, n, phase, p1, p2, p3, p4):
     l_v12Xv23 = np.linalg.norm(v12Xv23)
     v23Xv34 = np.cross(v23, v34)
     l_v23Xv34 = np.linalg.norm(v23Xv34)
-    # if dihedral angle is not defined, no forces
-    if l_v12Xv23 <= 1e-8 or l_v23Xv34 <= 1e-8:
-        return 0, 0, 0
-    else:
-        cosphi = np.abs(np.dot(v12Xv23, v23Xv34)) / (l_v12Xv23 * l_v23Xv34)
-        phi = np.arccos(cosphi)
-        f1 = (v23Xv34 - cosphi * v12Xv23) / l_v12Xv23
-        f4 = (v12Xv23 - cosphi * v23Xv34) / l_v23Xv34
 
-        v23Xf1 = np.cross(v23, f1)
-        v23Xf4 = np.cross(v23, f4)
-        v34Xf4 = np.cross(v34, f4)
-        v12Xf1 = np.cross(v12, f1)
+    phi = np.arctan2(np.dot(v23, np.cross(v12Xv23, v23Xv34)),
+                     np.dot(v23, v23) * np.dot(v12Xv23, v23Xv34))
+
+    f1 = (v23Xv34 - np.cos(phi) * v12Xv23) / l_v12Xv23
+    f4 = (v12Xv23 - np.cos(phi) * v23Xv34) / l_v23Xv34
 
+    v23Xf1 = np.cross(v23, f1)
+    v23Xf4 = np.cross(v23, f4)
+    v34Xf4 = np.cross(v34, f4)
+    v12Xf1 = np.cross(v12, f1)
+
+    # handle singularity near TINY_SIN_VALUE
+    if np.abs(np.sin(phi)) < 1e-10:
+        coeff = -k * n**2 * np.cos(n * phi - phase) / np.cos(phi)
+    else:
         coeff = -k * n * np.sin(n * phi - phase) / np.sin(phi)
 
-        force1 = coeff * v23Xf1
-        force2 = coeff * (v34Xf4 - v12Xf1 - v23Xf1)
-        force3 = coeff * (v12Xf1 - v23Xf4 - v34Xf4)
-        return force1, force2, force3
+    force1 = coeff * v23Xf1
+    force2 = coeff * (v34Xf4 - v12Xf1 - v23Xf1)
+    force3 = coeff * (v12Xf1 - v23Xf4 - v34Xf4)
+    force4 = coeff * v23Xf4
+    potential = k * (1 - np.cos(n * phi - phase))
+    return (potential, (force1, force2, force3, force4))
 
 
 class InteractionsBondedTest(ut.TestCase):
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system = espressomd.System(box_l=[10.0, 10.0, 10.0])
+    system.cell_system.skin = 0.4
+    system.time_step = 0.1
     np.random.seed(seed=42)
 
-    box_l = 10.
-
-    start_pos = [5., 5., 5.]
-    axis = np.array([1., 0., 0.])
-    axis /= np.linalg.norm(axis)
-    rel_pos_1 = np.array([0., 1., 0.])
-    rel_pos_2 = np.array([0., 0., 1.])
-
-    def setUp(self):
-
-        self.system.box_l = [self.box_l] * 3
-        self.system.cell_system.skin = 0.4
-        self.system.time_step = .1
-
-        self.system.part.add(pos=4 * [self.start_pos], type=4 * [0])
-
     def tearDown(self):
         self.system.part.clear()
 
-    # Analytical Expression
-    def dihedral_angle(self, p1, p2, p3, p4):
-        """
-        Calculate the dihedral angle phi based on particles' position p1, p2, p3, p4.
-        """
-        v12 = p2 - p1
-        v23 = p3 - p2
-        v34 = p4 - p3
-
-        v12Xv23 = np.cross(v12, v23)
-        l_v12Xv23 = np.linalg.norm(v12Xv23)
-        v23Xv34 = np.cross(v23, v34)
-        l_v23Xv34 = np.linalg.norm(v23Xv34)
-
-        # if dihedral angle is not defined, phi := -1.
-        if l_v12Xv23 <= 1e-8 or l_v23Xv34 <= 1e-8:
-            return -1
-        else:
-            cosphi = np.abs(np.dot(v12Xv23, v23Xv34)) / (
-                l_v12Xv23 * l_v23Xv34)
-            return np.arccos(cosphi)
+    def check_values(self, E_ref, forces_ref, tol=1e-12):
+        E_sim = self.system.analysis.energy()["bonded"]
+        np.testing.assert_allclose(E_sim, E_ref, atol=tol)
+        if forces_ref:
+            f0, f1, f2, f3 = self.system.part.all().f
+            f0_ref, f1_ref, f2_ref, f3_ref = forces_ref
+            np.testing.assert_allclose(np.copy(f0), f0_ref, atol=tol)
+            np.testing.assert_allclose(np.copy(f1), f1_ref, atol=tol)
+            np.testing.assert_allclose(np.copy(f2), f2_ref, atol=tol)
+            np.testing.assert_allclose(np.copy(f3), f3_ref, atol=tol)
 
     # Test Dihedral Angle
     def test_dihedral(self):
-        p0, p1, p2, p3 = self.system.part.all()
-
-        dh_k = 1
-        dh_phase = np.pi / 6
-        dh_n = 1
-
-        dh = espressomd.interactions.Dihedral(
-            bend=dh_k, mult=dh_n, phase=dh_phase)
-        self.system.bonded_inter.add(dh)
-        p1.add_bond((dh, p0, p2, p3))
-        p2.pos = p1.pos + [1, 0, 0]
-
-        N = 111
-        d_phi = np.pi / (N * 4)
-        for i in range(N):
-            p0.pos = p1.pos + \
-                rotate_vector(self.rel_pos_1, self.axis, i * d_phi)
-            p3.pos = p2.pos + \
-                rotate_vector(self.rel_pos_2, self.axis, -i * d_phi)
-            self.system.integrator.run(recalc_forces=True, steps=0)
-
-            # Calculate energies
-            E_sim = self.system.analysis.energy()["bonded"]
-            phi = self.dihedral_angle(p0.pos, p1.pos, p2.pos, p3.pos)
-            E_ref = dihedral_potential(dh_k, phi, dh_n, dh_phase)
-
-            # Calculate forces
-            f2_sim = p1.f
-            _, f2_ref, _ = dihedral_force(dh_k, dh_n, dh_phase,
-                                          p0.pos, p1.pos, p2.pos, p3.pos)
-
-            # Check that energies match, ...
-            np.testing.assert_almost_equal(E_sim, E_ref)
-            # and has correct value.
-            f2_sim_copy = np.copy(f2_sim)
-            np.testing.assert_almost_equal(f2_sim_copy, f2_ref)
+        axis = np.array([1., 0., 0.])
+        dh_k = 2.
+        N = 100  # even number to get singularities at phi=0 and phi=pi
+        d_phi = 2 * np.pi / N
+        for dh_n, dh_phi0_div in [(2, 3), (3, 6)]:
+            with self.subTest(multiplicity=dh_n, phi_0=f"pi / {dh_phi0_div}"):
+                dh_phi0 = np.pi / dh_phi0_div
+                dihedral = espressomd.interactions.Dihedral(
+                    bend=dh_k, mult=dh_n, phase=dh_phi0)
+                self.system.bonded_inter.add(dihedral)
+                self.system.part.clear()
+                p0, p1, p2, p3 = self.system.part.add(pos=4 * [(0., 0., 0.)])
+                p1.add_bond((dihedral, p0, p2, p3))
+                p1.pos = [5., 5., 5.]
+                p2.pos = p1.pos + [1., 0., 0.]
+                p0.pos = p1.pos + [0., 1., 0.]
+
+                for i in range(N):
+                    phi = i * d_phi
+                    p3.pos = p2.pos + rotate_vector([0., 1., 0.], axis, phi)
+                    self.system.integrator.run(recalc_forces=True, steps=0)
+
+                    # Calculate expected forces and energies
+                    E_ref, forces_ref = dihedral_potential_and_forces(
+                        dh_k, dh_n, dh_phi0, p0.pos, p1.pos, p2.pos, p3.pos)
+
+                    self.check_values(E_ref, forces_ref)
 
     # Test Tabulated Dihedral Angle
     @utx.skipIfMissingFeatures(["TABULATED"])
     def test_tabulated_dihedral(self):
-        p0, p1, p2, p3 = self.system.part.all()
-
-        N = 111
+        axis = np.array([1., 0., 0.])
+        dh_k = 2.
+        N = 100  # even number to get singularities at phi=0 and phi=pi
         d_phi = 2 * np.pi / N
-        # tabulated values for the range [0, 2*pi]
-        tab_energy = [np.cos(i * d_phi) for i in range(N + 1)]
-        tab_force = [np.cos(i * d_phi) for i in range(N + 1)]
-
-        dihedral_tabulated = espressomd.interactions.TabulatedDihedral(
-            energy=tab_energy, force=tab_force)
-        self.system.bonded_inter.add(dihedral_tabulated)
-        p1.add_bond((dihedral_tabulated, p0, p2, p3))
-        p2.pos = p1.pos + [1, 0, 0]
-
-        # check stored parameters
-        interaction_id = len(self.system.bonded_inter) - 1
-        tabulated = self.system.bonded_inter[interaction_id]
-        np.testing.assert_allclose(tabulated.params['force'], tab_force)
-        np.testing.assert_allclose(tabulated.params['energy'], tab_energy)
-        np.testing.assert_almost_equal(tabulated.params['min'], 0.)
-        np.testing.assert_almost_equal(tabulated.params['max'], 2 * np.pi)
-
-        # measure at half the angular resolution to observe interpolation
-        for i in range(2 * N - 1):
-            # increase dihedral angle by d_phi (phi ~ 0 at i = 0)
-            p0.pos = p1.pos + \
-                rotate_vector(self.rel_pos_1, self.axis, -i * d_phi / 4)
-            p3.pos = p2.pos + \
-                rotate_vector(self.rel_pos_1, self.axis, i * d_phi / 4)
-            self.system.integrator.run(recalc_forces=True, steps=0)
-
-            # Calculate energies
-            E_sim = self.system.analysis.energy()["bonded"]
-
-            # Get tabulated values
-            j = i // 2
-            if i % 2 == 0:
-                E_ref = tab_energy[j]
-            else:
-                E_ref = (tab_energy[j] + tab_energy[j + 1]) / 2.0
-
-            # Check that energies match, ...
-            np.testing.assert_almost_equal(E_sim, E_ref)
+        for dh_n, dh_phi0_div in [(2, 3), (3, 6)]:
+            with self.subTest(multiplicity=dh_n, phi_0=f"pi / {dh_phi0_div}"):
+                dh_phi0 = np.pi / dh_phi0_div
+                # tabulated values for the range [0, 2*pi]
+                phi = d_phi * np.arange(N + 1)
+                tab_energy = dh_k * (1. - np.cos(dh_n * phi - dh_phi0))
+                div = np.sin(phi)
+                div[0] = div[N // 2] = div[N] = 1.
+                tab_force = -dh_k * dh_n * np.sin(dh_n * phi - dh_phi0) / div
+                tab_force[0] = tab_force[N // 2] = tab_force[N] = 0.
+                dihedral_tabulated = espressomd.interactions.TabulatedDihedral(
+                    energy=tab_energy, force=tab_force)
+                self.system.bonded_inter.add(dihedral_tabulated)
+                self.system.part.clear()
+                p0, p1, p2, p3 = self.system.part.add(pos=4 * [(0., 0., 0.)])
+                p1.add_bond((dihedral_tabulated, p0, p2, p3))
+                p1.pos = [5., 5., 5.]
+                p2.pos = p1.pos + [1., 0., 0.]
+                p0.pos = p1.pos + [0., 1., 0.]
+
+                # use half the angular resolution to observe interpolation
+                for i in range(2 * N - 1):
+                    phi = i * d_phi / 2.
+                    p3.pos = p2.pos + rotate_vector([0., 1., 0.], axis, phi)
+                    self.system.integrator.run(recalc_forces=True, steps=0)
+
+                    # Calculate expected forces and energies
+                    j = i // 2
+                    if i % 2 == 0:
+                        E_ref = tab_energy[j]
+                        _, forces_ref = dihedral_potential_and_forces(
+                            dh_k, dh_n, dh_phi0, p0.pos, p1.pos, p2.pos, p3.pos)
+                    else:
+                        E_ref = (tab_energy[j] + tab_energy[j + 1]) / 2.0
+                        forces_ref = None
+
+                    self.check_values(E_ref, forces_ref)
 
 
 if __name__ == '__main__':
diff --git a/testsuite/python/lb_vtk.py b/testsuite/python/lb_vtk.py
index 2e4e28f8088..d48d3a1148e 100644
--- a/testsuite/python/lb_vtk.py
+++ b/testsuite/python/lb_vtk.py
@@ -19,16 +19,14 @@
 import unittest as ut
 import unittest_decorators as utx
 
-import os
+import pathlib
+import tempfile
+import contextlib
 import numpy as np
 
-try:
+with contextlib.suppress(ImportError):
     import vtk
-    from vtk.util import numpy_support as VN
-    skipIfMissingPythonPackage = utx.no_skip
-except ImportError:
-    skipIfMissingPythonPackage = ut.skip(
-        "Python module vtk not available, skipping test!")
+    import vtk.util.numpy_support
 
 import espressomd
 import espressomd.lb
@@ -61,7 +59,7 @@ def set_lbf(self):
 
     def parse_vtk(self, filepath, name, shape):
         reader = vtk.vtkStructuredPointsReader()
-        reader.SetFileName(filepath)
+        reader.SetFileName(str(filepath))
         reader.ReadAllVectorsOn()
         reader.ReadAllScalarsOn()
         reader.Update()
@@ -69,146 +67,130 @@ def parse_vtk(self, filepath, name, shape):
         data = reader.GetOutput()
         points = data.GetPointData()
 
-        return VN.vtk_to_numpy(points.GetArray(name)).reshape(shape, order='F')
+        return vtk.util.numpy_support.vtk_to_numpy(
+            points.GetArray(name)).reshape(shape, order='F')
 
     def test_vtk(self):
         '''
         Check VTK files.
         '''
 
-        os.makedirs('vtk_out', exist_ok=True)
-        filepaths = ['vtk_out/boundary.vtk', 'vtk_out/velocity.vtk',
-                     'vtk_out/velocity_bb.vtk']
-
-        # cleanup action
-        for filepath in filepaths:
-            if os.path.exists(filepath):
-                os.remove(filepath)
-
-        shape = [10, 11, 12]
-        lbf = self.set_lbf()
-        self.system.integrator.run(100)
-
-        # write VTK files
-        with self.assertRaises(RuntimeError):
-            lbf.write_vtk_velocity('non_existent_folder/file')
-        with self.assertRaises(RuntimeError):
-            lbf.write_vtk_boundary('non_existent_folder/file')
-        lbf.write_vtk_boundary('vtk_out/boundary.vtk')
-        lbf.write_vtk_velocity('vtk_out/velocity.vtk')
-        with self.assertRaises(ValueError):
-            lbf.write_vtk_velocity('vtk_out/delme', 3 * [0], None)
-        with self.assertRaises(ValueError):
-            lbf.write_vtk_velocity('vtk_out/delme', None, 3 * [0])
-        with self.assertRaises(RuntimeError):
-            lbf.write_vtk_velocity('vtk_out/delme', [-2, 1, 1], 3 * [1])
-        with self.assertRaises(RuntimeError):
-            lbf.write_vtk_velocity('vtk_out/delme', 3 * [0], [1, 2, 16])
-        with self.assertRaises(ValueError):
-            lbf.write_vtk_velocity('vtk_out/delme', [1, 1], 3 * [1])
-        with self.assertRaises(ValueError):
-            lbf.write_vtk_velocity('vtk_out/delme', 3 * [1], np.array([2, 3]))
-        bb1, bb2 = ([1, 2, 3], [9, 10, 11])
-        lbf.write_vtk_velocity('vtk_out/velocity_bb.vtk', bb1, bb2)
-
-        # check VTK files exist
-        for filepath in filepaths:
-            self.assertTrue(
-                os.path.exists(filepath),
-                f'VTK file "{filepath}" not written to disk')
-
-        # check VTK values match node values
-        node_velocity = np.zeros(shape + [3])
-        node_boundary = np.zeros(shape, dtype=int)
-        for i in range(shape[0]):
-            for j in range(shape[1]):
-                for k in range(shape[2]):
-                    node = lbf[i, j, k]
-                    node_velocity[i, j, k] = node.velocity
-                    node_boundary[i, j, k] = node.boundary
-        node_velocity_bb = node_velocity[bb1[0]:bb2[0],
-                                         bb1[1]:bb2[1],
-                                         bb1[2]:bb2[2]]
-
-        vtk_velocity = self.parse_vtk('vtk_out/velocity.vtk', 'velocity',
-                                      node_velocity.shape)
-        np.testing.assert_allclose(vtk_velocity, node_velocity, atol=5e-7)
-
-        vtk_velocity_bb = self.parse_vtk('vtk_out/velocity_bb.vtk', 'velocity',
-                                         node_velocity_bb.shape)
-        np.testing.assert_allclose(
-            vtk_velocity_bb, node_velocity_bb, atol=5e-7)
-
-        vtk_boundary = self.parse_vtk(
-            'vtk_out/boundary.vtk', 'boundary', shape)
-        np.testing.assert_equal(vtk_boundary, node_boundary.astype(int))
-        if self.system.lbboundaries is None:
-            np.testing.assert_equal(np.sum(node_boundary), 0.)
+        with tempfile.TemporaryDirectory() as tmp_directory:
+            path_vtk_root = pathlib.Path(tmp_directory)
+            path_vtk_boundary = path_vtk_root / 'boundary.vtk'
+            path_vtk_velocity = path_vtk_root / 'velocity.vtk'
+            path_vtk_velocity_bb = path_vtk_root / 'velocity_bb.vtk'
+            path_vtk_skip = path_vtk_root / 'skip.vtk'
+            path_vtk_invalid = path_vtk_root / 'non_existent_folder' / 'file'
+
+            shape = [10, 11, 12]
+            lbf = self.set_lbf()
+            self.system.integrator.run(100)
+
+            # write VTK files
+            with self.assertRaises(RuntimeError):
+                lbf.write_vtk_velocity(str(path_vtk_invalid))
+            with self.assertRaises(RuntimeError):
+                lbf.write_vtk_boundary(str(path_vtk_invalid))
+            lbf.write_vtk_boundary(str(path_vtk_boundary))
+            lbf.write_vtk_velocity(str(path_vtk_velocity))
+            with self.assertRaises(ValueError):
+                lbf.write_vtk_velocity(str(path_vtk_skip), 3 * [0], None)
+            with self.assertRaises(ValueError):
+                lbf.write_vtk_velocity(str(path_vtk_skip), None, 3 * [0])
+            with self.assertRaises(RuntimeError):
+                lbf.write_vtk_velocity(str(path_vtk_skip), [-2, 1, 1], 3 * [1])
+            with self.assertRaises(RuntimeError):
+                lbf.write_vtk_velocity(str(path_vtk_skip), 3 * [0], [1, 2, 16])
+            with self.assertRaises(ValueError):
+                lbf.write_vtk_velocity(str(path_vtk_skip), [1, 1], 3 * [1])
+            with self.assertRaises(ValueError):
+                lbf.write_vtk_velocity(
+                    str(path_vtk_skip), 3 * [1], np.array([2, 3]))
+            bb1, bb2 = ([1, 2, 3], [9, 10, 11])
+            lbf.write_vtk_velocity(str(path_vtk_velocity_bb), bb1, bb2)
+
+            # check VTK values match node values
+            node_velocity = np.zeros(shape + [3])
+            node_boundary = np.zeros(shape, dtype=int)
+            for i in range(shape[0]):
+                for j in range(shape[1]):
+                    for k in range(shape[2]):
+                        node = lbf[i, j, k]
+                        node_velocity[i, j, k] = node.velocity
+                        node_boundary[i, j, k] = node.boundary
+            node_velocity_bb = node_velocity[bb1[0]:bb2[0],
+                                             bb1[1]:bb2[1],
+                                             bb1[2]:bb2[2]]
+
+            vtk_velocity = self.parse_vtk(path_vtk_velocity, 'velocity',
+                                          node_velocity.shape)
+            np.testing.assert_allclose(vtk_velocity, node_velocity, atol=5e-7)
+
+            vtk_velocity_bb = self.parse_vtk(path_vtk_velocity_bb, 'velocity',
+                                             node_velocity_bb.shape)
+            np.testing.assert_allclose(
+                vtk_velocity_bb, node_velocity_bb, atol=5e-7)
+
+            vtk_boundary = self.parse_vtk(path_vtk_boundary, 'boundary', shape)
+            np.testing.assert_equal(vtk_boundary, node_boundary.astype(int))
 
     def test_print(self):
         '''
         Check data files.
         '''
 
-        os.makedirs('vtk_out', exist_ok=True)
-        filepaths = ['vtk_out/boundary.dat', 'vtk_out/velocity.dat']
-
-        # cleanup action
-        for filepath in filepaths:
-            if os.path.exists(filepath):
-                os.remove(filepath)
-
-        shape = [10, 11, 12]
-        lbf = self.set_lbf()
-        self.system.integrator.run(100)
-
-        # write data files
-        with self.assertRaises(RuntimeError):
-            lbf.write_velocity('non_existent_folder/file')
-        with self.assertRaises(RuntimeError):
-            lbf.write_boundary('non_existent_folder/file')
-        lbf.write_boundary('vtk_out/boundary.dat')
-        lbf.write_velocity('vtk_out/velocity.dat')
-
-        # check data files exist
-        for filepath in filepaths:
-            self.assertTrue(
-                os.path.exists(filepath),
-                f'data file "{filepath}" not written to disk')
-
-        # check data values match node values
-        node_velocity = np.zeros(shape + [3])
-        node_boundary = np.zeros(shape, dtype=int)
-        for i in range(shape[0]):
-            for j in range(shape[1]):
-                for k in range(shape[2]):
-                    node = lbf[i, j, k]
-                    node_velocity[i, j, k] = node.velocity
-                    node_boundary[i, j, k] = node.boundary
-
-        ref_coord = np.array([
-            np.tile(np.arange(shape[0]), shape[1] * shape[2]),
-            np.tile(np.repeat(np.arange(shape[1]), shape[0]), shape[2]),
-            np.repeat(np.arange(shape[2]), shape[0] * shape[1])]).T
-
-        dat_velocity = np.loadtxt('vtk_out/velocity.dat')
-        dat_coord = (dat_velocity[:, 0:3] - 0.5).astype(int)
-        np.testing.assert_equal(dat_coord, ref_coord)
-        dat_vel = dat_velocity[:, 3:]
-        ref_vel = np.swapaxes(node_velocity, 0, 2).reshape((-1, 3))
-        np.testing.assert_allclose(dat_vel, ref_vel, atol=5e-7)
-
-        dat_boundary = np.loadtxt('vtk_out/boundary.dat')
-        dat_coord = (dat_boundary[:, 0:3] - 0.5).astype(int)
-        np.testing.assert_equal(dat_coord, ref_coord)
-        dat_bound = dat_boundary[:, 3].astype(int)
-        ref_bound = np.swapaxes(node_boundary, 0, 2).reshape(-1)
-        if isinstance(lbf, espressomd.lb.LBFluid):
-            ref_bound = (ref_bound != 0).astype(int)
-        np.testing.assert_equal(dat_bound, ref_bound)
-
-
-@skipIfMissingPythonPackage
+        with tempfile.TemporaryDirectory() as tmp_directory:
+            path_dat_root = pathlib.Path(tmp_directory)
+            path_dat_boundary = path_dat_root / 'boundary.vtk'
+            path_dat_velocity = path_dat_root / 'velocity.vtk'
+            path_dat_invalid = path_dat_root / 'non_existent_folder' / 'file'
+
+            shape = [10, 11, 12]
+            lbf = self.set_lbf()
+            self.system.integrator.run(100)
+
+            # write data files
+            with self.assertRaises(RuntimeError):
+                lbf.write_velocity(str(path_dat_invalid))
+            with self.assertRaises(RuntimeError):
+                lbf.write_boundary(str(path_dat_invalid))
+            lbf.write_boundary(str(path_dat_boundary))
+            lbf.write_velocity(str(path_dat_velocity))
+
+            # check data values match node values
+            node_velocity = np.zeros(shape + [3])
+            node_boundary = np.zeros(shape, dtype=int)
+            for i in range(shape[0]):
+                for j in range(shape[1]):
+                    for k in range(shape[2]):
+                        node = lbf[i, j, k]
+                        node_velocity[i, j, k] = node.velocity
+                        node_boundary[i, j, k] = node.boundary
+
+            ref_coord = np.array([
+                np.tile(np.arange(shape[0]), shape[1] * shape[2]),
+                np.tile(np.repeat(np.arange(shape[1]), shape[0]), shape[2]),
+                np.repeat(np.arange(shape[2]), shape[0] * shape[1])]).T
+
+            dat_velocity = np.loadtxt(path_dat_velocity)
+            dat_coord = (dat_velocity[:, 0:3] - 0.5).astype(int)
+            np.testing.assert_equal(dat_coord, ref_coord)
+            dat_vel = dat_velocity[:, 3:]
+            ref_vel = np.swapaxes(node_velocity, 0, 2).reshape((-1, 3))
+            np.testing.assert_allclose(dat_vel, ref_vel, atol=5e-7)
+
+            dat_boundary = np.loadtxt(path_dat_boundary)
+            dat_coord = (dat_boundary[:, 0:3] - 0.5).astype(int)
+            np.testing.assert_equal(dat_coord, ref_coord)
+            dat_bound = dat_boundary[:, 3].astype(int)
+            ref_bound = np.swapaxes(node_boundary, 0, 2).reshape(-1)
+            if isinstance(lbf, espressomd.lb.LBFluid):
+                ref_bound = (ref_bound != 0).astype(int)
+            np.testing.assert_equal(dat_bound, ref_bound)
+
+
+@utx.skipIfMissingModules("vtk")
 class TestLBWriteCPU(TestLBWrite, ut.TestCase):
 
     def setUp(self):
@@ -216,7 +198,7 @@ def setUp(self):
 
 
 @utx.skipIfMissingGPU()
-@skipIfMissingPythonPackage
+@utx.skipIfMissingModules("vtk")
 class TestLBWriteGPU(TestLBWrite, ut.TestCase):
 
     def setUp(self):
diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py
index f0e47174bbb..731eb09c4ac 100644
--- a/testsuite/python/lees_edwards.py
+++ b/testsuite/python/lees_edwards.py
@@ -30,8 +30,9 @@
 
 
 np.random.seed(42)
-params_lin = {'shear_velocity': 1.2, 'initial_pos_offset': 0.1, 'time_0': 0.1}
-params_osc = {'amplitude': 2.3, 'omega': 2.51, 'time_0': -2.1}
+params_lin = {'initial_pos_offset': 0.1, 'time_0': 0.1, 'shear_velocity': 1.2}
+params_osc = {'initial_pos_offset': 0.1, 'time_0': -2.1, 'amplitude': 2.3,
+              'omega': 2.51}
 lin_protocol = espressomd.lees_edwards.LinearShear(**params_lin)
 osc_protocol = espressomd.lees_edwards.OscillatoryShear(**params_osc)
 off_protocol = espressomd.lees_edwards.Off()
@@ -113,7 +114,7 @@ def test_protocols(self):
         # check that LE offsets are recalculated on simulation time change
         for time in [0., 2.3]:
             system.time = time
-            expected_pos = params_osc['amplitude'] * \
+            expected_pos = params_osc['initial_pos_offset'] + params_osc['amplitude'] * \
                 np.sin(params_osc['omega'] * (time - params_osc['time_0']))
             expected_vel = params_osc['amplitude'] * params_osc['omega'] * \
                 np.cos(params_osc['omega'] * (time - params_osc['time_0']))
@@ -124,7 +125,7 @@ def test_protocols(self):
         # Check that time change during integration updates offsets
         system.integrator.run(1)
         time = system.time
-        expected_pos = params_osc['amplitude'] * \
+        expected_pos = params_osc['initial_pos_offset'] + params_osc['amplitude'] * \
             np.sin(params_osc['omega'] * (time - params_osc['time_0']))
         expected_vel = params_osc['amplitude'] * params_osc['omega'] * \
             np.cos(params_osc['omega'] * (time - params_osc['time_0']))
diff --git a/testsuite/python/lj.py b/testsuite/python/lj.py
index aaa6ad46e52..1fdd104beda 100644
--- a/testsuite/python/lj.py
+++ b/testsuite/python/lj.py
@@ -26,7 +26,7 @@
 @utx.skipIfMissingFeatures(["LENNARD_JONES"])
 class LennardJonesTest(ut.TestCase):
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
-    data = np.loadtxt(tests_common.abspath('data/lj_system.dat'))
+    data = np.loadtxt(tests_common.data_path('lj_system.dat'))
     pos = data[:, 1:4]
     forces = data[:, 4:7]
 
diff --git a/testsuite/python/mdanalysis.py b/testsuite/python/mdanalysis.py
index 8a2659f10c4..8d07e58df78 100644
--- a/testsuite/python/mdanalysis.py
+++ b/testsuite/python/mdanalysis.py
@@ -25,16 +25,14 @@
 import numpy as np
 import unittest as ut
 import unittest_decorators as utx
-try:
-    import MDAnalysis as mda
+import contextlib
+
+with contextlib.suppress(ImportError):
+    import MDAnalysis
     import espressomd.MDA_ESP
-    skipIfMissingPythonPackage = utx.no_skip
-except ImportError:
-    skipIfMissingPythonPackage = ut.skip(
-        "Python module MDAnalysis not available, skipping test!")
 
 
-@skipIfMissingPythonPackage
+@utx.skipIfMissingModules("MDAnalysis")
 class TestMDAnalysis(ut.TestCase):
     system = espressomd.System(box_l=[10.0, 10.0, 10.0])
     system.time_step = 0.001
@@ -58,7 +56,7 @@ def test_universe(self):
         system = self.system
         partcls = system.part.all()
         eos = espressomd.MDA_ESP.Stream(system)
-        u = mda.Universe(eos.topology, eos.trajectory)
+        u = MDAnalysis.Universe(eos.topology, eos.trajectory)
         # check atoms
         self.assertEqual(len(u.atoms), 10)
         np.testing.assert_equal(u.atoms.ids, np.arange(10) + 1)
diff --git a/testsuite/python/mmm1d.py b/testsuite/python/mmm1d.py
index a63c56afa09..c0c7f6f624b 100644
--- a/testsuite/python/mmm1d.py
+++ b/testsuite/python/mmm1d.py
@@ -35,7 +35,7 @@ class ElectrostaticInteractionsTests:
     system.cell_system.set_n_square()
     system.thermostat.set_langevin(kT=0, gamma=1, seed=8)
 
-    data = np.loadtxt(tests_common.abspath("data/mmm1d_data.txt"))
+    data = np.loadtxt(tests_common.data_path("mmm1d_data.txt"))
     p_pos = data[:, 1:4]
     p_q = data[:, 4]
     forces_target = data[:, 5:8]
@@ -46,36 +46,30 @@ class ElectrostaticInteractionsTests:
     def setUp(self):
         self.system.periodicity = [0, 0, 1]
         self.system.cell_system.set_n_square()
-        self.system.part.add(pos=self.p_pos, q=self.p_q)
-        self.mmm1d = self.MMM1D(prefactor=1.0, maxPWerror=1e-20)
-        self.system.actors.add(self.mmm1d)
-        self.system.integrator.run(steps=0)
 
     def tearDown(self):
         self.system.part.clear()
         self.system.actors.clear()
 
-    def test_forces(self):
+    def test_forces_and_energy(self):
+        self.system.part.add(pos=self.p_pos, q=self.p_q)
+        mmm1d = self.MMM1D(prefactor=1.0, maxPWerror=1e-20)
+        self.system.actors.add(mmm1d)
+        self.system.integrator.run(steps=0)
         measured_f = np.copy(self.system.part.all().f)
         np.testing.assert_allclose(measured_f, self.forces_target,
                                    atol=self.allowed_error)
-
-    def test_energy(self):
         measured_el_energy = self.system.analysis.energy()["coulomb"]
         self.assertAlmostEqual(
             measured_el_energy, self.energy_target, delta=self.allowed_error,
             msg="Measured energy deviates too much from stored result")
 
-    def test_with_analytical_result(self, prefactor=1.0, accuracy=1e-4):
-        self.system.part.clear()
-        p = self.system.part.add(pos=[0, 0, 0], q=1)
-        self.system.part.add(pos=[0, 0, 1], q=1)
-
-        self.system.integrator.run(steps=0)
+    def check_with_analytical_result(self, prefactor, accuracy):
+        p = self.system.part.by_id(0)
         f_measured = p.f
         energy_measured = self.system.analysis.energy()["total"]
-        target_energy_config = 1.00242505606 * prefactor
-        target_force_z_config = -0.99510759 * prefactor
+        target_energy_config = -1.00242505606 * prefactor
+        target_force_z_config = 0.99510759 * prefactor
 
         self.assertAlmostEqual(
             f_measured[0], 0, delta=self.allowed_error,
@@ -90,17 +84,23 @@ def test_with_analytical_result(self, prefactor=1.0, accuracy=1e-4):
             energy_measured, target_energy_config, delta=self.allowed_error,
             msg="Measured energy deviates too much from analytical result")
 
+    def test_with_analytical_result(self):
+        self.system.part.add(pos=[0, 0, 0], q=1)
+        self.system.part.add(pos=[0, 0, 1], q=-1)
+        mmm1d = self.MMM1D(prefactor=1.0, maxPWerror=1e-20)
+        self.system.actors.add(mmm1d)
+        self.system.integrator.run(steps=0)
+        self.check_with_analytical_result(prefactor=1.0, accuracy=0.0004)
+
     def test_bjerrum_length_change(self):
-        self.system.part.clear()
-        self.system.actors.clear()
-        prefactor = 2
-        mmm1d = self.MMM1D(prefactor=prefactor, maxPWerror=1e-20)
+        self.system.part.add(pos=[0, 0, 0], q=1)
+        self.system.part.add(pos=[0, 0, 1], q=-1)
+        mmm1d = self.MMM1D(prefactor=2.0, maxPWerror=1e-20)
         self.system.actors.add(mmm1d)
-        self.test_with_analytical_result(prefactor=prefactor, accuracy=0.0017)
+        self.system.integrator.run(steps=0)
+        self.check_with_analytical_result(prefactor=2.0, accuracy=0.0017)
 
     def test_exceptions(self):
-        self.system.actors.clear()
-        del self.mmm1d
         # check periodicity exceptions
         for periodicity in itertools.product(range(2), range(2), range(2)):
             if periodicity == (0, 0, 1):
diff --git a/testsuite/python/npt_thermostat.py b/testsuite/python/npt_thermostat.py
index 7a76573c661..83e27d3fece 100644
--- a/testsuite/python/npt_thermostat.py
+++ b/testsuite/python/npt_thermostat.py
@@ -117,7 +117,7 @@ def reset_particle_and_box():
     def test_02__direction(self):
         """Test for NpT constrained in one direction."""
 
-        data = np.genfromtxt(tests_common.abspath("data/npt_lj_system.data"))
+        data = np.genfromtxt(tests_common.data_path("npt_lj_system.data"))
         ref_box_l = 1.01 * np.max(data[:, 0:3])
 
         system = self.system
diff --git a/testsuite/python/oif_volume_conservation.py b/testsuite/python/oif_volume_conservation.py
index 7169b718280..5a74175235a 100644
--- a/testsuite/python/oif_volume_conservation.py
+++ b/testsuite/python/oif_volume_conservation.py
@@ -37,8 +37,9 @@ def test(self):
 
         # creating the template for OIF object
         cell_type = oif.OifCellType(
-            nodes_file=tests_common.abspath("data/sphere393nodes.dat"),
-            triangles_file=tests_common.abspath("data/sphere393triangles.dat"),
+            nodes_file=str(tests_common.data_path("sphere393nodes.dat")),
+            triangles_file=str(
+                tests_common.data_path("sphere393triangles.dat")),
             system=system, ks=1.0, kb=1.0, kal=1.0, kag=0.1, kv=0.1,
             check_orientation=False, resize=(3.0, 3.0, 3.0))
 
diff --git a/testsuite/python/p3m_electrostatic_pressure.py b/testsuite/python/p3m_electrostatic_pressure.py
index 26927ed9422..0e0cc065660 100644
--- a/testsuite/python/p3m_electrostatic_pressure.py
+++ b/testsuite/python/p3m_electrostatic_pressure.py
@@ -105,10 +105,14 @@ def setUp(self):
             max_displacement=0.01)
 
         # warmup
-        while self.system.analysis.energy()["total"] > 10 * num_part:
-            print("minimization: {:.1f}".format(
-                self.system.analysis.energy()["total"]))
+        energy = self.system.analysis.energy()["total"]
+        print(f"minimization: {energy:.1f}")
+        for _ in range(10):
             self.system.integrator.run(10)
+            energy = self.system.analysis.energy()["total"]
+            print(f"minimization: {energy:.1f}")
+            if energy < 2 * num_part:
+                break
         self.system.integrator.set_vv()
         self.system.thermostat.set_langevin(kT=self.kT, gamma=1.0, seed=41)
 
@@ -121,8 +125,9 @@ def test_p3m_pressure(self):
             cao=6,
             r_cut=1.4941e-01 * self.system.box_l[0])
         self.system.actors.add(p3m)
-        print("Tune skin: {}".format(self.system.cell_system.tune_skin(
-            min_skin=0.0, max_skin=2.5, tol=0.05, int_steps=100)))
+        skin = self.system.cell_system.tune_skin(
+            min_skin=0.0, max_skin=2.5, tol=0.05, int_steps=100)
+        print(f"Tuned skin: {skin}")
 
         num_samples = 25
         pressure_via_volume_scaling = pressureViaVolumeScaling(
diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py
index 9c762382b57..57ac7e6a0b6 100644
--- a/testsuite/python/save_checkpoint.py
+++ b/testsuite/python/save_checkpoint.py
@@ -17,7 +17,7 @@
 import unittest as ut
 import unittest_generator as utg
 import numpy as np
-import os
+import pathlib
 
 import espressomd
 import espressomd.checkpointing
@@ -52,10 +52,11 @@
 # create checkpoint folder
 checkpoint = espressomd.checkpointing.Checkpoint(
     **config.get_checkpoint_params())
+path_cpt_root = pathlib.Path(checkpoint.checkpoint_dir)
 
 # cleanup old checkpoint files
-for filepath in os.listdir(checkpoint.checkpoint_dir):
-    os.remove(os.path.join(checkpoint.checkpoint_dir, filepath))
+for filepath in path_cpt_root.iterdir():
+    filepath.unlink(missing_ok=True)
 
 n_nodes = system.cell_system.get_state()["n_nodes"]
 
@@ -245,7 +246,7 @@
     h5_units = espressomd.io.writer.h5md.UnitSystem(
         time="ps", mass="u", length="m", charge="e")
     h5 = espressomd.io.writer.h5md.H5md(
-        file_path=os.path.join(checkpoint.checkpoint_dir, "test.h5"),
+        file_path=str(path_cpt_root / "test.h5"),
         unit_system=h5_units)
     h5.write()
     h5.flush()
@@ -328,8 +329,8 @@
             for k in range(nz):
                 lbf[i, j, k].population = grid_3D[i, j, k] * np.arange(1, 20)
     # save LB checkpoint file
-    lbf_cpt_path = checkpoint.checkpoint_dir + "/lb.cpt"
-    lbf.save_checkpoint(lbf_cpt_path, lbf_cpt_mode)
+    lbf_cpt_path = path_cpt_root / "lb.cpt"
+    lbf.save_checkpoint(str(lbf_cpt_path), lbf_cpt_mode)
 
 # save checkpoint file
 checkpoint.save(0)
@@ -341,15 +342,15 @@ def test_checkpointing(self):
         '''
         Check for the presence of the checkpoint files.
         '''
-        self.assertTrue(os.path.isdir(checkpoint.checkpoint_dir),
+        self.assertTrue(path_cpt_root.is_dir(),
                         "checkpoint directory not created")
 
-        checkpoint_filepath = checkpoint.checkpoint_dir + "/0.checkpoint"
-        self.assertTrue(os.path.isfile(checkpoint_filepath),
+        checkpoint_filepath = path_cpt_root / "0.checkpoint"
+        self.assertTrue(checkpoint_filepath.is_file(),
                         "checkpoint file not created")
 
         if lbf_actor:
-            self.assertTrue(os.path.isfile(lbf_cpt_path),
+            self.assertTrue(lbf_cpt_path.is_file(),
                             "LB checkpoint file not created")
 
     @ut.skipIf(lbf_actor is None, "Skipping test due to missing mode.")
@@ -361,25 +362,23 @@ def test_lb_checkpointing_exceptions(self):
 
         # check exception mechanism
         with self.assertRaisesRegex(RuntimeError, 'could not open file'):
-            dirname, filename = os.path.split(lbf_cpt_path)
-            invalid_path = os.path.join(dirname, 'unknown_dir', filename)
-            lbf.save_checkpoint(invalid_path, lbf_cpt_mode)
+            invalid_path = lbf_cpt_path.parent / "unknown_dir" / "lb.cpt"
+            lbf.save_checkpoint(str(invalid_path), lbf_cpt_mode)
         system.actors.remove(lbf)
         with self.assertRaisesRegex(RuntimeError, 'one needs to have already initialized the LB fluid'):
-            lbf.load_checkpoint(lbf_cpt_path, lbf_cpt_mode)
+            lbf.load_checkpoint(str(lbf_cpt_path), lbf_cpt_mode)
 
         # read the valid LB checkpoint file
-        with open(lbf_cpt_path, "rb") as f:
-            lbf_cpt_str = f.read()
-        cpt_path = checkpoint.checkpoint_dir + "/lb{}.cpt"
+        lbf_cpt_data = lbf_cpt_path.read_bytes()
+        cpt_path = str(path_cpt_root / "lb") + "{}.cpt"
         # write checkpoint file with missing data
         with open(cpt_path.format("-missing-data"), "wb") as f:
-            f.write(lbf_cpt_str[:len(lbf_cpt_str) // 2])
+            f.write(lbf_cpt_data[:len(lbf_cpt_data) // 2])
         # write checkpoint file with extra data
         with open(cpt_path.format("-extra-data"), "wb") as f:
-            f.write(lbf_cpt_str + lbf_cpt_str[-8:])
+            f.write(lbf_cpt_data + lbf_cpt_data[-8:])
         if lbf_cpt_mode == 0:
-            boxsize, data = lbf_cpt_str.split(b"\n", 1)
+            boxsize, data = lbf_cpt_data.split(b"\n", 1)
             # write checkpoint file with incorrectly formatted data
             with open(cpt_path.format("-wrong-format"), "wb") as f:
                 f.write(boxsize + b"\ntext string\n" + data)
diff --git a/testsuite/python/scafacos_dipoles_1d_2d.py b/testsuite/python/scafacos_dipoles_1d_2d.py
index 9d1cde69274..b32b6a30aa7 100644
--- a/testsuite/python/scafacos_dipoles_1d_2d.py
+++ b/testsuite/python/scafacos_dipoles_1d_2d.py
@@ -62,19 +62,19 @@ def test_scafacos(self):
 
             # Read reference data
             if dim == 2:
-                file_prefix = "data/mdlc"
+                file_prefix = "mdlc"
                 s.periodicity = [1, 1, 0]
             else:
                 s.periodicity = [1, 0, 0]
-                file_prefix = "data/scafacos_dipoles_1d"
+                file_prefix = "scafacos_dipoles_1d"
 
-            ref_E_path = tests_common.abspath(
-                file_prefix + "_reference_data_energy.dat")
+            ref_E_path = tests_common.data_path(
+                f"{file_prefix}_reference_data_energy.dat")
             ref_E = float(np.genfromtxt(ref_E_path))
 
             # Particles
-            data = np.genfromtxt(tests_common.abspath(
-                file_prefix + "_reference_data_forces_torques.dat"))
+            data = np.genfromtxt(tests_common.data_path(
+                f"{file_prefix}_reference_data_forces_torques.dat"))
             s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
             s.part.all().rotation = 3 * [True]
 
diff --git a/testsuite/python/script_interface.py b/testsuite/python/script_interface.py
index 9940654e417..b762063db21 100644
--- a/testsuite/python/script_interface.py
+++ b/testsuite/python/script_interface.py
@@ -102,11 +102,9 @@ def test_variant_exceptions(self):
         constraint = espressomd.constraints.ShapeBasedConstraint()
         # check conversion of unsupported types
         with self.assertRaisesRegex(TypeError, "No conversion from type 'module' to 'Variant'"):
-            espressomd.constraints.ShapeBasedConstraint(unknown=ut)
+            espressomd.constraints.ShapeBasedConstraint(shape=ut)
         with self.assertRaisesRegex(TypeError, "No conversion from type 'module' to 'Variant'"):
             constraint.set_params(shape=ut)
-        with self.assertRaisesRegex(TypeError, "No conversion from type 'module' to 'Variant'"):
-            constraint.call_method('unknown', unknown=ut)
         # check restrictions on the dict type
         with self.assertRaisesRegex(TypeError, r"No conversion from type 'dict_item\(\[\(str, int\)\]\)' to 'Variant\[std::(__1::)?unordered_map<int, Variant>\]'"):
             constraint.shape = {'1': 2}
diff --git a/testsuite/python/sigint.py b/testsuite/python/sigint.py
index 1627fba0633..17fc6eed6fa 100644
--- a/testsuite/python/sigint.py
+++ b/testsuite/python/sigint.py
@@ -21,13 +21,13 @@
 import subprocess
 import time
 import sys
-import os
+import pathlib
 
 
 class SigintTest(ut.TestCase):
 
     def setUp(self):
-        script = os.path.join(os.path.dirname(__file__), 'sigint_child.py')
+        script = str(pathlib.Path(__file__).parent / 'sigint_child.py')
         self.process = subprocess.Popen([sys.executable, script])
 
     def test_signal_handling(self):
diff --git a/testsuite/python/stokesian_dynamics.py b/testsuite/python/stokesian_dynamics.py
index 0712455d351..deb95245844 100644
--- a/testsuite/python/stokesian_dynamics.py
+++ b/testsuite/python/stokesian_dynamics.py
@@ -35,7 +35,7 @@ class StokesianDynamicsTest(ut.TestCase):
     # Digitized reference data of Figure 5b from
     # Durlofsky et al., J. Fluid Mech. 180, 21 (1987)
     # https://doi.org/10.1017/S002211208700171X
-    data = np.loadtxt(tests_common.abspath('data/dancing.txt'))
+    data = np.loadtxt(tests_common.data_path('dancing.txt'))
 
     def setUp(self):
         self.system.box_l = [10] * 3
diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py
index 5012ff24da3..b298de69af0 100644
--- a/testsuite/python/test_checkpoint.py
+++ b/testsuite/python/test_checkpoint.py
@@ -20,6 +20,8 @@
 import unittest_decorators as utx
 import unittest_generator as utg
 import numpy as np
+import contextlib
+import pathlib
 
 import espressomd
 import espressomd.checkpointing
@@ -32,11 +34,8 @@
 import espressomd.shapes
 import espressomd.constraints
 
-import os
-try:
+with contextlib.suppress(ImportError):
     import h5py  # h5py has to be imported *after* espressomd (MPI)
-except ImportError:
-    h5py = None
 
 config = utg.TestGenerator()
 is_gpu_available = espressomd.gpu_available()
@@ -50,6 +49,7 @@ class CheckpointTest(ut.TestCase):
     checkpoint = espressomd.checkpointing.Checkpoint(
         **config.get_checkpoint_params())
     checkpoint.load(0)
+    path_cpt_root = pathlib.Path(checkpoint.checkpoint_dir)
     n_nodes = system.cell_system.get_state()["n_nodes"]
 
     @classmethod
@@ -69,6 +69,13 @@ def get_active_actor_of_type(self, actor_type):
         self.fail(
             f"system doesn't have an actor of type {actor_type.__name__}")
 
+    def test_get_active_actor_of_type(self):
+        if system.actors.active_actors:
+            actor = system.actors.active_actors[0]
+            self.assertEqual(self.get_active_actor_of_type(type(actor)), actor)
+        with self.assertRaisesRegex(AssertionError, "system doesn't have an actor of type Wall"):
+            self.get_active_actor_of_type(espressomd.shapes.Wall)
+
     @ut.skipIf(not has_lb_mode, "Skipping test due to missing mode.")
     def test_lb_fluid(self):
         '''
@@ -81,7 +88,8 @@ def test_lb_fluid(self):
         lbf = self.get_active_actor_of_type(
             espressomd.lb.HydrodynamicInteraction)
         cpt_mode = 0 if 'LB.ASCII' in modes else 1
-        cpt_path = self.checkpoint.checkpoint_dir + "/lb{}.cpt"
+        cpt_root = pathlib.Path(self.checkpoint.checkpoint_dir)
+        cpt_path = str(cpt_root / "lb") + "{}.cpt"
 
         # check exception mechanism with corrupted LB checkpoint files
         with self.assertRaisesRegex(RuntimeError, 'EOF found'):
@@ -415,16 +423,15 @@ def test_correlator(self):
                 expected)
 
     @utx.skipIfMissingFeatures('H5MD')
-    @ut.skipIf(h5py is None,
-               "Skipping test due to missing python module 'h5py'.")
+    @utx.skipIfMissingModules("h5py")
     def test_h5md(self):
         # check attributes
-        file_path = os.path.join(self.checkpoint.checkpoint_dir, "test.h5")
-        script_path = os.path.join(
-            os.path.dirname(__file__), "save_checkpoint.py")
+        file_path = self.path_cpt_root / "test.h5"
+        script_path = pathlib.Path(
+            __file__).resolve().parent / "save_checkpoint.py"
         self.assertEqual(h5.fields, ['all'])
-        self.assertEqual(h5.script_path, script_path)
-        self.assertEqual(h5.file_path, file_path)
+        self.assertEqual(h5.script_path, str(script_path))
+        self.assertEqual(h5.file_path, str(file_path))
 
         # write new frame
         h5.write()
diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py
index f57327f01ca..d8a7b80437c 100644
--- a/testsuite/python/tests_common.py
+++ b/testsuite/python/tests_common.py
@@ -14,7 +14,7 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-import os
+import pathlib
 import numpy as np
 
 
@@ -128,8 +128,11 @@ def verify_lj_forces(system, tolerance, ids_to_skip=()):
             err_msg=f"LJ force verification failed on particle {p.id}.")
 
 
-def abspath(path):
-    return os.path.join(os.path.dirname(os.path.abspath(__file__)), path)
+def data_path(filename):
+    """
+    Resolve abosulte path to a resource.
+    """
+    return pathlib.Path(__file__).resolve().parent / "data" / filename
 
 
 def transform_pos_from_cartesian_to_polar_coordinates(pos):
diff --git a/testsuite/python/unittest_generator.py b/testsuite/python/unittest_generator.py
index 6e16257d772..b399f8eac28 100644
--- a/testsuite/python/unittest_generator.py
+++ b/testsuite/python/unittest_generator.py
@@ -17,9 +17,9 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
-import os
 import sys
 import inspect
+import pathlib
 
 
 class TestGenerator:
@@ -141,4 +141,4 @@ def get_checkpoint_params(self):
         Generate parameters to instantiate an ESPResSo checkpoint file.
         """
         return {"checkpoint_id": f"checkpoint_{self.test_idx}",
-                "checkpoint_path": os.path.dirname(__file__)}
+                "checkpoint_path": str(pathlib.Path(__file__).parent)}
diff --git a/testsuite/scripts/importlib_wrapper.py b/testsuite/scripts/importlib_wrapper.py
index 87ec2229749..f7377ead943 100644
--- a/testsuite/scripts/importlib_wrapper.py
+++ b/testsuite/scripts/importlib_wrapper.py
@@ -22,13 +22,10 @@
 import ast
 import tokenize
 import unittest
+import unittest.mock
 import importlib
+import pathlib
 import espressomd
-from unittest.mock import MagicMock
-
-
-def _id(x):
-    return x
 
 
 # global variable: if one import failed, all subsequent imports will be skipped,
@@ -40,7 +37,7 @@ def configure_and_import(filepath,
                          gpu=False,
                          substitutions=lambda x: x,
                          cmd_arguments=None,
-                         script_suffix=None,
+                         script_suffix="",
                          move_to_script_dir=True,
                          mock_visualizers=True,
                          **parameters):
@@ -55,23 +52,23 @@ def configure_and_import(filepath,
 
     Parameters
     ----------
-    filepath : str
+    filepath : :obj:`str`
         python script to import
-    gpu : bool
+    gpu : :obj:`bool`
         whether GPU is necessary or not
-    substitutions : function
+    substitutions : :obj:`function`
         custom text replacement operation (useful to edit out calls to the
         OpenGL or Mayavi visualizers' ``run()`` method)
-    cmd_arguments : list
+    cmd_arguments : :obj:`list`
         command line arguments, i.e. sys.argv without the script path
-    script_suffix : str
+    script_suffix : :obj:`str`
         suffix to append to the configured script (useful when a single
         module is being tested by multiple tests in parallel)
-    mock_visualizers : bool
-        if ``True``, substitute ES visualizers with `Mock()` classes in case
-        of `ImportError()` (use ``False`` if an `ImportError()` is relevant
+    mock_visualizers : :obj:`bool`
+        if ``True``, substitute ES visualizers with ``Mock`` classes in case
+        of ``ImportError`` (use ``False`` if an ``ImportError`` is relevant
         to your test)
-    move_to_script_dir : bool
+    move_to_script_dir : :obj:`bool`
         if ``True``, move to the script's directory (useful when the script
         needs to load files hardcoded as relative paths, or when files are
         generated and need cleanup); this is enabled by default
@@ -79,22 +76,18 @@ def configure_and_import(filepath,
         global variables to replace
 
     """
+    filepath = pathlib.Path(filepath).resolve()
     if skip_future_imports:
-        module = MagicMock()
+        module = unittest.mock.MagicMock()
         skipIfMissingImport = skip_future_imports_dependency(filepath)
         return module, skipIfMissingImport
     if gpu and not espressomd.gpu_available():
         skip_future_imports_dependency(filepath)
         skipIfMissingGPU = unittest.skip("gpu not available, skipping test!")
-        module = MagicMock()
+        module = unittest.mock.MagicMock()
         return module, skipIfMissingGPU
-    filepath = os.path.abspath(filepath)
     # load original script
-    # read in binary mode, then decode as UTF-8 to avoid this python3.5 error:
-    # UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 915:
-    #                      ordinal not in range(128)
-    with open(filepath, "rb") as f:
-        code = f.read().decode(encoding="utf-8")
+    code = filepath.read_text()
     # custom substitutions
     code = substitutions(code)
     assert code.strip()
@@ -109,31 +102,25 @@ def configure_and_import(filepath,
     if mock_visualizers:
         code = mock_es_visualization(code)
     # save changes to a new file
-    if script_suffix:
-        if script_suffix[0] != "_":
-            script_suffix = "_" + script_suffix
-    else:
-        script_suffix = ""
-    script_suffix += "_processed.py"
-    output_filepath = os.path.splitext(filepath)[0] + script_suffix
-    assert os.path.isfile(output_filepath) is False, \
+    output_filepath = filepath.parent / \
+        f"{filepath.stem}_{script_suffix}_processed.py"
+    assert not output_filepath.exists(), \
         f"File {output_filepath} already processed, cannot overwrite"
-    with open(output_filepath, "wb") as f:
-        f.write(code.encode(encoding="utf-8"))
+    output_filepath.write_bytes(code.encode(encoding="utf-8"))
     # import
-    dirname, basename = os.path.split(output_filepath)
+    dirname = output_filepath.parent
     if move_to_script_dir:
         os.chdir(dirname)
-    sys.path.insert(0, dirname)
-    module_name = os.path.splitext(basename)[0]
+    sys.path.insert(0, str(dirname))
+    module_name = output_filepath.stem
     try:
         module = importlib.import_module(module_name)
     except espressomd.FeaturesError as err:
         skip_future_imports_dependency(filepath)
         skipIfMissingFeatures = unittest.skip(f"{err}, skipping test!")
-        module = MagicMock()
+        module = unittest.mock.MagicMock()
     else:
-        skipIfMissingFeatures = _id
+        def skipIfMissingFeatures(x): return x
     if cmd_arguments is not None:
         # restore original command line arguments
         sys.argv = old_sys_argv
@@ -165,7 +152,7 @@ def set_cmd(code, filepath, cmd_arguments):
     """
     assert isinstance(cmd_arguments, (list, tuple))
     sys_argv = list(map(str, cmd_arguments))
-    sys_argv.insert(0, os.path.basename(filepath))
+    sys_argv.insert(0, filepath.name)
     visitor = GetSysArgparseImports()
     visitor.visit(ast.parse(protect_ipython_magics(code)))
     assert visitor.linenos, "module sys (or argparse) is not imported"
@@ -216,12 +203,12 @@ def substitute_variable_values(code, strings_as_is=False, keep_original=True,
 
     Parameters
     ----------
-    code : str
+    code : :obj:`str`
         Source code to edit.
-    strings_as_is : bool
+    strings_as_is : :obj:`bool`
         If ``True``, consider all values in \*\*parameters are strings and
         substitute them in-place without formatting by ``repr()``.
-    keep_original : bool
+    keep_original : :obj:`bool`
         Keep the original value (e.g. ``N = 10; _N__original = 1000``), helps
         with debugging.
     \*\*parameters :
@@ -245,8 +232,8 @@ def substitute_variable_values(code, strings_as_is=False, keep_original=True,
             if keep_original:
                 lines[lineno - 1] += "; _" + varname + "__original" + old_value
             else:
-                for lineno in range(lineno + 1, mapping[lineno]):
-                    lines[lineno - 1] = ""
+                for lineno in range(lineno, mapping[lineno]):
+                    lines[lineno] = ""
     return "\n".join(lines)
 
 
@@ -582,9 +569,9 @@ def visit_ImportFrom(self, node):
 
 def mock_es_visualization(code):
     """
-    Replace ``import espressomd.visualization_<backend>`` by a ``MagicMock()``
+    Replace ``import espressomd.visualization_<backend>`` by a ``MagicMock``
     when the visualization module is unavailable, by catching the
-    ``ImportError()`` exception. Please note that ``espressomd.visualization``
+    ``ImportError`` exception. Please note that ``espressomd.visualization``
     is deferring the exception, thus requiring additional checks.
 
     Import aliases are supported, however please don't use
@@ -596,9 +583,9 @@ def mock_es_visualization(code):
 try:
     {0}{1}
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    {2} = MagicMock()
+    {2} = unittest.mock.MagicMock()
 """.lstrip()
 
     def check_for_deferred_ImportError(line, alias):
@@ -639,7 +626,7 @@ def skip_future_imports_dependency(filepath):
     """
     global skip_future_imports
     if not skip_future_imports:
-        module_name = os.path.splitext(os.path.basename(filepath))[0]
+        module_name = filepath.stem
         assert module_name != ""
         skip_future_imports = module_name
     return unittest.skip(
diff --git a/testsuite/scripts/samples/test_ekboundaries.py b/testsuite/scripts/samples/test_ekboundaries.py
index 0b2eecf914b..ab10be10f2e 100644
--- a/testsuite/scripts/samples/test_ekboundaries.py
+++ b/testsuite/scripts/samples/test_ekboundaries.py
@@ -17,7 +17,7 @@
 
 import unittest as ut
 import importlib_wrapper
-import os
+import pathlib
 
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@SAMPLES_DIR@/ekboundaries.py", gpu=True, n_int_cycles=50)
@@ -28,13 +28,12 @@ class Sample(ut.TestCase):
     system = sample.system
 
     def test_file_generation(self):
-        # test .checkpoint files exist
+        # test .vtk files exist
+        path_vtk_root = pathlib.Path("ek")
         for basename in ["pos_dens_0.vtk", "pos_flux_0.vtk", "ekv_0.vtk",
                          "neg_dens_0.vtk", "neg_flux_0.vtk", "ekb_0.vtk"]:
-            filepath = os.path.join("ek", basename)
-            self.assertTrue(
-                os.path.isfile(filepath),
-                filepath + " not created")
+            filepath = path_vtk_root / basename
+            self.assertTrue(filepath.is_file(), f"File {filepath} not created")
 
 
 if __name__ == "__main__":
diff --git a/testsuite/scripts/samples/test_immersed_boundary.py b/testsuite/scripts/samples/test_immersed_boundary.py
index e2b1aae57fe..7c32d17b294 100644
--- a/testsuite/scripts/samples/test_immersed_boundary.py
+++ b/testsuite/scripts/samples/test_immersed_boundary.py
@@ -17,7 +17,7 @@
 
 import unittest as ut
 import importlib_wrapper
-import os
+import pathlib
 
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@SAMPLES_DIR@/immersed_boundary/sampleImmersedBoundary.py",
@@ -31,10 +31,8 @@ class Sample(ut.TestCase):
     def test_file_generation(self):
         # test .vtk files exist
         for i in range(sample.numSteps + 1):
-            filepath = os.path.join("outputVolParaCUDA", "cell_%i.vtk" % i)
-            self.assertTrue(
-                os.path.isfile(filepath),
-                filepath + " not created")
+            filepath = pathlib.Path("outputVolParaCUDA") / f"cell_{i}.vtk"
+            self.assertTrue(filepath.is_file(), f"File {filepath} not created")
 
 
 if __name__ == "__main__":
diff --git a/testsuite/scripts/samples/test_object_in_fluid__motivation.py b/testsuite/scripts/samples/test_object_in_fluid__motivation.py
index 558d57941d5..191346001e3 100644
--- a/testsuite/scripts/samples/test_object_in_fluid__motivation.py
+++ b/testsuite/scripts/samples/test_object_in_fluid__motivation.py
@@ -18,6 +18,7 @@
 import unittest as ut
 import importlib_wrapper
 import os
+import pathlib
 
 os.chdir("@SAMPLES_DIR@/object_in_fluid")
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
@@ -39,13 +40,13 @@ def test_file_generation(self):
             "wallFront.vtk"]
         for i in [0, 1]:
             for j in range(sample.maxCycle):
-                basenames.append("cell{}_{}.vtk".format(i, j))
+                basenames.append(f"cell{i}_{j}.vtk")
+
         # test .vtk files exist
+        path_vtk_root = pathlib.Path("output")
         for name in basenames:
-            filepath = os.path.join("output", "sim0", name)
-            self.assertTrue(
-                os.path.isfile(filepath),
-                filepath + " not created")
+            filepath = path_vtk_root / "sim0" / name
+            self.assertTrue(filepath.is_file(), f"File {filepath} not created")
 
 
 if __name__ == "__main__":
diff --git a/testsuite/scripts/samples/test_save_checkpoint.py b/testsuite/scripts/samples/test_save_checkpoint.py
index 55b61ac5405..b06e265ee41 100644
--- a/testsuite/scripts/samples/test_save_checkpoint.py
+++ b/testsuite/scripts/samples/test_save_checkpoint.py
@@ -17,7 +17,7 @@
 
 import unittest as ut
 import importlib_wrapper
-import os
+import pathlib
 
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@SAMPLES_DIR@/save_checkpoint.py")
@@ -29,8 +29,8 @@ class Sample(ut.TestCase):
 
     def test_file_generation(self):
         # test .checkpoint files exist
-        filepath = os.path.join("mycheckpoint", "0.checkpoint")
-        self.assertTrue(os.path.isfile(filepath), filepath + " not created")
+        filepath = pathlib.Path("mycheckpoint") / "0.checkpoint"
+        self.assertTrue(filepath.exists(), f"File {filepath} not created")
 
 
 if __name__ == "__main__":
diff --git a/testsuite/scripts/test_importlib_wrapper.py b/testsuite/scripts/test_importlib_wrapper.py
index e21904ecb6f..67d800a450c 100644
--- a/testsuite/scripts/test_importlib_wrapper.py
+++ b/testsuite/scripts/test_importlib_wrapper.py
@@ -19,21 +19,25 @@
 import importlib_wrapper as iw
 import sys
 import ast
+import pathlib
+import tempfile
+import espressomd
 
 
 class importlib_wrapper(ut.TestCase):
 
     def test_substitute_variable_values(self):
-        str_inp = "n_steps=5000\nnsteps == 5\n"
-        str_exp = "n_steps = 10; _n_steps__original=5000\nnsteps == 5\n"
+        str_inp = "n_steps=(\n5000)\nnsteps == 5\n"
+        str_exp = "n_steps = 10; _n_steps__original=(\n5000)\nnsteps == 5\n"
         str_out = iw.substitute_variable_values(str_inp, n_steps=10)
         self.assertEqual(str_out, str_exp)
         str_out = iw.substitute_variable_values(str_inp, n_steps='10',
                                                 strings_as_is=True)
         self.assertEqual(str_out, str_exp)
-        str_inp = "N=5000\nnsteps == 5\n"
-        str_exp = "N = 10\nnsteps == 5\n"
-        str_out = iw.substitute_variable_values(str_inp, N=10, keep_original=0)
+        str_inp = "N=(\n5000)\nnsteps == 5\n"
+        str_exp = "N = 10\n\nnsteps == 5\n"
+        str_out = iw.substitute_variable_values(str_inp, N=10,
+                                                keep_original=False)
         self.assertEqual(str_out, str_exp)
         # test exceptions
         str_inp = "n_steps=5000\nnsteps == 5\n"
@@ -49,20 +53,21 @@ def test_substitute_variable_values(self):
     def test_set_cmd(self):
         original_sys_argv = list(sys.argv)
         sys.argv = [0, "test"]
+        path = pathlib.Path("a.py")
         # test substitutions
         str_inp = "import sys\nimport argparse"
-        str_exp = "import sys;sys.argv = ['a.py', '1', '2'];" + str_inp
-        str_out, sys_argv = iw.set_cmd(str_inp, "a.py", (1, 2))
+        str_exp = f"import sys;sys.argv = ['{path.name}', '1', '2'];" + str_inp
+        str_out, sys_argv = iw.set_cmd(str_inp, path, (1, 2))
         self.assertEqual(str_out, str_exp)
         self.assertEqual(sys_argv, [0, "test"])
         str_inp = "import argparse"
-        str_exp = "import sys;sys.argv = ['a.py', '1', '2'];" + str_inp
-        str_out, sys_argv = iw.set_cmd(str_inp, "a.py", ["1", 2])
+        str_exp = f"import sys;sys.argv = ['{path.name}', '1', '2'];" + str_inp
+        str_out, sys_argv = iw.set_cmd(str_inp, path, ["1", 2])
         self.assertEqual(str_out, str_exp)
         self.assertEqual(sys_argv, [0, "test"])
         # test exceptions
         str_inp = "import re"
-        self.assertRaises(AssertionError, iw.set_cmd, str_inp, "a.py", (1, 2))
+        self.assertRaises(AssertionError, iw.set_cmd, str_inp, path, (1, 2))
         # restore sys.argv
         sys.argv = original_sys_argv
 
@@ -91,9 +96,9 @@ def test_mock_es_visualization(self):
        hasattr(espressomd.visualization.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    espressomd.visualization = MagicMock()
+    espressomd.visualization = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -105,9 +110,9 @@ def test_mock_es_visualization(self):
        hasattr(test.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -119,18 +124,18 @@ def test_mock_es_visualization(self):
        hasattr(espressomd.visualization.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    espressomd.visualization = MagicMock()
+    espressomd.visualization = unittest.mock.MagicMock()
 try:
     import espressomd.visualization as test
     if hasattr(test.mayaviLive, 'deferred_ImportError') or \\
        hasattr(test.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -142,9 +147,9 @@ def test_mock_es_visualization(self):
        hasattr(visualization.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    visualization = MagicMock()
+    visualization = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -156,9 +161,9 @@ def test_mock_es_visualization(self):
        hasattr(test.openGLLive, 'deferred_ImportError'):
         raise ImportError()
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -167,9 +172,9 @@ def test_mock_es_visualization(self):
 try:
     from espressomd import visualization_mayavi
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    visualization_mayavi = MagicMock()
+    visualization_mayavi = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -178,9 +183,9 @@ def test_mock_es_visualization(self):
 try:
     from espressomd import visualization_mayavi as test
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -189,9 +194,9 @@ def test_mock_es_visualization(self):
 try:
     from espressomd.visualization_mayavi import mayaviLive
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    mayaviLive = MagicMock()
+    mayaviLive = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -200,9 +205,9 @@ def test_mock_es_visualization(self):
 try:
     from espressomd.visualization_mayavi import mayaviLive as test
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -211,15 +216,15 @@ def test_mock_es_visualization(self):
 try:
     from espressomd.visualization_mayavi import a as b
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    b = MagicMock()
+    b = unittest.mock.MagicMock()
 try:
     from espressomd.visualization_mayavi import mayaviLive
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    mayaviLive = MagicMock()
+    mayaviLive = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -230,9 +235,9 @@ def test_mock_es_visualization(self):
     if hasattr(openGLLive, 'deferred_ImportError'):
         raise openGLLive.deferred_ImportError
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    openGLLive = MagicMock()
+    openGLLive = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -243,9 +248,9 @@ def test_mock_es_visualization(self):
     if hasattr(test, 'deferred_ImportError'):
         raise test.deferred_ImportError
 except ImportError:
-    from unittest.mock import MagicMock
+    import unittest.mock
     import espressomd
-    test = MagicMock()
+    test = unittest.mock.MagicMock()
 """
         self.assertEqual(iw.mock_es_visualization(statement), expected[1:])
 
@@ -311,6 +316,7 @@ def test_prng_seed_espressomd_system_visitor(self):
             'import espressomd.System as s1, espressomd.system.System as s2',
             'from espressomd import System as s3, electrostatics',
             'from espressomd.system import System as s4',
+            'from espressomd import system as es5',
             'sys1 = es1.System()',
             'sys2 = es1.system.System()',
             'sys3 = es2.System()',
@@ -318,6 +324,7 @@ def test_prng_seed_espressomd_system_visitor(self):
             'sys5 = s2()',
             'sys6 = s3()',
             'sys7 = s4()',
+            'sys8 = es5.System()',
             'import numpy as np',
             'import numpy.random as npr1',
             'from numpy import random as npr2',
@@ -330,13 +337,14 @@ def test_prng_seed_espressomd_system_visitor(self):
         v.visit(tree)
         # find all aliases for espressomd.system.System
         expected_es_sys_aliases = {'es1.System', 'es1.system.System',
-                                   'es2.System', 's1', 's2', 's3', 's4'}
+                                   'es2.System', 's1', 's2', 's3', 's4',
+                                   'es5.System'}
         self.assertEqual(v.es_system_aliases, expected_es_sys_aliases)
         # find all variables of type espressomd.system.System
-        expected_es_sys_objs = set(f'sys{i}' for i in range(1, 8))
+        expected_es_sys_objs = set(f'sys{i}' for i in range(1, 9))
         self.assertEqual(v.variable_system_aliases, expected_es_sys_objs)
         # find all seeds setup
-        self.assertEqual(v.numpy_seeds, [17, 18, 19])
+        self.assertEqual(v.numpy_seeds, [19, 20, 21])
         # test exceptions
         str_es_sys_list = [
             'import espressomd.System',
@@ -395,6 +403,49 @@ def test_ipython_magics(self):
         code_deprotected = iw.deprotect_ipython_magics(code_protected)
         self.assertEqual(code_deprotected, code)
 
+    def test_configure_and_import(self):
+        with tempfile.TemporaryDirectory() as temp_dir:
+            path_in = pathlib.Path(temp_dir) / "sample.py"
+            path_out = pathlib.Path(temp_dir) / "sample_test_processed.py"
+            path_features = pathlib.Path(temp_dir) / "sample_impossible.py"
+
+            # test importing a simple sample
+            sys.argv.append("42")
+            sys_argv_ref = list(sys.argv)
+            path_in.write_text("""
+import sys
+from argparse import ArgumentParser
+value = 42
+argv = list(sys.argv)
+import espressomd.visualization_opengl
+""")
+            sample, _ = iw.configure_and_import(
+                path_in, move_to_script_dir=False, cmd_arguments=["TestCase"],
+                gpu=False, script_suffix="test", value=43)
+            self.assertEqual(sys.argv, sys_argv_ref)
+            self.assertEqual(sample.argv, [path_in.name, "TestCase"])
+            self.assertTrue(path_out.exists(), f"File {path_out} not found")
+            self.assertEqual(sample.value, 43)
+            self.assertIn(
+                "espressomd.visualization_opengl = unittest.mock.MagicMock()",
+                path_out.read_text())
+
+            # test importing a sample that relies on features not compiled in
+            inactive_features = espressomd.all_features() - set(espressomd.features())
+            if inactive_features:
+                path_features.write_text(f"""
+import espressomd
+espressomd.assert_features({list(inactive_features)})
+""")
+                module, _ = iw.configure_and_import(path_features)
+                self.assertIsInstance(module, ut.mock.MagicMock)
+
+                # importing a valid sample is impossible after a failed import
+                sample, _ = iw.configure_and_import(
+                    path_in, move_to_script_dir=True, script_suffix="test",
+                    cmd_arguments=["TestCase"], gpu=False, value=43)
+                self.assertIsInstance(sample, ut.mock.MagicMock)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/scripts/tutorials/test_convert.py b/testsuite/scripts/tutorials/test_convert.py
index 0f24519cdc2..b35aab8482f 100644
--- a/testsuite/scripts/tutorials/test_convert.py
+++ b/testsuite/scripts/tutorials/test_convert.py
@@ -16,10 +16,10 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 import unittest as ut
-import os
 import sys
 import nbformat
 import traceback
+import pathlib
 
 sys.path.insert(0, '@CMAKE_BINARY_DIR@/doc/tutorials')
 import convert
@@ -96,17 +96,16 @@ def run_command(self, cmd, output=None):
             traceback.print_exc()
             self.fail(error_msg)
         if output is not None:
-            self.assertTrue(os.path.isfile(output), f"{output} not created")
+            self.assertTrue(output.exists(), f"File {output} not created")
 
     def test_html_wrapper(self):
-        f_input = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_notebook.ipynb'
-        f_output = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_notebook.run.ipynb'
-        f_script = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_script.py'
+        root = pathlib.Path("@CMAKE_CURRENT_BINARY_DIR@")
+        f_input = root / "test_convert_notebook.ipynb"
+        f_output = root / "test_convert_notebook.run.ipynb"
+        f_script = root / "test_convert_script.py"
         # setup
-        if os.path.isfile(f_output):
-            os.remove(f_output)
-        with open(f_script, 'w') as f:
-            f.write('global_var = 5')
+        f_output.unlink(missing_ok=True)
+        f_script.write_text('global_var = 5')
         with open(f_input, 'w', encoding='utf-8') as f:
             nb = nbformat.v4.new_notebook(metadata=self.nb_metadata)
             cell_md = nbformat.v4.new_markdown_cell(source=self.cell_md_src)
@@ -116,9 +115,9 @@ def test_html_wrapper(self):
             nbformat.write(nb, f)
         # run command and check for errors
         cmd = ['ci',
-               '--input', f_input,
-               '--output', f_output,
-               '--scripts', f_script,
+               '--input', str(f_input),
+               '--output', str(f_output),
+               '--scripts', str(f_script),
                '--substitutions', 'global_var=20',
                '--execute']
         self.run_command(cmd, f_output)
@@ -155,11 +154,11 @@ def test_html_wrapper(self):
         self.assertEqual(nb_output['cells'][4]['source'], 'global_var = 20')
 
     def test_exercise2_plugin(self):
-        f_input = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_exercise2.ipynb'
-        f_output = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_exercise2.run.ipynb'
+        root = pathlib.Path("@CMAKE_CURRENT_BINARY_DIR@")
+        f_input = root / "test_convert_exercise2.ipynb"
+        f_output = root / "test_convert_exercise2.run.ipynb"
         # setup
-        if os.path.isfile(f_output):
-            os.remove(f_output)
+        f_output.unlink(missing_ok=True)
         with open(f_input, 'w', encoding='utf-8') as f:
             nb = nbformat.v4.new_notebook(metadata=self.nb_metadata)
             # question with 2 answers and an empty cell
@@ -191,8 +190,8 @@ def test_exercise2_plugin(self):
             nbformat.write(nb, f)
         # run command and check for errors
         cmd = ['ci',
-               '--input', f_input,
-               '--output', f_output,
+               '--input', str(f_input),
+               '--output', str(f_output),
                '--substitutions', 'global_var=20',
                '--exercise2', '--remove-empty-cells']
         self.run_command(cmd, f_output)
@@ -226,7 +225,8 @@ def test_exercise2_plugin(self):
         self.assertEqual(next(cells, 'EOF'), 'EOF')
 
     def test_exercise2_conversion(self):
-        f_input = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_exercise2_conversion.ipynb'
+        root = pathlib.Path("@CMAKE_CURRENT_BINARY_DIR@")
+        f_input = root / "test_convert_exercise2_conversion.ipynb"
         # setup
         with open(f_input, 'w', encoding='utf-8') as f:
             nb = nbformat.v4.new_notebook(metadata=self.nb_metadata)
@@ -242,7 +242,7 @@ def test_exercise2_conversion(self):
             nb['cells'].append(cell_md)
             nbformat.write(nb, f)
         # run command and check for errors
-        cmd = ['exercise2', '--to-py', f_input]
+        cmd = ['exercise2', '--to-py', str(f_input)]
         self.run_command(cmd, f_input)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f:
@@ -259,7 +259,7 @@ def test_exercise2_conversion(self):
         self.assertEqual(cell['metadata']['key'], 'value')
         self.assertEqual(next(cells, 'EOF'), 'EOF')
         # run command and check for errors
-        cmd = ['exercise2', '--to-md', f_input]
+        cmd = ['exercise2', '--to-md', str(f_input)]
         self.run_command(cmd, f_input)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f:
@@ -278,7 +278,8 @@ def test_exercise2_conversion(self):
 
     @skipIfMissingModules
     def test_exercise2_autopep8(self):
-        f_input = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_exercise2_autopep8.ipynb'
+        root = pathlib.Path("@CMAKE_CURRENT_BINARY_DIR@")
+        f_input = root / "test_convert_exercise2_autopep8.ipynb"
         # setup
         with open(f_input, 'w', encoding='utf-8') as f:
             nb = nbformat.v4.new_notebook(metadata=self.nb_metadata)
@@ -293,7 +294,7 @@ def test_exercise2_autopep8(self):
             nb['cells'].append(cell_md)
             nbformat.write(nb, f)
         # run command and check for errors
-        cmd = ['exercise2', '--pep8', f_input]
+        cmd = ['exercise2', '--pep8', str(f_input)]
         self.run_command(cmd)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f: