Increasing pyMBE interoperability (part I)

pyMBE/exceptions/pmb_warnings.py

@@ -0,0 +1,17 @@
+import warnings
+from functools import wraps
+
+def deprecated(new_function):
+    """Wrapper function to display a warning deprecation message to all functions that are going to be removed in a future
+    """
+    def decorator(func):
+        @wraps(func)
+        def wrapper(*args, **kwargs):
+            warnings.warn(
+                f"{func.__name__} is deprecated and it will be removed in the future. Please use"+ new_function +"instead.",
+                DeprecationWarning,
+                stacklevel=2
+            )
+            return func(*args, **kwargs)
+        return wrapper
+    return decorator

pyMBE/lib/handy_functions.py

@@ -20,6 +20,7 @@
 import re
 import numpy as np
 import scipy
+from pyMBE.exceptions.pmb_warnings import deprecated
 
 def calculate_initial_bond_length(bond_parameters, bond_type, lj_parameters):
     """
@@ -281,7 +282,7 @@ def define_peptide_AA_residues(sequence,model, pmb):
             pmb.define_residue(name = residue_name, 
                                     central_bead = central_bead,
                                     side_chains = side_chains)
-
+@deprecated('pmb.simulation_engine.do_reaction')
 def do_reaction(algorithm, steps):
     """
     Executes reaction steps using an ESPResSo reaction algorithm with
@@ -311,6 +312,7 @@ def do_reaction(algorithm, steps):
     else:
         algorithm.reaction(steps=steps)
 
+@deprecated('pmb.simulation_engine.get_number_of_particles')
 def get_number_of_particles(espresso_system, ptype):
     """
     Returns the number of particles of a given ESPResSo particle type.
@@ -508,7 +510,7 @@ def protein_sequence_parser(sequence):
             clean_sequence.append(residue_ok)
     return clean_sequence
 
-
+@deprecated('pmb.simulation_engine.relax_system')
 def relax_espresso_system(espresso_system, 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:
@@ -571,7 +573,7 @@ def relax_espresso_system(espresso_system, seed, gamma=1e-3, Nsteps_steepest_des
     logging.info(f"*** Minimum particle distance after relaxation: {espresso_system.analysis.min_dist()} ***")
     logging.debug("*** Relaxation finished ***")
     return espresso_system.analysis.min_dist()
-
+@deprecated('pmb.simulation_engine.setup_langevin_dynamics')
 def setup_langevin_dynamics(espresso_system, kT, seed,time_step=1e-2, gamma=1, tune_skin=True, min_skin=1, max_skin=None, tolerance=1e-3, int_steps=200, adjust_max_skin=True):
     """
     Sets up Langevin Dynamics for an ESPResSo simulation system.
@@ -634,7 +636,7 @@ def setup_langevin_dynamics(espresso_system, kT, seed,time_step=1e-2, gamma=1, t
                                               int_steps=int_steps, 
                                               adjust_max_skin=adjust_max_skin)
         logging.info(f"Optimized skin value: {espresso_system.cell_system.skin}")
-
+@deprecated('pmb.simulation_engine.setup_electrostatic_interactions')
 def setup_electrostatic_interactions(units, espresso_system, kT, c_salt=None, solvent_permittivity=78.5, method='p3m', tune_p3m=True, accuracy=1e-3, params=None, verbose=False):
     """
     Sets up electrostatic interactions in an ESPResSo system.

pyMBE/lib/lattice.py

@@ -333,5 +333,5 @@ def __init__(self,mpc,bond_l):
             raise ValueError("mpc must be a non-zero positive integer.")
         self.mpc = mpc
         self.bond_l = bond_l
-        self.box_l = (self.mpc+2)*self.bond_l.magnitude / (np.sqrt(3)*0.25)
+        self.box_l = (self.mpc+1)*self.bond_l.magnitude / (np.sqrt(3)*0.25)       
 

pyMBE/simulation_builder/base_engine.py

@@ -0,0 +1,144 @@
+from abc import ABC,abstractmethod
+import numpy as np
+
+class SimulationEngine(ABC):
+    """Base Class for simulation engines contains methods """
+    def __init__(self):
+        pass
+    @abstractmethod
+    def _check_bond_inputs(self):
+        return
+    @abstractmethod
+    def _create_bond_instance(self):
+        return
+    @abstractmethod
+    def _get_bond_instance(self):
+        return
+    @abstractmethod
+    def add_instances_to_engine(self):
+        return
+    def calculate_center_of_mass(self, instance_id, pmb_type):
+        """
+        Calculates the center of mass of a pyMBE object instance in an ESPResSo system.
+
+        Args:
+            instance_id ('int'):
+                pyMBE instance ID of the object whose center of mass is calculated.
+
+            pmb_type ('str'):
+                Type of the pyMBE object. Must correspond to a particle-aggregating
+                template type (e.g. '"molecule"', '"residue"', '"peptide"', '"protein"').
+
+        Returns:
+            ('numpy.ndarray'):
+                Array of shape '(3,)' containing the Cartesian coordinates of the
+                center of mass.
+
+        Notes:
+            - This method assumes equal mass for all particles.
+            - Periodic boundary conditions are *not* unfolded; positions are taken
+            directly from ESPResSo particle coordinates.
+        """
+        center_of_mass = np.zeros(3)
+        axis_list = [0,1,2]
+        inst = self.db.get_instance(pmb_type=pmb_type,
+                                    instance_id=instance_id)
+        particle_id_list = self.db.get_particle_id_map(object_name=inst.name)["all"]
+        for pid in particle_id_list:
+            for axis in axis_list:
+                center_of_mass[axis] += self.db.get_instance(pmb_type='particle',
+                                    instance_id=pid).position[axis]
+        center_of_mass = center_of_mass / len(particle_id_list)
+        return center_of_mass
+    def determine_reservoir_concentrations(self, pH_res, c_salt_res, activity_coefficient_monovalent_pair, max_number_sc_runs=200):
+        """
+        Determines ionic concentrations in the reservoir at fixed pH and salt concentration.
+
+        Args:
+            pH_res ('float'):
+                Target pH value in the reservoir.
+
+            c_salt_res ('pint.Quantity'):
+                Concentration of monovalent salt (e.g., NaCl) in the reservoir.
+
+            activity_coefficient_monovalent_pair ('callable'):
+                Function returning the activity coefficient of a monovalent ion pair
+                as a function of ionic strength:
+                'gamma = activity_coefficient_monovalent_pair(I)'.
+
+            max_number_sc_runs ('int', optional):
+                Maximum number of self-consistent iterations allowed before
+                convergence is enforced. Defaults to 200.
+
+        Returns:
+            tuple:
+                (cH_res, cOH_res, cNa_res, cCl_res)
+                - cH_res ('pint.Quantity'): Concentration of H⁺ ions.
+                - cOH_res ('pint.Quantity'): Concentration of OH⁻ ions.
+                - cNa_res ('pint.Quantity'): Concentration of Na⁺ ions.
+                - cCl_res ('pint.Quantity'): Concentration of Cl⁻ ions.
+
+        Notess:
+            - The algorithm enforces electroneutrality in the reservoir.
+            - Water autodissociation is included via the equilibrium constant 'Kw'.
+            - Non-ideal effects enter through activity coefficients depending on
+            ionic strength.
+            - The implementation follows the self-consistent scheme described in
+            Landsgesell (PhD thesis, Sec. 5.3, doi:10.18419/opus-10831), adapted
+            from the original code (doi:10.18419/darus-2237).
+        """
+        def determine_reservoir_concentrations_selfconsistently(cH_res, c_salt_res):
+            """
+            Iteratively determines reservoir ion concentrations self-consistently.
+
+            Args:
+                cH_res ('pint.Quantity'):
+                    Current estimate of the H⁺ concentration.
+                c_salt_res ('pint.Quantity'):
+                    Concentration of monovalent salt in the reservoir.
+
+            Returns:
+                'tuple':
+                    (cH_res, cOH_res, cNa_res, cCl_res)
+            """
+            # Initial ideal estimate
+            cOH_res = self.Kw / cH_res
+            if cOH_res >= cH_res:
+                cNa_res = c_salt_res + (cOH_res - cH_res)
+                cCl_res = c_salt_res
+            else:
+                cCl_res = c_salt_res + (cH_res - cOH_res)
+                cNa_res = c_salt_res
+            # Self-consistent iteration
+            for _ in range(max_number_sc_runs):
+                ionic_strength_res = 0.5 * (cNa_res + cCl_res + cOH_res + cH_res)
+                cOH_new = self.Kw / (cH_res * activity_coefficient_monovalent_pair(ionic_strength_res))
+                if cOH_new >= cH_res:
+                    cNa_new = c_salt_res + (cOH_new - cH_res)
+                    cCl_new = c_salt_res
+                else:
+                    cCl_new = c_salt_res + (cH_res - cOH_new)
+                    cNa_new = c_salt_res
+                # Update values
+                cOH_res = cOH_new
+                cNa_res = cNa_new
+                cCl_res = cCl_new
+            return cH_res, cOH_res, cNa_res, cCl_res
+        # Initial guess for H+ concentration from target pH
+        cH_res = 10 ** (-pH_res) * self.units.mol / self.units.l
+        # First self-consistent solve
+        cH_res, cOH_res, cNa_res, cCl_res = (determine_reservoir_concentrations_selfconsistently(cH_res, 
+                                                                                                 c_salt_res))
+        ionic_strength_res = 0.5 * (cNa_res + cCl_res + cOH_res + cH_res)
+        determined_pH = -np.log10(cH_res.to("mol/L").magnitude* np.sqrt(activity_coefficient_monovalent_pair(ionic_strength_res)))
+        # Outer loop to enforce target pH
+        while abs(determined_pH - pH_res) > 1e-6:
+            if determined_pH > pH_res:
+                cH_res *= 1.005
+            else:
+                cH_res /= 1.003
+            cH_res, cOH_res, cNa_res, cCl_res = (determine_reservoir_concentrations_selfconsistently(cH_res, 
+                                                                                                     c_salt_res))
+            ionic_strength_res = 0.5 * (cNa_res + cCl_res + cOH_res + cH_res)
+            determined_pH = -np.log10(cH_res.to("mol/L").magnitude * np.sqrt(activity_coefficient_monovalent_pair(ionic_strength_res)))
+        return cH_res, cOH_res, cNa_res, cCl_res

pyMBE/simulation_builder/engine_protocol.py

@@ -0,0 +1,25 @@
+from typing import Protocol,runtime_checkable
+
+class EspressoParticleProtocol(Protocol):
+    """Class that emulates the estructure of the Espresso Particle class"""
+    def add():
+        return
+    def by_id():
+        return
+class EspressoBondedInterProtocol(Protocol):
+    """Class that emulates the structure of the EspressoBondedInterProtocol"""
+    def add():
+        return
+
+@runtime_checkable
+class EspressoSystemProtocol(Protocol):
+    """ Class that emulates the structure of the methods employed by Pymbe from the espressomd.System class
+        . The decorator @runtime_checkable allows to only check for the structure not the types"""
+    part: EspressoParticleProtocol
+    bonded_inter: EspressoBondedInterProtocol
+
+@runtime_checkable
+class LammpsProtocol(Protocol):
+    """ Class that emulates the structure of the methods employed by Pymbe from the Lammps class
+        . The decorator @runtime_checkable allows to only check for the structure not the types"""
+