pm-blanco · Aug 28, 2026 · Aug 26, 2026 · Aug 26, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -68,8 +68,8 @@ git clone [email protected]: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
@@ -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
@@ -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
@@ -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,