diff --git a/.github/actions/dependencies/action.yml b/.github/actions/dependencies/action.yml
index 03953a1..95f3372 100644
--- a/.github/actions/dependencies/action.yml
+++ b/.github/actions/dependencies/action.yml
@@ -11,6 +11,7 @@ runs:
   using: "composite"
   steps:
     - run: |
+        test "${RUNNER_ARCH}" = "X64" && module use /cvmfs/dev.eessi.io/espresso/versions/${EESSI_VERSION}/software/linux/x86_64/amd/zen2/modules/all
         module load ${{ inputs.modules }}
         module save pymbe
         python3 -m venv --system-site-packages venv
diff --git a/.github/workflows/samples.yml b/.github/workflows/samples.yml
index b24527f..ab79356 100644
--- a/.github/workflows/samples.yml
+++ b/.github/workflows/samples.yml
@@ -17,28 +17,29 @@ jobs:
       OMPI_MCA_mtl_ofi_provider_exclude: psm3
     strategy:
       matrix:
-        espresso:
-          - version: "4.2.2"
-            eessi_modules: ESPResSo/4.2.2-foss-2023b
-            eessi_stack_version: "2023.06"
-            upload_artifact: true
-          - version: "5.0.1"
+        software:
+          - label: "ESPResSo 5.0.1"
             eessi_modules: ESPResSo/5.0.1-foss-2025a
             eessi_stack_version: "2025.06"
+            upload_artifact: true
+          - label: "ESPResSo 5.1-dev"
+            eessi_modules: ESPResSo/cd7547c43b2dae8e96487bf149d049d701dddd43-foss-2025a
+            eessi_stack_version: "2025.06"
             upload_artifact: false
-    name: ubuntu - ESPResSo ${{ matrix.espresso.version }}
+    name: ubuntu - ${{ matrix.software.label }}
     steps:
+      - name: Checkout repository
+        uses: actions/checkout@main
       - name: Setup EESSI
         uses: eessi/github-action-eessi@v3
         with:
-          eessi_stack_version: ${{ matrix.espresso.eessi_stack_version }}
-      - name: Checkout repository
-        uses: actions/checkout@main
+          eessi_stack_version: ${{ matrix.software.eessi_stack_version }}
+          use_eessi_module: false
       - name: Install dependencies
         uses: ./.github/actions/dependencies
         with:
           modules: |-
-            ${{ matrix.espresso.eessi_modules }}
+            ${{ matrix.software.eessi_modules }}
       - name: Run testsuite
         run: |
           export NUM_PROC=$(nproc)
diff --git a/.github/workflows/testsuite.yml b/.github/workflows/testsuite.yml
index ecb2ef9..04b29f3 100644
--- a/.github/workflows/testsuite.yml
+++ b/.github/workflows/testsuite.yml
@@ -17,13 +17,13 @@ jobs:
     strategy:
       matrix:
         software:
-          - label: "ESPResSo 4.2.2, LAMMPS 2024"
-            eessi_modules: ESPResSo/4.2.2-foss-2023b LAMMPS/29Aug2024-foss-2023b-kokkos
-            eessi_stack_version: "2023.06"
-            upload_artifact: true
           - label: "ESPResSo 5.0.1"
             eessi_modules: ESPResSo/5.0.1-foss-2025a
             eessi_stack_version: "2025.06"
+            upload_artifact: true
+          - label: "ESPResSo 5.1-dev"
+            eessi_modules: ESPResSo/cd7547c43b2dae8e96487bf149d049d701dddd43-foss-2025a
+            eessi_stack_version: "2025.06"
             upload_artifact: false
     name: ubuntu - ${{ matrix.software.label }}
     steps:
@@ -37,6 +37,7 @@ jobs:
         uses: eessi/github-action-eessi@v3
         with:
           eessi_stack_version: ${{ matrix.software.eessi_stack_version }}
+          use_eessi_module: false
       - name: Install dependencies
         uses: ./.github/actions/dependencies
         with:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 05025f7..a15004b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -45,6 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Methods that interact directly with the pyMBE dataframe. These methods have been replaced by private methods that instead interact with the new canonical pyMBE database in (`pyMBE/storage/manager`). This includes the methods: `add_bond_in_df`, `add_value_to_df`, `assign_molecule_id`, `check_if_df_cell_has_a_value`, `check_if_name_is_defined_in_df`, `check_if_multiple_pmb_types_for_name`, `clean_df_row`, `clean_ids_in_df_row`, `copy_df_entry`, `create_variable_with_units`, `convert_columns_to_original_format`, `convert_str_to_bond_object`, `delete_entries_in_df`, `find_bond_key`, `setup_df`, `define_particle_entry_in_df`, custom `NumpyEncoder`. (#145,#147)
 - Method `add_bonds_to_espresso` has been removed from the API. pyMBE now adds bonds internally to ESPResSo when molecule instances are created into ESPResSo. (#147)
 - Tutorial `lattice_builder.ipynb` has been removed because its content is redundant with sample script `build_hydrogel.py`. (#147)
+- Legacy ESPResSo 4.2-specific helper functions `do_reaction` and `get_number_of_particles` were removed from the API. ESPResSo 4.2 is no longer officially supported. (#155)
 
 
 ## [1.0.0] - 2025-10-08
diff --git a/README.md b/README.md
index 0212541..e4bac0f 100644
--- a/README.md
+++ b/README.md
@@ -68,8 +68,8 @@ git clone git@github.com:pyMBE-dev/pyMBE.git
 
 Please, be aware that pyMBE is intended to be a supporting tool to setup simulations with ESPResSo.
 Thus, for most of its functionalities ESPResSo must also be available.
-pyMBE supports ESPResSo 4.2 and ESPResSo 4.3-dev.
-Following the NEP29 guidelines, we recommend using Python3.10+.
+pyMBE supports ESPResSo 5.0 and ESPResSo 5.1-dev.
+Following the NEP29 guidelines, we recommend using Python3.11+.
 Both NumPy 1 and NumPy 2 are supported.
 
 The pyMBE module needs a Python virtual environment to avoid compatibility issues with its dependencies.
@@ -99,7 +99,7 @@ python3 -m pip install -r requirements.txt "numpy>=2.1" "pandas>=2.0"
 We highlight that the path `/home/user/espresso/build` is just an example of a possible path to the ESPResSo build folder. 
 The user should change this path to match the local absolute path where ESPResSo was built.
 Also, ESPResSo must be built with the same NumPy version as the one installed in the environment to avoid API version mismatch.
-For more details on how to install ESPResSo, please consult the [ESPResSo installation guide](https://espressomd.github.io/doc4.2.2/installation.html).
+For more details on how to install ESPResSo, please consult the [ESPResSo installation guide](https://espressomd.github.io/doc5.0.1/installation.html).
 
 The pyMBE virtual environment can be deactivated at any moment as follows:
 
@@ -111,7 +111,7 @@ Cluster users who rely on module files to load dependencies should opt for the
 following alternative:
 
 ```sh
-module load ESPResSo/4.2.2-foss-2023a # adapt release if needed
+module load ESPResSo/5.0.1-foss-2025a # adapt release if needed
 python3 -m venv --system-site-packages pymbe
 source pymbe/bin/activate
 python3 maintainer/configure_venv.py
@@ -128,7 +128,7 @@ Now you can use pyMBE and ESPResSo by activating the virtual environment:
 ```sh
 $ source pymbe/bin/activate
 (pymbe) $ python3 -c "import espressomd.version; print(espressomd.version.friendly())"
-4.2
+5.0.1
 (pymbe) $ python3 -c "import pyMBE; print(pyMBE.__file__)"
 /home/user/Documents/pyMBE/pyMBE/__init__.py
 $ deactivate
diff --git a/pyMBE/simulation_builder/espresso_engine.py b/pyMBE/simulation_builder/espresso_engine.py
index 1cbf013..68fe830 100644
--- a/pyMBE/simulation_builder/espresso_engine.py
+++ b/pyMBE/simulation_builder/espresso_engine.py
@@ -37,7 +37,6 @@ def __init__(self,box_l,db,espresso_system,units,kT,Kw,seed):
         self.kT=kT
         self.Kw=Kw
         self.seed=seed
-        pass
 
     def _add_angle(self,particle_id1,particle_id2,particle_id3, angle_inst):
         """ 
@@ -365,36 +364,6 @@ def change_volume_and_rescale_particles(self, d_new, dir="xyz"):
             
         self.espresso_system.change_volume_and_rescale_particles(d_new=d_new,
                                                                  dir=dir)
-    
-    
-    def do_reaction(self,algorithm, steps):
-        """
-        Executes reaction steps using an ESPResSo reaction algorithm with
-        version-compatible calling semantics.
-
-        This function wraps the `reaction` method of an ESPResSo reaction
-        algorithm to account for differences in the method signature between
-        ESPResSo versions.
-
-        Args:
-            algorithm ('espressomd.reaction_methods'):
-                ESPResSo reaction algorithm object (e.g. constant pH,
-                reaction ensemble, or similar).
-            steps ('int'):
-                Number of reaction steps to perform.
-
-        Notes:
-            - In ESPResSo 4.2, the `reaction` method expects the number of steps
-            to be passed as the keyword argument `reaction_steps`.
-            - In newer ESPResSo versions, the keyword argument is `steps`.
-            - This helper function provides a stable interface across ESPResSo
-            versions by dispatching to the appropriate keyword internally.
-        """
-        import espressomd.version
-        if espressomd.version.friendly() == '4.2':
-            algorithm.reaction(reaction_steps=steps)
-        else:
-            algorithm.reaction(steps=steps)
 
     def enable_motion_of_rigid_object(self, instance_id, pmb_type):
         """
@@ -442,35 +411,6 @@ def enable_motion_of_rigid_object(self, instance_id, pmb_type):
             pid = self.espresso_system.part.by_id(particle_id)
             pid.vs_auto_relate_to(rigid_object_center.id)
 
-    def get_number_of_particles(self, ptype):
-        """
-        Returns the number of particles of a given ESPResSo particle type.
-
-        Args:
-            ptype ('int'):
-                ESPResSo particle type identifier.
-
-        Returns:
-            ('int'):
-                Number of particles in `espresso_system` with particle type `ptype`.
-
-        Notes:
-            - In ESPResSo 4.2, `number_of_particles` expects the particle type
-            as a positional argument.
-            - In later ESPResSo versions, the particle type must be passed as a
-            keyword argument (`type=ptype`).
-            - This helper function hides these API differences and provides
-            a uniform interface across ESPResSo versions.
-        """
-        import espressomd.version
-        if espressomd.version.friendly() == "4.2":
-            args = (ptype,)
-            kwargs = {}
-        else:
-            args = ()
-            kwargs = {"type": ptype}
-        return self.espresso_system.number_of_particles(*args, **kwargs)
-    
     def relax_espresso_system(self, seed, gamma=1e-3, Nsteps_steepest_descent=5000, max_displacement=0.01, Nsteps_iter_relax=500):
         """
         Relaxes the energy of the given ESPResSo system by performing the following steps:
@@ -610,19 +550,13 @@ def setup_electrostatic_interactions(self,units, kT, c_salt=None, solvent_permit
 
             if tune_p3m:
                 self.espresso_system.time_step=0.01
-                if espressomd.version.friendly() == "4.2":
-                    self.espresso_system.actors.add(coulomb)
-                else:
-                    self.espresso_system.electrostatics.solver = coulomb
+                self.espresso_system.electrostatics.solver = coulomb
 
 
                 # save the optimal parameters and add them by hand
 
                 p3m_params = coulomb.get_params()
-                if espressomd.version.friendly() == "4.2":
-                    self.espresso_system.actors.remove(coulomb)
-                else:
-                    self.espresso_system.electrostatics.solver = None
+                self.espresso_system.electrostatics.solver = None
                 coulomb = espressomd.electrostatics.P3M(prefactor = COULOMB_PREFACTOR.m_as("reduced_length * reduced_energy"),
                                                         accuracy = accuracy,
                                                         mesh = p3m_params['mesh'],
@@ -641,10 +575,7 @@ def setup_electrostatic_interactions(self,units, kT, c_salt=None, solvent_permit
             coulomb = espressomd.electrostatics.DH(prefactor = COULOMB_PREFACTOR.m_as("reduced_length * reduced_energy"), 
                                                 kappa = (1./KAPPA).to('1/ reduced_length').magnitude, 
                                                 r_cut = r_cut)
-        if espressomd.version.friendly() == "4.2":
-            self.espresso_system.actors.add(coulomb)
-        else:
-            self.espresso_system.electrostatics.solver = coulomb
+        self.espresso_system.electrostatics.solver = coulomb
         logging.debug("*** Electrostatics successfully added to the system ***")
 
     def setup_cpH (self, counter_ion, constant_pH, exclusion_range=None, use_exclusion_radius_per_type = False):
@@ -676,11 +607,15 @@ def setup_cpH (self, counter_ion, constant_pH, exclusion_range=None, use_exclusi
             exclusion_radius_per_type = self.db.get_radius_map()
         else:
             exclusion_radius_per_type = {}
+        kwargs = {}
+        if espressomd.version.version() >= (5, 1, 0):
+            kwargs["system"] = self.espresso_system
         RE = reaction_methods.ConstantpHEnsemble(kT=self.kT.to('reduced_energy').magnitude,
                                                 exclusion_range=exclusion_range, 
                                                 seed=self.seed, 
                                                 constant_pH=constant_pH,
-                                                exclusion_radius_per_type = exclusion_radius_per_type)
+                                                exclusion_radius_per_type = exclusion_radius_per_type,
+                                                **kwargs)
         conterion_tpl = self.db.get_template(name=counter_ion,
                                              pmb_type="particle")
         conterion_state = self.db.get_template(name=conterion_tpl.initial_state,
@@ -749,10 +684,14 @@ def setup_gcmc(self, c_salt_res, salt_cation_name, salt_anion_name, activity_coe
             exclusion_radius_per_type = self.db.get_radius_map()
         else:
             exclusion_radius_per_type = {}
+        kwargs = {}
+        if espressomd.version.version() >= (5, 1, 0):
+            kwargs["system"] = self.espresso_system
         RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
                                                exclusion_range=exclusion_range, 
                                                seed=self.seed, 
-                                               exclusion_radius_per_type = exclusion_radius_per_type)
+                                               exclusion_radius_per_type = exclusion_radius_per_type,
+                                               **kwargs)
         # Determine the concentrations of the various species in the reservoir and the equilibrium constants
         determined_activity_coefficient = activity_coefficient(c_salt_res)
         K_salt = (c_salt_res.to('1/(N_A * reduced_length**3)')**2) * determined_activity_coefficient
@@ -846,10 +785,14 @@ def setup_grxmc_reactions(self, pH_res, c_salt_res, proton_name, hydroxide_name,
             exclusion_radius_per_type = self.db.get_radius_map()
         else:
             exclusion_radius_per_type = {}
+        kwargs = {}
+        if espressomd.version.version() >= (5, 1, 0):
+            kwargs["system"] = self.espresso_system
         RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
                                                exclusion_range=exclusion_range, 
                                                seed=self.seed, 
-                                               exclusion_radius_per_type = exclusion_radius_per_type)
+                                               exclusion_radius_per_type = exclusion_radius_per_type,
+                                               **kwargs)
         # Determine the concentrations of the various species in the reservoir and the equilibrium constants
         cH_res, cOH_res, cNa_res, cCl_res = self.determine_reservoir_concentrations(pH_res, c_salt_res, activity_coefficient)
         ionic_strength_res = 0.5*(cNa_res+cCl_res+cOH_res+cH_res)
@@ -1142,10 +1085,14 @@ def setup_grxmc_unified(self, pH_res, c_salt_res, cation_name, anion_name, activ
             exclusion_radius_per_type = self.db.get_radius_map()
         else:
             exclusion_radius_per_type = {}
+        kwargs = {}
+        if espressomd.version.version() >= (5, 1, 0):
+            kwargs["system"] = self.espresso_system
         RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
                                                exclusion_range=exclusion_range, 
                                                seed=self.seed, 
-                                               exclusion_radius_per_type = exclusion_radius_per_type)
+                                               exclusion_radius_per_type = exclusion_radius_per_type,
+                                               **kwargs)
         # Determine the concentrations of the various species in the reservoir and the equilibrium constants
         cH_res, cOH_res, cNa_res, cCl_res = self.determine_reservoir_concentrations(pH_res, c_salt_res, activity_coefficient)
         ionic_strength_res = 0.5*(cNa_res+cCl_res+cOH_res+cH_res)
@@ -1319,7 +1266,7 @@ def setup_lj_interactions(self, shift_potential=True, combining_rule='Lorentz-Be
 
         Notes:
             - Currently, the only 'combining_rule' supported is Lorentz-Berthelot.
-            - Check the documentation of ESPResSo for more info about the potential https://espressomd.github.io/doc4.2.0/inter_non-bonded.html
+            - Check the documentation of ESPResSo for more info about the potential https://espressomd.github.io/doc5.0.1/inter_non-bonded.html
 
         """
         from itertools import combinations_with_replacement
diff --git a/samples/Beyer2024/globular_protein.py b/samples/Beyer2024/globular_protein.py
index ea07025..eeec6db 100644
--- a/samples/Beyer2024/globular_protein.py
+++ b/samples/Beyer2024/globular_protein.py
@@ -19,6 +19,7 @@
 from pathlib import Path
 import tqdm
 import espressomd
+import espressomd.version
 import argparse
 import numpy as np
 import pandas as pd 
@@ -268,8 +269,8 @@
     print(pmb.get_reactions_df())
 
 type_map = pmb.get_type_map()
-types = list (type_map.values())
-espresso_system.setup_type_map( type_list = types)
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = type_map.values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
 non_interacting_type = max(type_map.values())+1
@@ -320,7 +321,7 @@
 
 for step in tqdm.trange(N_samples, disable=not verbose):
     espresso_system.integrator.run (steps = integ_steps)
-    pmb.simulation_engine.do_reaction(cpH, steps=total_ionisable_groups)
+    cpH.reaction(steps=total_ionisable_groups)
     protein_net_charge = pmb.calculate_net_charge(
                                                 object_name=protein_name,
                                                 pmb_type="protein",
diff --git a/samples/Beyer2024/peptide.py b/samples/Beyer2024/peptide.py
index 8fdda11..ca8009f 100644
--- a/samples/Beyer2024/peptide.py
+++ b/samples/Beyer2024/peptide.py
@@ -19,6 +19,7 @@
 # Load espresso, pyMBE and other necessary libraries
 from pathlib import Path
 import espressomd
+import espressomd.version
 import pandas as pd
 import argparse
 import tqdm
@@ -184,11 +185,11 @@
     print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
-type_map =pmb.get_type_map()
-espresso_system.setup_type_map(type_list = list(type_map.values()))
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
-non_interacting_type = max(type_map.values())+1
+non_interacting_type = max(pmb.get_type_map().values())+1
 cpH.set_non_interacting_type (type=non_interacting_type)
 if verbose:
     print(f"The non-interacting type is set to {non_interacting_type}")
@@ -224,7 +225,7 @@
     # Run LD
     espresso_system.integrator.run(steps=MD_steps_per_sample)
     # Run MC
-    pmb.simulation_engine.do_reaction(cpH, steps=len(sequence))
+    cpH.reaction(steps=len(sequence))
     # Sample observables
     charge_dict=pmb.calculate_net_charge(
                                         object_name=sequence,
diff --git a/samples/Beyer2024/weak_polyelectrolyte_dialysis.py b/samples/Beyer2024/weak_polyelectrolyte_dialysis.py
index aa24618..c1e3fd0 100644
--- a/samples/Beyer2024/weak_polyelectrolyte_dialysis.py
+++ b/samples/Beyer2024/weak_polyelectrolyte_dialysis.py
@@ -16,12 +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/>.
 
-#######################################################
-# Loading modules 
-#######################################################
-
 # Load python modules
 import espressomd
+import espressomd.version
 from pathlib import Path
 import numpy as np
 import pandas as pd
@@ -36,8 +33,6 @@
 # Create an instance of pyMBE library
 pmb = pyMBE.pymbe_library(seed=42)
 
-# Load some functions from the handy_scripts library for convenience
-
 
 #######################################################
 # Setting parameters for the simulation
@@ -209,13 +204,12 @@
     print("The acid-base reaction has been successfully set up for:")
     print(pmb.get_reactions_df())
 
-# Setup espresso to track the ionization of the acid groups
-type_map = pmb.get_type_map()
-types = list(type_map.values())
-espresso_system.setup_type_map(type_list = types)
+# Setup espresso to track the ionization of the acid/basic groups in peptide
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
-non_interacting_type = max(type_map.values())+1
+non_interacting_type = max(pmb.get_type_map().values())+1
 grxmc.set_non_interacting_type (type=non_interacting_type)
 
 #Set up the interactions
@@ -233,7 +227,7 @@
     print("Running warmup without electrostatics")
 for i in tqdm.trange(100, disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc, steps=1000)
+    grxmc.reaction(steps=1000)
 
 pmb.simulation_engine.setup_electrostatic_interactions(units=pmb.units,
                                 kT=pmb.kT,
@@ -257,7 +251,7 @@
     N_warmup_loops = 100
 for i in tqdm.trange(N_warmup_loops, disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc, steps=100)
+    grxmc.reaction(steps=100)
 
 # Main loop
 print("Started production run.")
@@ -273,7 +267,7 @@
     N_production_loops = 100
 for i in tqdm.trange(N_production_loops, disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc, steps=100)
+    grxmc.reaction(steps=100)
     # Measure time
     time_series["time"].append(espresso_system.time)
     # Measure degree of ionization
diff --git a/samples/branched_polyampholyte.py b/samples/branched_polyampholyte.py
index ed9b0d9..29353a8 100644
--- a/samples/branched_polyampholyte.py
+++ b/samples/branched_polyampholyte.py
@@ -16,16 +16,14 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-# Load espresso, pyMBE and other necessary libraries
 from pathlib import Path
 import espressomd
+import espressomd.version
 import argparse
 import tqdm
 import pandas as pd
 from espressomd.io.writer import vtf
 import pyMBE
-
-# Load some functions from the handy_scripts library for convenience
 from pyMBE.lib.analysis import built_output_name
 
 
@@ -181,12 +179,11 @@
     print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups 
-type_map = pmb.get_type_map()
-types = list(type_map.values())
-espresso_system.setup_type_map(type_list = types)
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
-non_interacting_type = max(type_map.values())+1
+non_interacting_type = max(pmb.get_type_map().values())+1
 cpH.set_non_interacting_type (type=non_interacting_type)
 if verbose:
     print(f"The non interacting type is set to {non_interacting_type}")
@@ -232,9 +229,9 @@
 
 # Production loop
 N_frame=0
-for step in tqdm.trange(N_samples):
+for step in tqdm.trange(N_samples, disable=args.test):
     espresso_system.integrator.run(steps=MD_steps_per_sample)        
-    pmb.simulation_engine.do_reaction(cpH, steps=total_ionisable_groups)   
+    cpH.reaction(steps=total_ionisable_groups)
     # Get polyampholyte net charge
     charge_dict=pmb.calculate_net_charge(
                                         object_name="polyampholyte",
diff --git a/samples/peptide_cpH.py b/samples/peptide_cpH.py
index 78076a5..c62edee 100644
--- a/samples/peptide_cpH.py
+++ b/samples/peptide_cpH.py
@@ -18,6 +18,7 @@
 
 # Load espresso, pyMBE and other necessary libraries
 import espressomd
+import espressomd.version
 import pandas as pd
 import tqdm
 from espressomd.io.writer import vtf
@@ -72,7 +73,7 @@
 if args.test:
     MD_steps_per_sample = 1
     ideal=True
-    N_samples = 2000 # improve sampling for testing
+    N_samples = 3200 # improve sampling for testing
 
 # Peptide parameters
 sequence = args.sequence
@@ -191,12 +192,11 @@
     print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
-type_map =pmb.get_type_map()
-types = list (type_map.values())
-espresso_system.setup_type_map( type_list = types)
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
-non_interacting_type = max(type_map.values())+1
+non_interacting_type = max(pmb.get_type_map().values())+1
 cpH.set_non_interacting_type (type=non_interacting_type)
 if verbose:
     print('The non-interacting type is set to ', non_interacting_type)
@@ -240,11 +240,11 @@
 
 # Main loop for performing simulations at different pH-values
 N_frame=0
-for sample in tqdm.trange(N_samples):
+for sample in tqdm.trange(N_samples, disable=args.test):
     # LD sampling of the configuration space
     espresso_system.integrator.run(steps=MD_steps_per_sample)        
     # cpH sampling of the reaction space
-    pmb.simulation_engine.do_reaction(cpH, steps=total_ionisable_groups) # rule of thumb: one reaction step per titratable group (on average)
+    cpH.reaction(steps=total_ionisable_groups) # rule of thumb: one reaction step per titratable group (on average)
     # Get peptide net charge
     charge_dict=pmb.calculate_net_charge(
                                         object_name=peptide_name,
diff --git a/samples/peptide_mixture_grxmc_ideal.py b/samples/peptide_mixture_grxmc_ideal.py
index 3cb6a48..d26cc13 100644
--- a/samples/peptide_mixture_grxmc_ideal.py
+++ b/samples/peptide_mixture_grxmc_ideal.py
@@ -16,8 +16,8 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#Load espresso, pyMBE and other necessary libraries
 import espressomd
+import espressomd.version
 from pathlib import Path
 import pandas as pd
 import argparse
@@ -269,11 +269,11 @@
     print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
-type_map =pmb.get_type_map()
-types = list (type_map.values())
-espresso_system.setup_type_map(type_list = types)
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
+type_map = pmb.get_type_map()
 non_interacting_type = max(type_map.values())+1
 grxmc.set_non_interacting_type (type=non_interacting_type)
 if verbose:
@@ -302,7 +302,7 @@
 N_frame=0
 for step in range(N_samples):
     espresso_system.integrator.run(steps=MD_steps_per_sample)        
-    pmb.simulation_engine.do_reaction(grxmc, steps=total_ionisable_groups)
+    grxmc.reaction(steps=total_ionisable_groups)
     time_series["time"].append(espresso_system.time)
     # Get net charge of peptide1 and peptide2
     charge_dict_peptide1=pmb.calculate_net_charge(
diff --git a/samples/salt_solution_gcmc.py b/samples/salt_solution_gcmc.py
index b3900ab..5f6dc7c 100644
--- a/samples/salt_solution_gcmc.py
+++ b/samples/salt_solution_gcmc.py
@@ -19,6 +19,7 @@
 # Load python modules
 from pathlib import Path
 import espressomd
+import espressomd.version
 import numpy as np
 import pandas as pd
 from scipy import interpolate
@@ -28,8 +29,6 @@
 # Import pyMBE
 import pyMBE
 from pyMBE.lib import analysis
-#Import functions from handy_functions script 
-
 
 # Create an instance of pyMBE library
 pmb = pyMBE.pymbe_library(seed=42)
@@ -129,11 +128,11 @@
     print("Set up GCMC...")
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
-type_map = pmb.get_type_map()
-types = list (type_map.values())
-espresso_system.setup_type_map(type_list = types)
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
+type_map = pmb.get_type_map()
 non_interacting_type = max(type_map.values())+1
 RE.set_non_interacting_type(type=non_interacting_type)
 if verbose:
@@ -159,7 +158,7 @@
     print("Running warmup without electrostatics")
 for i in tqdm.trange(100, disable=not verbose):
     espresso_system.integrator.run(steps=100)
-    pmb.simulation_engine.do_reaction(RE, steps=100)
+    RE.reaction(steps=100)
 
 if args.mode == "interacting":
     pmb.simulation_engine.setup_electrostatic_interactions(units=pmb.units,
@@ -181,7 +180,7 @@
 N_warmup_loops = 100
 for i in tqdm.trange(N_warmup_loops, disable=not verbose):
     espresso_system.integrator.run(steps=100)
-    pmb.simulation_engine.do_reaction(RE, steps=100)
+    RE.reaction(steps=100)
 
 # Main loop
 print("Started production run.")
@@ -195,13 +194,13 @@
 N_production_loops = 100
 for i in tqdm.trange(N_production_loops, disable=not verbose):
     espresso_system.integrator.run(steps=100)
-    pmb.simulation_engine.do_reaction(RE, steps=100)
+    RE.reaction(steps=100)
 
     # Measure time
     time_series["time"].append(espresso_system.time)
 
     # Measure degree of ionization
-    number_of_ion_pairs = pmb.simulation_engine.get_number_of_particles(type_map[cation_name])
+    number_of_ion_pairs = espresso_system.number_of_particles(type=type_map[cation_name])
     time_series["c_salt"].append((number_of_ion_pairs/(volume * pmb.N_A)).magnitude)
 
 data_path = args.output
diff --git a/samples/weak_polyacid_hydrogel_grxmc.py b/samples/weak_polyacid_hydrogel_grxmc.py
index 944118e..270e87d 100644
--- a/samples/weak_polyacid_hydrogel_grxmc.py
+++ b/samples/weak_polyacid_hydrogel_grxmc.py
@@ -18,6 +18,7 @@
 #
 
 import espressomd
+import espressomd.version
 from espressomd.io.writer import vtf
 from pathlib import Path
 import numpy as np
@@ -235,18 +236,17 @@
                                                     salt_anion_name=chloride_name, 
                                                     activity_coefficient=activity_coefficient_monovalent_pair)
 
-# Setup espresso to track the ionization of the acid groups
-type_map = pmb.get_type_map()
-types = list(type_map.values())
-espresso_system.setup_type_map(type_list = types)
+# Setup espresso to track the ionization of the acid/basic groups in peptide
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
-non_interacting_type = max(type_map.values())+1
+non_interacting_type = max(pmb.get_type_map().values())+1
 grxmc.set_non_interacting_type (type=non_interacting_type)
 
 for i in tqdm.trange(100):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc,1000)
+    grxmc.reaction(steps=1000)
 
 pmb.simulation_engine.setup_electrostatic_interactions(units=pmb.units,
                                  kT=pmb.kT,
@@ -260,7 +260,7 @@
 print("*** Running warmup with electrostatics... ***")
 for i in tqdm.trange(N_warmup_loops):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc,100)
+    grxmc.reaction(steps=100)
 
 # Main loop
 print("*** Starting production run... ***")
@@ -280,7 +280,7 @@
 
 for i in tqdm.trange(N_production_loops):
     espresso_system.integrator.run(steps=1000)
-    pmb.simulation_engine.do_reaction(grxmc,200)
+    grxmc.reaction(steps=200)
 
     # Measure time
     time_series["time"].append(espresso_system.time)
diff --git a/testsuite/setup_salt_ions_unit_tests.py b/testsuite/setup_salt_ions_unit_tests.py
index 67caad9..f3c83f5 100644
--- a/testsuite/setup_salt_ions_unit_tests.py
+++ b/testsuite/setup_salt_ions_unit_tests.py
@@ -44,7 +44,6 @@
                     sigma=0.3*pmb.units.nm,
                     epsilon=1*pmb.units.Quantity(1,"reduced_energy"))
 
-type_map=pmb.get_type_map()
 # System parameters
 c_salt_input = 0.01 * pmb.units.mol/ pmb.units.L
 N_SALT_ION_PAIRS = 50
@@ -53,7 +52,8 @@
 box_l=[L.to('reduced_length').magnitude]*3
 # Create an instance of an espresso system
 espresso_system=espressomd.System (box_l = box_l )
-espresso_system.setup_type_map(type_list=type_map.values())
+if espressomd.version.version() < (5, 1, 0):
+    espresso_system.setup_type_map(type_list = pmb.get_type_map().values())
 
 pmb.define_particle(name='0P',
                     z=0,
@@ -124,14 +124,13 @@ def test_salt_addition(self):
         def check_salt_concentration(espresso_system,cation_name,anion_name,c_salt,N_SALT_ION_PAIRS):
             charge_number_map=pmb.get_charge_number_map()
             type_map=pmb.get_type_map()
-            espresso_system.setup_type_map(type_list=type_map.values())
             c_salt_calculated = pmb.create_added_salt(box_l=box_l,
                                                     cation_name=cation_name,
                                                     anion_name=anion_name,
                                                     c_salt=c_salt)
             pmb.add_instances_to_engine()
-            self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map[cation_name]),N_SALT_ION_PAIRS*abs(charge_number_map[type_map[anion_name]]))
-            self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map[anion_name]),N_SALT_ION_PAIRS*abs(charge_number_map[type_map[cation_name]]))
+            self.assertEqual(espresso_system.number_of_particles(type=type_map[cation_name]),N_SALT_ION_PAIRS*abs(charge_number_map[type_map[anion_name]]))
+            self.assertEqual(espresso_system.number_of_particles(type=type_map[anion_name]),N_SALT_ION_PAIRS*abs(charge_number_map[type_map[cation_name]]))
             self.assertAlmostEqual(c_salt_calculated.m_as("mol/L"), c_salt.m_as("mol/L"))
             cation_ids = pmb.get_particle_id_map(object_name=cation_name)["all"]
             anion_ids  = pmb.get_particle_id_map(object_name=anion_name)["all"]
@@ -168,15 +167,15 @@ def test_salt_addition_concentration_units(self):
         """
         Unit test: check that create_added_salt works for an input c_salt in [particle/lenght**3]. 
         """
+        type_map=pmb.get_type_map()
         c_salt_part=c_salt_input*pmb.N_A
-        espresso_system.setup_type_map(type_list=type_map.values())
         c_salt_calculated = pmb.create_added_salt(box_l=box_l,
                                                     cation_name="Na",
                                                     anion_name="Cl",
                                                     c_salt=c_salt_part)
         pmb.add_instances_to_engine()
-        self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map["Na"]),N_SALT_ION_PAIRS)
-        self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map["Cl"]),N_SALT_ION_PAIRS)
+        self.assertEqual(espresso_system.number_of_particles(type=type_map["Na"]),N_SALT_ION_PAIRS)
+        self.assertEqual(espresso_system.number_of_particles(type=type_map["Cl"]),N_SALT_ION_PAIRS)
         self.assertAlmostEqual(c_salt_calculated.m_as("reduced_length**-3"), c_salt_part.m_as("reduced_length**-3"))
         cation_ids = pmb.get_particle_id_map(object_name="Na")["all"]
         anion_ids = pmb.get_particle_id_map(object_name="Cl")["all"]
@@ -233,13 +232,13 @@ def test_counterions(molecule_name, cation_name, anion_name, espresso_system, ex
                                     anion_name=anion_name,
                                     box_l=box_l)
             pmb.add_instances_to_engine()
-            espresso_system.setup_type_map(type_list=type_map.values())
-            self.assertEqual(pmb.simulation_engine.get_number_of_particles(
-                                                     type_map[cation_name]),
-                                                     expected_numbers[cation_name])
-            self.assertEqual(pmb.simulation_engine.get_number_of_particles(
-                                                     type_map[anion_name]),
-                                                     expected_numbers[anion_name])
+            type_map=pmb.get_type_map()
+            self.assertEqual(
+                espresso_system.number_of_particles(type=type_map[cation_name]),
+                expected_numbers[cation_name])
+            self.assertEqual(
+                espresso_system.number_of_particles(type=type_map[anion_name]),
+                expected_numbers[anion_name])
             molecule_ids = list(pmb.get_particle_id_map(object_name=molecule_name)["molecule_map"].keys())
             for mol_id in molecule_ids:
                 pmb.delete_instances_in_system(instance_id=mol_id,
@@ -289,6 +288,7 @@ def test_sanity_create_counterions(self):
                                 box_l=box_l,
                                 use_default_bond=True)
         pmb.add_instances_to_engine()
+        type_map=pmb.get_type_map()
         input_parameters={"cation_name":"Ca",
                             "anion_name":"Cl",
                             "object_name":'isoelectric_polyampholyte',
@@ -312,10 +312,9 @@ def test_sanity_create_counterions(self):
                                     anion_name="Cl",
                                     box_l=box_l)
         pmb.add_instances_to_engine()
-        espresso_system.setup_type_map(type_list=type_map.values())
 
-        self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map["Na"]),0)
-        self.assertEqual(pmb.simulation_engine.get_number_of_particles(type_map["Cl"]),0)
+        self.assertEqual(espresso_system.number_of_particles(type=type_map["Na"]),0)
+        self.assertEqual(espresso_system.number_of_particles(type=type_map["Cl"]),0)
         # Assert that no counterions are created if the wrong object names are provided
         inputs = {"object_name":'test',
                 "cation_name":"Na",
diff --git a/testsuite/test_handy_functions.py b/testsuite/test_handy_functions.py
index 6939e25..f5ec153 100644
--- a/testsuite/test_handy_functions.py
+++ b/testsuite/test_handy_functions.py
@@ -104,12 +104,9 @@ def test_langevin_setup(self):
                          second=espresso_system.time_step,
                          msg="The input time step in `lib.handy_functions.setup_langevin_dynamics` is not consistent with the one in the espresso simulation System")
         ## Test setup of the thermostat
-        if espressomd.version.friendly() == "4.2":
-            thermostat_setup=espresso_system.thermostat.get_state()[0]
-        else:
-            thermostat_setup=espresso_system.thermostat.langevin.get_params()
-            thermostat_setup["kT"] = espresso_system.thermostat.kT
-            thermostat_setup["type"] = "LANGEVIN"
+        thermostat_setup=espresso_system.thermostat.langevin.get_params()
+        thermostat_setup["kT"] = espresso_system.thermostat.kT
+        thermostat_setup["type"] = "LANGEVIN"
         self.assertEqual(first="LANGEVIN",
                          second=thermostat_setup["type"],
                          msg="`lib.handy_functions.setup_langevin_dynamics` is setting a different thermostat than Langevin")
@@ -210,10 +207,7 @@ def test_setup_electrostatics(self):
         coulomb_prefactor=Bjerrum_length*electrostatics_inputs["kT"]
         # Test the P3M setup
         pmb.simulation_engine.setup_electrostatic_interactions(**electrostatics_inputs)
-        if espressomd.version.friendly() == "4.2":
-            coulomb = espresso_system.actors.active_actors.copy()[0]
-        else:
-            coulomb = espresso_system.electrostatics.solver
+        coulomb = espresso_system.electrostatics.solver
         coulomb_params = coulomb.get_params()
         self.assertEqual(first=coulomb.name(),
                          second='Coulomb::CoulombP3M',
@@ -227,10 +221,7 @@ def test_setup_electrostatics(self):
         self.assertEqual(first=electrostatics_inputs["tune_p3m"],
                          second=coulomb_params["is_tuned"],
                          msg="lib.handy_functions.setup_electrostatic_interactions does not tune the P3M method")
-        if espressomd.version.friendly() == "4.2":
-            espresso_system.actors.remove(coulomb)
-        else:
-            coulomb = espresso_system.electrostatics.solver = None
+        coulomb = espresso_system.electrostatics.solver = None
         ## Test the setup of the P3M method without tuning it with some input parameters
         electrostatics_inputs["tune_p3m"] = False
         electrostatics_inputs["params"] = {"mesh": [8, 8, 8], 
@@ -238,30 +229,21 @@ def test_setup_electrostatics(self):
                                            "alpha": 1.1265e+01,
                                            "r_cut": 1}
         pmb.simulation_engine.setup_electrostatic_interactions(**electrostatics_inputs)
-        if espressomd.version.friendly() == "4.2":
-            coulomb = espresso_system.actors.active_actors.copy()[0]
-        else:
-            coulomb = espresso_system.electrostatics.solver
+        coulomb = espresso_system.electrostatics.solver
         coulomb_params = coulomb.get_params()
         for param in electrostatics_inputs["params"]:
             np.testing.assert_allclose(
                  np.copy(electrostatics_inputs["params"][param]),
                  np.copy(coulomb_params[param]),
                  err_msg="lib.handy_functions.setup_electrostatic_interactions sets up the wrong P3M parameters")
-        if espressomd.version.friendly() == "4.2":
-            espresso_system.actors.remove(coulomb)
-        else:
-            coulomb = espresso_system.electrostatics.solver = None
+        coulomb = espresso_system.electrostatics.solver = None
         electrostatics_inputs["params"] = None
         # Test the Debye–Hückel setup
         electrostatics_inputs["method"] = "dh"
         electrostatics_inputs["c_salt"] = pmb.units.Quantity(1, "mol/L")
         kappa=1./np.sqrt(8*pmb.units.pi*Bjerrum_length*pmb.N_A*electrostatics_inputs["c_salt"])
         pmb.simulation_engine.setup_electrostatic_interactions(**electrostatics_inputs)
-        if espressomd.version.friendly() == "4.2":
-            dh = espresso_system.actors.active_actors.copy()[0]
-        else:
-            dh = espresso_system.electrostatics.solver
+        dh = espresso_system.electrostatics.solver
         dh_params = dh.get_params()
         self.assertEqual(first=dh.name(),
                          second='Coulomb::DebyeHueckel',
@@ -275,16 +257,10 @@ def test_setup_electrostatics(self):
         self.assertAlmostEqual(first=dh_params["r_cut"],
                                 second=3*kappa.m_as('reduced_length'),
                                 msg="lib.handy_functions.setup_electrostatic_interactions sets up the wrong cut-off for the DH method")
-        if espressomd.version.friendly() == "4.2":
-            espresso_system.actors.remove(dh)
-        else:
-            coulomb = espresso_system.electrostatics.solver = None
+        coulomb = espresso_system.electrostatics.solver = None
         electrostatics_inputs["c_salt"] = pmb.units.Quantity(1, "mol/L")*pmb.N_A
         pmb.simulation_engine.setup_electrostatic_interactions(**electrostatics_inputs)
-        if espressomd.version.friendly() == "4.2":
-            dh = espresso_system.actors.active_actors.copy()[0]
-        else:
-            dh = espresso_system.electrostatics.solver
+        dh = espresso_system.electrostatics.solver
         dh_params = dh.get_params()
         self.assertAlmostEqual(first=dh_params["kappa"],
                                 second=(1./kappa).m_as('1/ reduced_length'),
@@ -292,17 +268,11 @@ def test_setup_electrostatics(self):
         self.assertAlmostEqual(first=dh_params["r_cut"],
                                 second=3*kappa.m_as('reduced_length'),
                                 msg="lib.handy_functions.setup_electrostatic_interactions sets up the wrong cut-off for the DH method")
-        if espressomd.version.friendly() == "4.2":
-            espresso_system.actors.remove(dh)
-        else:
-            coulomb = espresso_system.electrostatics.solver = None
+        coulomb = espresso_system.electrostatics.solver = None
         # Test a non-default cut-off
         electrostatics_inputs["params"] = {"r_cut": 3}
         pmb.simulation_engine.setup_electrostatic_interactions(**electrostatics_inputs)
-        if espressomd.version.friendly() == "4.2":
-            dh = espresso_system.actors.active_actors.copy()[0]
-        else:
-            dh = espresso_system.electrostatics.solver
+        dh = espresso_system.electrostatics.solver
         dh_params = dh.get_params()
         self.assertAlmostEqual(first=dh_params["r_cut"],
                                 second=electrostatics_inputs["params"]["r_cut"],