Add a hydrogel builder

CHANGELOG.md

@@ -19,7 +19,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Switched to CTest for testing, allowing to run the tests on paralel (#87)
 
 ### Added
-- Optional argument `params` in `lib.handy_functions.setup_electrostatic_interactions` enabling the user to easilly tune parameters for the electrostatic solvers. (#121)
+- New benchmark for hydrogels, including scripts to reproduce the data `samples/Landsgesell2022/run_simulations.py` and `samples/Landsgesell2022/plot_pH_vs_alpha.py` and `samples/Landsgesell2022/plot_P_vs_V.py` (#103)
+- New sample scripts for hydrogels `samples/build_hydrogel.py` and  `samples/weak_polyacid_hydrogel_grxmc.py` (#103)
+- New methods to support building hydrogels with pyMBE `pmb.define_hydrogel`, `pmb.create_hydrogel`, `pmb.initialize_lattice_builder`, `pmb.create_hydrogel_chain`, `pmb.create_hydrogel_node`. (#103)
 - CI testing for functions in `lib.handy_functions`. (#118)
 - sanity tests for `lib.handy_functions`. (#118)
 - Use of `logging`  in `lib.handy_functions` to handle output and error logs. (#118)

README.md

@@ -12,7 +12,7 @@
 [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](LICENSE.txt)
 [![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.1-4baaaa.svg)](CODE_OF_CONDUCT.md) 
 
-pyMBE provides tools to facilitate building up molecules with complex architectures in the Molecular Dynamics software [ESPResSo](https://espressomd.org/wordpress/). Some examples of molecules that can be set up with pyMBE are polyelectrolytes, peptides and proteins. pyMBE bookkeeps all the information about the molecule topology, permitting to link each particle to its corresponding residue and molecule. pyMBE uses the [Pint](https://pint.readthedocs.io/en/stable/) library to enable input parameters in any arbitrary unit system, which is later transformed in the reduced unit system used in ESPResSo.
+pyMBE provides tools to facilitate building up molecules with complex architectures in the Molecular Dynamics software [ESPResSo](https://espressomd.org/wordpress/). Some examples of molecules that can be set up with pyMBE are polyelectrolytes, hydrogels, peptides and globular proteins. pyMBE bookkeeps all the information about the molecule topology, permitting to link each particle to its corresponding residue and molecule. pyMBE uses the [Pint](https://pint.readthedocs.io/en/stable/) library to enable input parameters in any arbitrary unit system, which is later transformed in the reduced unit system used in ESPResSo.
 
 An up-to-date documentation of all methods of the library can be found [here](pymbe-dev.github.io/pyMBE/pyMBE.html) and in the source code.
 

lib/analysis.py

@@ -153,27 +153,33 @@ def block_analyze(full_data, n_blocks=16, time_col = "time", equil=0.1,  columns
 
     return result
 
+import re
+
 def built_output_name(input_dict):
     """
     Builts the output name for a given set of input parameters.
-
+    
     Args:
         input_dict(`dict`): dictionary with all terminal inputs.
 
     Returns:
         output_name(`str`): name used for the output files
-
+    
     Note:
         - The standard formatting rule is parametername1_parametervalue1_parametername2_parametervalue2
     """
     output_name=""
     for label in input_dict:
-        if type(input_dict[label]) in [str,bool]:
-            formatted_variable=f"{input_dict[label]:}"
-        else:
-            formatted_variable=f"{input_dict[label]:.3g}"
-        output_name+=f"{label}_{formatted_variable}_"
-    return output_name[:-1]
+        if isinstance(input_dict[label], (str, bool)):
+            formatted_variable = f"{input_dict[label]}"
+        else:  
+            formatted_variable = f"{input_dict[label]:.3g}"
+
+        formatted_variable = formatted_variable.replace(" ", "_")
+        formatted_variable = re.sub(r"[^\w\-.]", "", formatted_variable)
+        output_name += f"{label}_{formatted_variable}_"
+
+    return output_name.rstrip("_")  # Remove trailing underscore
 
 def get_dt(data, time_col = "time", relative_tolerance = 0.01, verbose = False):
     """ 

lib/handy_functions.py

@@ -85,6 +85,7 @@ def setup_electrostatic_interactions(units, espresso_system, kT, c_salt=None, so
             else:
                 espresso_system.electrostatics.solver = coulomb
 
+
             # save the optimal parameters and add them by hand
 
             p3m_params = coulomb.get_params()

lib/lattice.py

@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2024 pyMBE-dev team
+# Copyright (C) 2024-2025 pyMBE-dev team
 #
 # This file is part of pyMBE.
 #
@@ -15,22 +15,18 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
 
 import itertools
 import numpy as np
 
-
 class LatticeBuilder:
     """
     Generic lattice builder.
-
     Args:
         lattice(`object`): Data class that represents a lattice, for example a :class:`DiamondLattice`.
         strict(`bool`): If set to `True`, the lattice connectivity cannot be altered.
             For example, a chain cannot be created between two nodes that are
             not explicitly connected according to the definition in `lattice`.
-
     Attributes:
         kwargs_node_labels(`dict`): Keyword arguments passed to the matplotlib
             `Axes3D.text() <https://matplotlib.org/stable/api/_as_gen/mpl_toolkits.mplot3d.axes3d.Axes3D.text.html>`_
@@ -51,13 +47,17 @@ def __init__(self, lattice, strict=True):
         self.strict = strict
         self.node_labels = {str([int(x) for x in indices]).replace(",", ""): i
                             for i, indices in enumerate(self.lattice.indices)}
+        self.chain_labels = {f"({start}, {end})": idx for idx, (start, end) in enumerate(lattice.connectivity)}
         self.nodes = {label: "default_linker" for label in self.node_labels}
         self.colormap = {}
         self.chains = {}
         self.kwargs_node_labels = {}
         self.kwargs_monomers = {}
         self.kwargs_bonds = {}
         self.kwargs_box = {}
+        self.mpc = lattice.mpc
+        self.box_l = lattice.box_l
+
 
     def _get_node_by_label(self, node):
         assert node in self.node_labels, f"node '{node}' doesn't exist in a {self.lattice.name} lattice"
@@ -93,7 +93,6 @@ def _make_sphere(cls, radius, resolution):
     def add_default_chains(self, mpc):
         """
         Build a default lattice network. Skip node pairs that already have a chain defined.
-
         Args:
             mpc(`int`): Number of monomers per chain.
         """
@@ -104,7 +103,6 @@ def add_default_chains(self, mpc):
     def draw_lattice(self, ax):
         """
         Draw the lattice in an `Axes3D <https://matplotlib.org/stable/api/toolkits/mplot3d/axes3d.html>`_ canvas.
-
         Args:
             ax: Axes.
         """
@@ -150,8 +148,10 @@ def draw_lattice(self, ax):
             if self.colormap:
                 color = self.colormap[node_type]
                 kwargs_monomers["c"] = color
+
             node_positions = scatter_data[node_type]
             pos = np.array(node_positions)
+            # plotting nodes and monomers
             ax_data = ax.scatter(pos[:,0], pos[:,1], pos[:,2], edgecolor="none",
                                  zorder=2, label=node_type, s=12**2, **kwargs_monomers)
             color = ax_data.get_facecolors()[0]
@@ -163,7 +163,6 @@ def draw_lattice(self, ax):
     def draw_simulation_box(self, ax):
         """
         Draw the simulation box in an `Axes3D <https://matplotlib.org/stable/api/toolkits/mplot3d/axes3d.html>`_ canvas.
-
         Args:
             ax: Axes.
         """
@@ -182,11 +181,9 @@ def draw_simulation_box(self, ax):
     def get_chain(self, node_start, node_end):
         """
         Get a chain between two nodes.
-
         Args:
             node_start(`str`): Start node label, e.g. ``[0 0 0]`` for the node at the origin.
             node_end(`str`): End node label, e.g. ``[1 1 1]`` for the first node along the diagonal.
-
         Returns:
             `list`: Sequence of residue labels.
         """
@@ -202,10 +199,8 @@ def get_monomer_color(self, label):
         """
         Get the color corresponding to a monomer label.
         Only works if a custom color map was set via :meth:`set_colormap`!
-
         Args:
             label(`str`): Monomer label.
-
         Returns:
             The color.
         """
@@ -216,10 +211,8 @@ def get_monomer_color_index(self, label):
         """
         Get the color wheel index corresponding to a monomer label.
         Only works if a custom color map was set via :meth:`set_colormap`!
-
         Args:
             label(`str`): Monomer label.
-
         Returns:
             `int`: The color index.
         """
@@ -230,56 +223,24 @@ def get_monomer_color_index(self, label):
     def get_node(self, node):
         """
         Get a node residue label.
-
         Args:
             node(`str`): Node label, e.g. ``[0 0 0]`` for the node at the origin.
-
         Returns:
             `str`: Residue label.
         """
         key = self._get_node_by_label(node)
         return self.nodes[key]
 
-    def set_chain(self, node_start, node_end, sequence):
-        """
-        Set a chain between two nodes.
-
-        Args:
-            node_start(`str`): Start node label, e.g. ``[0 0 0]`` for the node at the origin.
-            node_end(`str`): End node label, e.g. ``[1 1 1]`` for the first node along the diagonal.
-            sequence(`list`): Sequence of residue labels.
-        """
-        assert len(sequence) != 0 and not isinstance(sequence, str)
-        key, reverse = self._get_node_vector_pair(node_start, node_end)
-        assert node_start != node_end or sequence == sequence[::-1], \
-            (f"chain cannot be defined between '{node_start}' and '{node_end}' since it "
-            "would form a loop with a non-symmetric sequence (under-defined stereocenter)")
-        if reverse:
-            sequence = sequence[::-1]
-        self.chains[key] = sequence
-
     def set_colormap(self, colormap):
         """
         Set a discrete color map. By default, the standard matplotlib color wheel will be used.
-
         Args:
             colormap(`dict`): Discrete color wheel that maps monomer labels to colors.
         """
         assert isinstance(colormap, dict)
         self.colormap = colormap.copy()
         self.colormap_sorted_names = list(self.colormap.keys())
 
-    def set_node(self, node, residue):
-        """
-        Set a node residue type.
-
-        Args:
-            node(`str`): Node label, e.g. ``[0 0 0]`` for the node at the origin.
-            residue(`str`): Node residue label.
-        """
-        key = self._get_node_by_label(node)
-        self.nodes[key] = residue
-
 
 class DiamondLattice:
     """
@@ -294,3 +255,11 @@ class DiamondLattice:
                     (2, 5), (3, 6), (4, 7), (5, 0),
                     (5, 3), (5, 4), (6, 0), (6, 2),
                     (6, 4), (7, 0), (7, 2), (7, 3)}
+
+    def __init__(self,mpc,bond_l):
+        if not isinstance(mpc, int) or mpc <= 0:
+            raise ValueError("mpc must be a non-zero positive integer.")
+        self.mpc = mpc
+        self.bond_l = bond_l
+        self.box_l = (self.mpc+1)*self.bond_l.magnitude / (np.sqrt(3)*0.25)
+

maintainer/standarize_data.py

@@ -31,7 +31,8 @@
                      "Lys-AspMSDE.csv",
                      "histatin5_SoftMatter.txt",
                      "1beb-10mM-torres.dat",
-                     "1f6s-10mM-torres.dat"]
+                     "1f6s-10mM-torres.dat",
+                     "landsgesell2022"]
 
 
 parser = argparse.ArgumentParser(description='Script to standarize the data from various authors')
@@ -50,40 +51,46 @@
                   "Glu-HisMSDE.csv": "Lunkad2021b.csv",
                   "histatin5_SoftMatter.txt": "Blanco2020a.csv",
                   "1beb-10mM-torres.dat": "Torres2017.csv",
-                  "1f6s-10mM-torres.dat": "Torres2022.csv"}
+                  "1f6s-10mM-torres.dat": "Torres2022.csv",
+                  "landsgesell2022": "Landsgesell2022a.csv"}
 
 # Sanity checks
 if filename not in supported_filenames:
     ValueError(f"Filename {filename} not supported, supported files are {supported_filenames}")
 
-# Extract the data from Ref.
-ref_path=pmb.get_resource(f"testsuite/data/src/{filename}")
-Refs_lunkad=["Glu-HisMSDE.csv","Lys-AspMSDE.csv"]
-Ref_blanco=["histatin5_SoftMatter.txt"]
+files_lunkad=["Glu-HisMSDE.csv","Lys-AspMSDE.csv"]
+files_blanco=["histatin5_SoftMatter.txt"]
+files_torres = ["1f6s-10mM-torres.dat","1beb-10mM-torres.dat" ]
+files_landsgesell2020=["data_landsgesell.csv"]
+files_landsgesell2022=["equilibrium_values_gel_MD.csv","weak-gel_total_data.csv"]
+
+# Get path to the data from publication
+if filename == "landsgesell2022":
+    ref_paths=[]
+    for file in files_landsgesell2022:
+        ref_paths.append(pmb.get_resource(f"testsuite/data/src/{file}"))
+else:
+    ref_path=pmb.get_resource(f"testsuite/data/src/{filename}")
 
-Ref_torres = ["1f6s-10mM-torres.dat","1beb-10mM-torres.dat" ]
-Ref_landsgesell=["data_landsgesell.csv"]
 
-if filename in Refs_lunkad:
+if filename in files_lunkad:
     data=pd.read_csv(ref_path)
     Z_ref = 5*-1*data['aaa']+5*data['aab']
     # Error propagation calculation
     # 1/4 factor added to correct for bug in the original calculation of the error reported by the authors       
     Z_ref_err = 5/4*np.sqrt((data['eaa'])**2+(data['eab'])**2)   
     pH_range = np.linspace(2, 12, num=21)
 
-elif filename in Ref_blanco:
+elif filename in files_blanco:
     data=np.loadtxt(ref_path, delimiter=",")
     Z_ref=data[:,1]         
     Z_ref_err=data[:,2]
     pH_range = np.linspace(2, 12, num=21)
 
-elif filename in Ref_torres:
-
+elif filename in files_torres:
     Z_ref = []
     Z_ref_err = []
     pH_range = []
-
     with open (ref_path,'r') as file:
         for line in file: 
             line_split = line.split ()
@@ -92,12 +99,18 @@
             Z_ref.append (float(line_split[1]))
             Z_ref_err.append(float(line_split[2]))
 
-elif filename in Ref_landsgesell:
+elif filename in files_landsgesell2020:
     data = pd.read_csv(ref_path, sep="\t", index_col=False)
+elif filename == "landsgesell2022":
+    dfs = []
+    for path in ref_paths:
+        dfs.append(pd.read_csv(path))
+    data = pd.concat(dfs, axis=0, ignore_index=True)
+
 else:
     raise RuntimeError()
 
-if filename in Refs_lunkad+Ref_blanco+Ref_torres:
+if filename in files_lunkad+files_blanco+files_torres:
     # Create the pandas DataFrame
     data=pd.DataFrame({"pH": pH_range,
                       "charge": Z_ref,

parameters/salt/excess_chemical_potential_excess_pressure.csv

@@ -0,0 +1,26 @@
+cs_bulk_[1/sigma^3],excess_chemical_potential_[kbT],bjerrum_length,error_excess_chemical_potential_[kbT],excess_pressure_[kT/sigma^3],std_error_excess_pressure_[kT/sigma^3]
+2.695417789759999814e-09,-4.351021137877512848e-05,2.000000000000000000e+00,6.921676554563794465e-07,-3.851827594318893624e-14,2.108526448153898900e-16
+8.523659461370000089e-09,-7.427122674414749989e-05,2.000000000000000000e+00,4.493516326433293331e-06,-1.807616112580177201e-13,1.146671935590645798e-15
+2.695417789760000145e-08,-1.035372639496550149e-04,2.000000000000000000e+00,5.063677845596300503e-06,-8.171432597348823871e-13,1.654129333888704547e-14
+2.695417789759999999e-06,-3.031940271292494007e-02,2.000000000000000000e+00,1.413240857466053363e-04,-2.636400004572970776e-08,2.636031040805764856e-11
+8.523659461370000539e-06,-4.127796588505479314e-02,2.000000000000000000e+00,2.449703657160429604e-04,-1.124972002996681094e-07,1.474783738054616950e-10
+2.695417789760000507e-05,-7.615421877488305624e-02,2.000000000000000000e+00,5.713996171758657476e-04,-6.566989095396179666e-07,1.170514215713901519e-09
+8.523659461370000200e-05,-1.252709992055541510e-01,2.000000000000000000e+00,6.564159273489764005e-04,-3.386738439846527169e-06,4.837470792608095765e-09
+2.695417699999999995e-04,-2.072884574302087946e-01,2.000000000000000000e+00,1.106112782639213227e-03,-1.745661018404891484e-05,2.851032310417292497e-08
+2.695417789760000101e-04,-2.085397242843187171e-01,2.000000000000000000e+00,1.076965644017058039e-03,-1.742957410003485755e-05,2.847577671665599596e-08
+5.390835499999999591e-04,-2.784248888038032832e-01,2.000000000000000000e+00,1.375397742180195981e-03,-4.550204871694554966e-05,7.791469246477961844e-08
+8.523659461370000200e-04,-3.365753246398797693e-01,2.000000000000000000e+00,1.590969027075880404e-03,-8.447465025296207596e-05,1.468765130415688624e-07
+1.347708889999999892e-03,-4.011228563537788672e-01,2.000000000000000000e+00,1.854691962095065616e-03,-1.536430422354193373e-04,2.881141243655981952e-07
+2.695417779999999784e-03,-5.081652726283076849e-01,2.000000000000000000e+00,2.422843717875428852e-03,-3.619845221956632117e-04,8.019199244451515603e-07
+5.390835570000000396e-03,-6.204356667380437340e-01,2.000000000000000000e+00,3.754857762950025105e-03,-7.862234220488582750e-04,2.633689671874538852e-06
+8.523659461370000634e-03,-6.738392222804889808e-01,2.000000000000000000e+00,4.207216501142822286e-03,-1.200338737953249085e-03,5.077061434219377856e-06
+9.973045822100000748e-03,-6.883749507113008370e-01,2.000000000000000000e+00,4.142739213683488755e-03,-1.347949799578286589e-03,5.757014473469371011e-06
+1.185983827489999962e-02,-6.998672021084888506e-01,2.000000000000000000e+00,5.425017468960176108e-03,-1.512824374308426555e-03,9.371259478143074167e-06
+1.374663072780000009e-02,-7.198982053501644662e-01,2.000000000000000000e+00,4.554926564004751696e-03,-1.616960989778337601e-03,8.968066952719475849e-06
+1.563342318060000069e-02,-7.154785956112332812e-01,2.000000000000000000e+00,6.054499564496755408e-03,-1.671669149869207732e-03,1.263119223688707693e-05
+1.752021563340000129e-02,-7.097015018783637830e-01,2.000000000000000000e+00,5.148834780668647484e-03,-1.629480842783949173e-03,1.342254352194756690e-05
+1.940700808629999830e-02,-7.298771071109483310e-01,2.000000000000000000e+00,6.608397430862909816e-03,-1.580985612615371221e-03,1.759721758743890897e-05
+2.129380053909999890e-02,-7.084622937719967650e-01,2.000000000000000000e+00,7.753353976806658283e-03,-1.389621207906211950e-03,2.121427526850919976e-05
+2.318059299189999950e-02,-7.002012770913444983e-01,2.000000000000000000e+00,6.554812199665331282e-03,-1.209649228197008956e-03,2.227643781670899859e-05
+2.506738544470000010e-02,-6.889219784851126072e-01,2.000000000000000000e+00,6.871396583625398111e-03,-1.003024820532205695e-03,2.468442377998216467e-05
+2.695417789760000057e-02,-6.609159955081040927e-01,2.000000000000000000e+00,6.475734836756414509e-03,-6.696675815468191477e-04,2.506221438124657882e-05