Drop support for ESPResSo 4.2

Author: jngrad

.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

.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)

.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:

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

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",

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,

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