new testuites for hydrogel builder

lib/lattice.py

@@ -260,6 +260,8 @@ class DiamondLattice:
                     (6, 4), (7, 0), (7, 2), (7, 3)}
 
     def __init__(self,MPC,BOND_LENGTH):
+        if not isinstance(MPC, int) or MPC <= 0:
+            raise ValueError("MPC must be a non-zero positive integer.")
         self.MPC = MPC
         self.BOND_LENGTH = BOND_LENGTH
         self.BOXL = (self.MPC+1)*self.BOND_LENGTH.magnitude / (np.sqrt(3)*0.25)

pyMBE.py

@@ -25,7 +25,7 @@
 import pandas as pd
 import scipy.constants
 import scipy.optimize
-from lib.lattice import LatticeBuilder
+from lib.lattice import LatticeBuilder, DiamondLattice
 
 
 class pymbe_library():
@@ -891,26 +891,30 @@ def format_node(node_list):
         self.check_if_name_is_defined_in_df(name=name,
                                             pmb_type_to_be_defined='hydrogel')
 
+        hydrogel_info={"name":name, "chains":{}, "nodes":{}}
         # placing nodes
         node_positions = {}
         node_topology = self.df[self.df["name"]==name]["node_map"].iloc[0]
         for node_info in node_topology.values():
             node_index = node_info["lattice_index"]
             node_name = node_info["particle_name"]
-
-            node_pos = self.set_node(format_node(node_index), node_name, espresso_system)
+            node_pos, node_id = self.set_node(format_node(node_index), node_name, espresso_system)
             node_label = self.lattice_builder.node_labels[format_node(node_index)]
+            hydrogel_info["nodes"][format_node(node_index)]=node_id
             node_positions[node_label]=node_pos
-
+        
         # Placing chains between nodes
         # Looping over all the 16 chains
         chain_topology = self.df[self.df["name"]==name]["chain_map"].iloc[0]
         for chain_info in chain_topology.values():
             node_s = chain_info["node_start"]
             node_e = chain_info["node_end"]
-            self.set_chain(node_s, node_e, node_positions, espresso_system)
-
-        return node_positions;
+            molecule_info = self.set_chain(node_s, node_e, node_positions, espresso_system)
+            for molecule_id in molecule_info:
+                hydrogel_info["chains"][molecule_id] = molecule_info[molecule_id]
+                hydrogel_info["chains"][molecule_id]["node_start"]=node_s
+                hydrogel_info["chains"][molecule_id]["node_end"]=node_e
+        return hydrogel_info;
 
 
     def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False):
@@ -1515,7 +1519,22 @@ def define_hydrogel(self, name, node_map, chain_map):
             chain_map(`dict`): {"chain_id": {"node_start": , "node_end": , "residue_list": , ... }
 
         """
+        expected_node_labels = set(range(8))  # Node labels 0 through 7
+        actual_node_labels = set(node_map.keys())
+        if actual_node_labels != expected_node_labels:
+            print(
+            f"Incomplete hydrogel: node_map must contain exactly 8 node_labels (0 through 7). "
+            f"Found: {sorted(actual_node_labels)}"
+        )
 
+        # Check chain_map contains exactly 16 chain_ids (0 through 15)
+        expected_chain_ids = set(range(16))  # Chain IDs 0 through 15
+        actual_chain_ids = set(chain_map.keys())
+        if actual_chain_ids != expected_chain_ids:
+            print(
+            f"Incomplete hydrogel: chain_map must contain exactly 16 chain_ids (0 through 15). "
+            f"Found: {sorted(actual_chain_ids)}"
+        )
 
         if self.check_if_name_is_defined_in_df(name=name,pmb_type_to_be_defined='hydrogel'):
             return
@@ -2297,6 +2316,8 @@ def initialize_lattice_builder(self, diamond_lattice):
         """
         Initialize the lattice builder with the DiamondLattice object.
         """
+        if not isinstance(diamond_lattice, DiamondLattice):
+            raise TypeError("Currently only DiamondLattice objects are supported.")
         self.lattice_builder = LatticeBuilder(lattice=diamond_lattice)
         print(f"LatticeBuilder initialized with MPC={diamond_lattice.MPC} and BOXL={diamond_lattice.BOXL}")
         return self.lattice_builder
@@ -2740,7 +2761,7 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system):
                 backbone_vector=np.array(backbone_vector),
                 use_default_bond=False  # Use default bonds between monomers
             )
-
+       
         # Collecting ids of beads of the chain/molecule
         chain_ids = []
         residue_ids = []
@@ -2793,6 +2814,7 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system):
                         new_value=residue_ids[-1],
                         verbose=False,
                         overwrite=True)
+        return chain_molecule_info
 
     def set_node(self, node, residue, espresso_system):
         """
@@ -2807,14 +2829,14 @@ def set_node(self, node, residue, espresso_system):
             raise ValueError("LatticeBuilder is not initialized. Use `initialize_lattice_builder` first.")
 
         node_position = np.array(list(int(x) for x in node.strip('[]').split()))*0.25*self.lattice_builder.BOXL
-        self.create_particle(name = residue,
+        p_id = self.create_particle(name = residue,
                          espresso_system=espresso_system,
                          number_of_particles=1,
                          position = [node_position])
         key = self.lattice_builder._get_node_by_label(node)
         self.lattice_builder.nodes[key] = residue
 
-        return node_position.tolist()
+        return node_position.tolist(), p_id
 
     def set_particle_acidity(self, name, acidity=None, default_charge_number=0, pka=None, verbose=True, overwrite=True):
         """

samples/hydrogel.py

@@ -73,6 +73,8 @@
                               'residue_list':residue_list}
     chain_id+=1
 
+del chain_topology[0]
+print(chain_topology)
 pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology)
 node_positions = pmb.create_hydrogel("my_hydrogel", espresso_system)
 

testsuite/hydrogel_builder.py

@@ -0,0 +1,225 @@
+import numpy as np
+import random
+import lib.lattice
+import unittest as ut
+import matplotlib
+import matplotlib.pyplot as plt
+import pyMBE
+from lib.lattice import DiamondLattice
+import espressomd
+
+pmb = pyMBE.pymbe_library(seed=42)
+MPC = 0 
+BOND_LENGTH = 0.355 * pmb.units.nm
+
+print("*** Unit test: check that any objects are other than DiamondLattice passed to initialize_lattice_builder raises a NameError ***")
+np.testing.assert_raises(TypeError, pmb.initialize_lattice_builder, None)
+print("*** Unit test passed ***")
+# Define node particle
+NodeType = "node_type"
+pmb.define_particle(name=NodeType, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
+# define monomers
+BeadType1 = "C"
+pmb.define_particle(name=BeadType1, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
+BeadType2 = "M"
+pmb.define_particle(name=BeadType2, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
+
+Res1 = "res_1"
+pmb.define_residue(
+    name=Res1,  # Name of the residue
+    central_bead=BeadType1,  # Define the central bead name
+    side_chains=[]  # Assuming no side chains for the monomer
+)
+
+Res2 = "res_2"
+pmb.define_residue(
+    name=Res2,  # Name of the residue
+    central_bead=BeadType2,  # Define the central bead name
+    side_chains=[]  # Assuming no side chains for the monomer
+)
+
+generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')
+generic_bond_length = 0.355*pmb.units.nm
+HARMONIC_parameters = {'r_0'    : generic_bond_length,
+                      'k'      : generic_harmonic_constant}
+
+pmb.define_bond(bond_type = 'harmonic',
+                        bond_parameters = HARMONIC_parameters, particle_pairs = [[BeadType1, BeadType1],
+                                                                                 [BeadType1, BeadType2],
+                                                                                 [BeadType2, BeadType2]])
+pmb.define_bond(bond_type = 'harmonic',
+                        bond_parameters = HARMONIC_parameters, particle_pairs = [[NodeType, BeadType1],
+                                                                                 [NodeType, BeadType2]])
+
+print("*** Unit Test: check that only non-negative values of monomers per chain are allowed ***")
+
+with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."):
+    diamond_lattice = DiamondLattice(MPC, generic_bond_length)
+
+print("*** Unit Test passed ***")
+
+MPC=8
+diamond_lattice = DiamondLattice(MPC, generic_bond_length)
+espresso_system = espressomd.System(box_l = [diamond_lattice.BOXL]*3)
+pmb.add_bonds_to_espresso(espresso_system = espresso_system)
+lattice_builder = pmb.initialize_lattice_builder(diamond_lattice)
+
+# Setting up node topology
+indices = diamond_lattice.indices
+node_topology = {}
+
+for index in range(len(indices)):
+    node_topology[index]={"particle_name": NodeType,
+                          "lattice_index": indices[index]}
+
+# Setting up chain topology
+connectivity = diamond_lattice.connectivity
+node_labels = lattice_builder.node_labels
+reverse_node_labels = {v: k for k, v in node_labels.items()}
+connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) for i, j in connectivity}
+chain_topology = {}
+residue_list = [Res1]*(MPC//2) + [Res2]*(MPC//2)
+chain_id = 0
+for node_s, node_e in connectivity_with_labels:
+    chain_topology[chain_id]={'node_start':node_s,
+                              'node_end': node_e,
+                              'residue_list':residue_list}
+    chain_id+=1
+
+#last_chain_id = max(chain_topology.keys())
+#chain_topology[last_chain_id]['residue_list'] = chain_topology[last_chain_id]['residue_list'][::-1]
+
+last_chain_id = max(chain_topology.keys())
+chain_topology[last_chain_id]['residue_list'] = ["res_1" if i % 2 == 0 else "res_2" for i in range(len(residue_list))]
+
+pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology)
+
+print("*** Unit Test: check node map and chain map parameters ***")
+
+random_chain_id = random.choice(list(chain_topology.keys())) # Choose a random chain of 0-15
+# extract node_start, node_end and residue_list of random chain
+chain_data = chain_topology[random_chain_id] 
+node_start = chain_data["node_start"] 
+node_end = chain_data["node_end"]
+residue_list = chain_data['residue_list']
+# choosing random residue in the residue_list
+random_res_in_res_list = random.choice(list(residue_list))
+res_indices = [i for i, residue in enumerate(residue_list) if residue == random_res_in_res_list]
+chosen_index = random.choice(res_indices)
+# Expected res_id of the randomly chosen residue
+random_res_id = random_chain_id*(len(residue_list)) + chosen_index
+
+hydrogel_data = pmb.df[pmb.df["name"] == "my_hydrogel"]
+stored_chain_map = hydrogel_data["chain_map"].values[0]
+
+ut.TestCase().assertEqual(stored_chain_map[str(random_chain_id)]["node_start"], node_start)
+ut.TestCase().assertEqual(stored_chain_map[str(random_chain_id)]["node_end"], node_end)
+ut.TestCase().assertEqual(stored_chain_map[str(random_chain_id)]["residue_list"],residue_list)
+molecule_name = f"chain_{node_start}_{node_end}"
+molecule_data = pmb.df[
+    (pmb.df["pmb_type"] == "molecule") & (pmb.df["name"] == molecule_name)
+]
+
+# Ensure that exactly one molecule matches
+ut.TestCase().assertEqual(len(molecule_data), 1, f"Molecule {molecule_name} not found in pmb.df")
+
+random_node = random.choice(list(node_topology.keys()))
+node_data = node_topology[random_node]
+node_name = node_data['particle_name']
+node_index = node_data['lattice_index']
+
+stored_node_map = hydrogel_data["node_map"].values[0]
+ut.TestCase().assertEqual(stored_node_map[str(random_node)]["particle_name"], node_name)
+ut.TestCase().assertEqual(all(stored_node_map[str(random_node)]["lattice_index"]), all(node_index))
+
+print("*** Unit Test passed ***")
+
+# Creating hydrogel
+hydrogel_info = pmb.create_hydrogel("my_hydrogel", espresso_system)
+print("*** Hydrogel created: Unit test to verify their name and positions ***")
+
+Node_name_in_espresso = pmb.df[(pmb.df["pmb_type"]=="particle") & (pmb.df["particle_id"]==random_node)]["name"].values[0]
+ut.TestCase().assertEqual(Node_name_in_espresso, node_name)
+ut.TestCase().assertEqual(all(espresso_system.part.by_id(random_node).pos), all(node_index*0.25*diamond_lattice.BOXL))
+
+Chain_name_in_espresso = pmb.df[(pmb.df["pmb_type"]=="molecule") & (pmb.df["molecule_id"]==random_chain_id)]["name"].values[0] 
+Residue_list_in_espresso = pmb.df[(pmb.df["pmb_type"]=="molecule") & (pmb.df["molecule_id"]==random_chain_id)]["residue_list"].values[0]
+Residue_random_name_in_espresso = pmb.df[(pmb.df["pmb_type"]=="residue") & (pmb.df["residue_id"]==random_res_id)]["name"].values[0]
+Residue_random_mol_id = pmb.df[(pmb.df["name"]==Residue_random_name_in_espresso) & (pmb.df["residue_id"]==random_res_id)]["molecule_id"].values[0]
+
+ut.TestCase().assertEqual(Chain_name_in_espresso, "chain_"+node_start+"_"+node_end)
+ut.TestCase().assertEqual(Residue_list_in_espresso, residue_list)
+ut.TestCase().assertEqual(Residue_random_name_in_espresso,  random_res_in_res_list)
+
+vec_between_nodes = (np.array(list(int(x) for x in node_end.strip('[]').split())) -  np.array(list(int(x) for x in node_start.strip('[]').split())))*0.25*diamond_lattice.BOXL
+vec_between_nodes = vec_between_nodes - diamond_lattice.BOXL * np.round(vec_between_nodes/diamond_lattice.BOXL)
+backbone_vector = np.array(list(vec_between_nodes/(diamond_lattice.MPC + 1)))
+random_res_pos_central_bead =  np.array(list(int(x) for x in node_start.strip('[]').split()))*0.25*diamond_lattice.BOXL + (chosen_index + 1)*backbone_vector
+random_res_central_bead_id = hydrogel_info["chains"][Residue_random_mol_id][random_res_id]["central_bead_id"]
+
+np.allclose(espresso_system.part.by_id(random_res_central_bead_id).pos, random_res_pos_central_bead, atol=1e-2)
+print("*** Unit Test passed ***")
+
+print("*** Checking if the ends of the randomly chosen chain is connected to node_start and node_end ***")
+
+molecule_random = hydrogel_info["chains"][random_chain_id]
+numeric_keys = {key: value for key, value in molecule_random.items() if isinstance(key, int)}
+
+# Extract the first and last elements
+keys = list(numeric_keys.keys())
+first_key = keys[0]
+last_key = keys[-1]
+
+Res_node_start = numeric_keys[first_key]
+Res_node_end = numeric_keys[last_key]
+
+central_bead_near_node_start = Res_node_start["central_bead_id"]
+central_bead_near_node_end = Res_node_end["central_bead_id"]
+
+node_ids = [0,1,2,3,4,5,6,7]
+bead_ids_in_random_molecule = [i for i in range(central_bead_near_node_start, central_bead_near_node_end+1)]
+#print(pmb.df[pmb.df["molecule_id"]==random_chain_id])
+filtered_df = pmb.df[
+    pmb.df["particle_id"].isin(node_ids) & 
+    pmb.df["particle_id2"].isin(bead_ids_in_random_molecule)
+    ]
+
+# Extract scalar values for central_bead_node_start and central_bead_node_end
+central_bead_node_start = filtered_df[filtered_df["particle_id2"] == central_bead_near_node_start]["particle_id"].iloc[0]
+central_bead_node_end = filtered_df[filtered_df["particle_id2"] == central_bead_near_node_end]["particle_id"].iloc[0]
+
+bond_name_node_start = filtered_df[
+    (filtered_df["particle_id"] == central_bead_node_start) & 
+    (filtered_df["particle_id2"] == central_bead_near_node_start)
+]["name"].iloc[0]
+
+bond_name_node_end = filtered_df[
+    (filtered_df["particle_id"] == central_bead_node_end) & 
+    (filtered_df["particle_id2"] == central_bead_near_node_end)
+]["name"].iloc[0]
+
+for _, row in filtered_df.iterrows():
+    bond_object = row["bond_object"]
+    if bond_object is None:
+        raise ValueError(f"Bond object is not defined near nodes")
+
+central_bead_name_near_node_start = pmb.df[pmb.df["particle_id"]==central_bead_near_node_start]["name"].values[0]
+central_bead_name_near_node_end = pmb.df[pmb.df["particle_id"]==central_bead_near_node_end]["name"].values[0]
+
+if central_bead_name_near_node_start == BeadType1:
+    possible_bond_names = [NodeType+"-"+BeadType1, BeadType1+"-"+NodeType]
+    assert bond_name_node_start in possible_bond_names
+
+elif central_bead_name_near_node_start == BeadType2:
+    possible_bond_names = [NodeType+"-"+BeadType2, BeadType2+"-"+NodeType]
+    assert bond_name_node_start in possible_bond_names
+
+if central_bead_name_near_node_end == BeadType1:
+    possible_bond_names = [NodeType+"-"+BeadType1, BeadType1+"-"+NodeType]
+    assert bond_name_node_end in possible_bond_names
+
+elif central_bead_name_near_node_end == BeadType2:
+    possible_bond_names = [NodeType+"-"+BeadType2, BeadType2+"-"+NodeType]
+    assert bond_name_node_end in possible_bond_names
+
+print("*** Unit Test passed ***")

testsuite/lattice_builder.py

@@ -125,11 +125,11 @@ def test_lattice_setup(self):
 
         node_positions={}
         node1_label = lattice.node_labels["[1 1 1]"]
-        node_positions[node1_label]=pos_node1
+        node_positions[node1_label]=pos_node1[0]
         node2_label = lattice.node_labels["[0 0 0]"]
-        node_positions[node2_label]=pos_node2
+        node_positions[node2_label]=pos_node2[0]
         node3_label = lattice.node_labels["[2 2 0]"]
-        node_positions[node3_label]=pos_node3
+        node_positions[node3_label]=pos_node3[0]
 
         # define molecule in forward direction
         molecule_name = "chain_[1 1 1]_[0 0 0]"