From 62bca141353cb7d0a4a5dd1bcfa5c93e3ab25ce5 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 17 Nov 2024 23:29:03 +0100 Subject: [PATCH 01/12] hydrogel feature in pyMBE.py adding hydrogel feature pyMBE.py - not sure why it was left out --- lib/lattice.py | 62 ++-- pyMBE.py | 303 +++++++++++++++++- samples/hydrogel.py | 110 +++++++ .../define_and_create_molecules_unit_tests.py | 2 + testsuite/read-write-df_test.py | 6 + 5 files changed, 434 insertions(+), 49 deletions(-) create mode 100644 samples/hydrogel.py diff --git a/lib/lattice.py b/lib/lattice.py index 8b1e942..6fba8be 100644 --- a/lib/lattice.py +++ b/lib/lattice.py @@ -1,4 +1,3 @@ -# # Copyright (C) 2024 pyMBE-dev team # # This file is part of pyMBE. @@ -19,18 +18,18 @@ import itertools import numpy as np +import espressomd +#pmb = pyMBE.pymbe_library(seed=42) 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() `_ @@ -46,11 +45,14 @@ class LatticeBuilder: function to draw the simulation box. """ +# Remove pmb when integrated to main pyMBE.py file + def __init__(self, lattice, strict=True): self.lattice = lattice self.strict = strict self.node_labels = {str([int(x) for x in indices]).replace(",", ""): i for i, indices in enumerate(self.lattice.indices)} + #self.node_labels_arr = {str(indices):i for i, indices in enumerate(self.lattice.indices)} self.nodes = {label: "default_linker" for label in self.node_labels} self.colormap = {} self.chains = {} @@ -58,6 +60,11 @@ def __init__(self, lattice, strict=True): self.kwargs_monomers = {} self.kwargs_bonds = {} self.kwargs_box = {} + #self.pmb = pmb + self.MPC = lattice.MPC + self.BOXL = lattice.BOXL + #self.BOND_LENGTH = BOND_LENGTH.to("reduced_length") if BOND_LENGTH is not None else self.pmb.units.Quantity(1, 'reduced_length') + 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 +100,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 +110,6 @@ def add_default_chains(self, mpc): def draw_lattice(self, ax): """ Draw the lattice in an `Axes3D `_ canvas. - Args: ax: Axes. """ @@ -150,8 +155,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 +170,6 @@ def draw_lattice(self, ax): def draw_simulation_box(self, ax): """ Draw the simulation box in an `Axes3D `_ canvas. - Args: ax: Axes. """ @@ -182,11 +188,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 +206,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 +218,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,38 +230,17 @@ 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. """ @@ -269,17 +248,6 @@ def set_colormap(self, colormap): 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 +262,9 @@ 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_LENGTH): + self.MPC = MPC + self.BOND_LENGTH = BOND_LENGTH + self.BOXL = (self.MPC+1)*self.BOND_LENGTH.magnitude / (np.sqrt(3)*0.25) + diff --git a/pyMBE.py b/pyMBE.py index f31d7fd..aca55bb 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -25,6 +25,7 @@ import pandas as pd import scipy.constants import scipy.optimize +from lib.lattice import LatticeBuilder class pymbe_library(): @@ -53,7 +54,7 @@ def default(self, obj): return obj.item() return super().default(obj) - def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, Kw=None): + def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, Kw=None, lattice_builder=None): """ Initializes the pymbe_library by setting up the reduced unit system with `temperature` and `reduced_length` and sets up the `pmb.df` for bookkeeping. @@ -83,6 +84,7 @@ def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, K Kw=Kw, verbose=False) self.setup_df() + self.lattice_builder = None return def add_bond_in_df(self, particle_id1, particle_id2, use_default_bond=False): @@ -872,8 +874,91 @@ def create_counterions(self, object_name, cation_name, anion_name, espresso_syst for name in counterion_number.keys(): print(f'Ion type: {name} created number: {counterion_number[name]}') return counterion_number + + def create_hydrogel(self, name, espresso_system): + """ + creates the hyrogel `name` in espresso_system - def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, use_default_bond=False): + Args: + name(`str`): Label of the hydrogel to be created. `name` must be defined in the `pmb.df` + espresso_system(`espressomd.system.System`): Instance of a system object from the espressomd library. + + Returns: + + """ + + def format_node(node_list): + return "[" + " ".join(map(str, node_list)) + "]" + + self.check_if_name_is_defined_in_df(name=name, + pmb_type_to_be_defined='hydrogel') + + # placing nodes + + node_positions = {} + node_topology = self.df[self.df["name"]==name]["node_map"].iloc[0] + for index, node_info in node_topology.items(): + node_index = node_info["lattice_index"] + node_name = node_info["particle_name"] + + node_id,node_pos = self.set_node(format_node(node_index), node_name, espresso_system) + node_label = self.lattice_builder.node_labels[format_node(node_index)] + 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_id,chain_info in chain_topology.items(): + 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; + + + def create_hydrogel(self, name, espresso_system): + """ + creates the hyrogel `name` in espresso_system + + Args: + name(`str`): Label of the hydrogel to be created. `name` must be defined in the `pmb.df` + espresso_system(`espressomd.system.System`): Instance of a system object from the espressomd library. + + Returns: + + """ + + def format_node(node_list): + return "[" + " ".join(map(str, node_list)) + "]" + + self.check_if_name_is_defined_in_df(name=name, + pmb_type_to_be_defined='hydrogel') + + # placing nodes + + node_positions = {} + node_topology = self.df[self.df["name"]==name]["node_map"].iloc[0] + for index, node_info in node_topology.items(): + node_index = node_info["lattice_index"] + node_name = node_info["particle_name"] + + node_id,node_pos = self.set_node(format_node(node_index), node_name, espresso_system) + node_label = self.lattice_builder.node_labels[format_node(node_index)] + 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_id,chain_info in chain_topology.items(): + 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; + + def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False): """ Creates `number_of_molecules` molecule of type `name` into `espresso_system` and bookkeeps them into `pmb.df`. @@ -926,10 +1011,12 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi for item in list_of_first_residue_positions: residue_position = [np.array(list_of_first_residue_positions[pos_index])] # Generate an arbitrary random unit vector - backbone_vector = self.generate_random_points_in_a_sphere(center=[0,0,0], - radius=1, + if backbone_vector is None: + backbone_vector = self.generate_random_points_in_a_sphere(center=[0,0,0], + radius=1, n_samples=1, on_surface=True)[0] + residues_info = self.create_residue(name=residue, espresso_system=espresso_system, central_bead_position=residue_position, @@ -1458,6 +1545,82 @@ def define_default_bond(self, bond_type, bond_parameters, epsilon=None, sigma=No new_value = bond_object.get_params(), non_standard_value=True) return + + def define_hydrogel(self, name, node_map, chain_map): + """ + Defines a pyMBE object of type `hydrogel` in `pymbe.df`. + + Args: + name(`str`): Unique label that identifies the `hydrogel`. + node_map(`dict`): {"node_label": {"particle_name": , "lattice_index": }, ... } + chain_map(`dict`): {"chain_id": {"node_start": , "node_end": , "residue_list": , ... } + + """ + + + if self.check_if_name_is_defined_in_df(name=name,pmb_type_to_be_defined='hydrogel'): + return + + index = len(self.df) + self.df.at [index, "name"] = name + self.df.at [index, "pmb_type"] = "hydrogel" + self.add_value_to_df(index = index, + key = ('node_map',''), + new_value = node_map, + non_standard_value=True) + self.add_value_to_df(index = index, + key = ('chain_map',''), + new_value = chain_map, + non_standard_value=True) + + for chain_id in chain_map: + node_start = chain_map[chain_id]["node_start"] + node_end = chain_map[chain_id]["node_end"] + residue_list = chain_map[chain_id]['residue_list'] + # Molecule name + molecule_name = "chain_"+node_start+"_"+node_end + self.define_molecule(name=molecule_name, residue_list=residue_list) + + return; + + + def define_hydrogel(self, name, node_map, chain_map): + """ + Defines a pyMBE object of type `hydrogel` in `pymbe.df`. + + Args: + name(`str`): Unique label that identifies the `hydrogel`. + node_map(`dict`): {"node_label": {"particle_name": , "lattice_index": }, ... } + chain_map(`dict`): {"chain_id": {"node_start": , "node_end": , "residue_list": , ... } + + """ + + + if self.check_if_name_is_defined_in_df(name=name,pmb_type_to_be_defined='hydrogel'): + return + + index = len(self.df) + self.df.at [index, "name"] = name + self.df.at [index, "pmb_type"] = "hydrogel" + self.add_value_to_df(index = index, + key = ('node_map',''), + new_value = node_map, + non_standard_value=True) + self.add_value_to_df(index = index, + key = ('chain_map',''), + new_value = chain_map, + non_standard_value=True) + + for chain_id in chain_map: + node_start = chain_map[chain_id]["node_start"] + node_end = chain_map[chain_id]["node_end"] + residue_list = chain_map[chain_id]['residue_list'] + # Molecule name + molecule_name = "chain_"+node_start+"_"+node_end + self.define_molecule(name=molecule_name, residue_list=residue_list) + + return; + def define_molecule(self, name, residue_list): """ @@ -2210,6 +2373,14 @@ def get_type_map(self): type_map = pd.concat([state_one,state_two],axis=0).to_dict() return type_map + def initialize_lattice_builder(self, diamond_lattice): + """ + Initialize the lattice builder with the DiamondLattice object. + """ + 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 + def load_interaction_parameters(self, filename, verbose=False, overwrite=False): """ Loads the interaction parameters stored in `filename` into `pmb.df` @@ -2600,6 +2771,124 @@ def search_particles_in_residue(self, residue_name): return list_of_particles_in_residue + def set_chain(self, node_start, node_end, node_positions, espresso_system): + """ + 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. + node_positions(`dict`): + espresso_system(`espressomd.system.System`): Instance of a system object from the espressomd library. + """ + if self.lattice_builder is None: + raise ValueError("LatticeBuilder is not initialized. Use `initialize_lattice_builder` first.") + + molecule_name = "chain_"+node_start+"_"+node_end + + sequence = self.df[self.df['name']==molecule_name].residue_list.values [0] + + assert len(sequence) != 0 and not isinstance(sequence, str) + assert len(sequence) == self.lattice_builder.MPC + + key, reverse = self.lattice_builder._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] + + node1 = self.lattice_builder.node_labels[node_start] + node2 = self.lattice_builder.node_labels[node_end] + + # Finding a backbone vector between node_start and node_end + vec_between_nodes = np.array(node_positions[node2]) - np.array(node_positions[node1]) + vec_between_nodes = vec_between_nodes - self.lattice_builder.BOXL * np.round(vec_between_nodes/self.lattice_builder.BOXL) + backbone_vector = list(vec_between_nodes/(self.lattice_builder.MPC + 1)) + chain_molecule_info = self.create_molecule( + name=molecule_name, # Use the name defined earlier + number_of_molecules=1, # Creating one chain + espresso_system=espresso_system, + list_of_first_residue_positions=[list(np.array(node_positions[node1]) + np.array(backbone_vector))], # Start at the first node + 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 = [] + for molecule_id in chain_molecule_info: + for residue_id in chain_molecule_info[molecule_id]: + residue_ids.append(residue_id) + bead_id = chain_molecule_info[molecule_id][residue_id]['central_bead_id'] + chain_ids.append(bead_id) + + self.lattice_builder.chains[key] = sequence + # Search bonds between nodes and chain ends + BeadType_near_to_node_start = self.df[(self.df["residue_id"] == residue_ids[0]) & (self.df["central_bead"].notnull())]["central_bead"].drop_duplicates().iloc[0] + BeadType_near_to_node_end = self.df[(self.df["residue_id"] == residue_ids[-1]) & (self.df["central_bead"].notnull())]["central_bead"].drop_duplicates().iloc[0] + + bond_node1_first_monomer = self.search_bond(particle_name1 = self.lattice_builder.nodes[node_start], + particle_name2 = BeadType_near_to_node_start, + hard_check=False, + use_default_bond=False) + bond_node2_last_monomer = self.search_bond(particle_name1 = self.lattice_builder.nodes[node_end], + particle_name2 = BeadType_near_to_node_end, + hard_check=False, + use_default_bond=False) + + espresso_system.part.by_id(node1).add_bond((bond_node1_first_monomer, chain_ids[0])) + espresso_system.part.by_id(node2).add_bond((bond_node2_last_monomer, chain_ids[-1])) + # Add bonds to data frame + bond_index1 = self.add_bond_in_df(particle_id1=node1, + particle_id2=chain_ids[0], + use_default_bond=False) + self.add_value_to_df(key=('molecule_id',''), + index=int(bond_index1), + new_value=molecule_id, + verbose=False, + overwrite=True) + self.add_value_to_df(key=('residue_id',''), + index=int (bond_index1), + new_value=residue_ids[0], + verbose=False, + overwrite=True) + bond_index2 = self.add_bond_in_df(particle_id1=node2, + particle_id2=chain_ids[-1], + use_default_bond=False) + self.add_value_to_df(key=('molecule_id',''), + index=int(bond_index2), + new_value=molecule_id, + verbose=False, + overwrite=True) + self.add_value_to_df(key=('residue_id',''), + index=int (bond_index2), + new_value=residue_ids[-1], + verbose=False, + overwrite=True) + + def set_node(self, node, residue, espresso_system): + """ + 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. + Returns: + + """ + if self.lattice_builder is None: + 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 + node_particle_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_particle_id,node_position.tolist() + def set_particle_acidity(self, name, acidity=None, default_charge_number=0, pka=None, verbose=True, overwrite=True): """ Sets the particle acidity including the charges in each of its possible states. @@ -3221,7 +3510,11 @@ def setup_df (self): '': object}, 'l0': { '': float}, - } + 'node_map':{ + '':object}, + 'chain_map':{ + '':object}, + } self.df = pd.DataFrame(columns=pd.MultiIndex.from_tuples([(col_main, col_sub) for col_main, sub_cols in columns_dtypes.items() for col_sub in sub_cols.keys()])) diff --git a/samples/hydrogel.py b/samples/hydrogel.py new file mode 100644 index 0000000..50394f3 --- /dev/null +++ b/samples/hydrogel.py @@ -0,0 +1,110 @@ +import pyMBE +import espressomd +import numpy as np +from tqdm import tqdm +import re +from pathlib import Path +import matplotlib.pyplot as plt +import json +import pandas as pd +from lib.handy_functions import minimize_espresso_system_energy +from lib.handy_functions import setup_langevin_dynamics +from lib.analysis import block_analyze +from lib.lattice import DiamondLattice + +# Create an instance of pyMBE library +pmb = pyMBE.pymbe_library(seed=42) +# Creating instance from diamond lattice +reduced_unit_set = pmb.get_reduced_units() +# Monomers per chain +MPC = 8 +# 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 +) + +residue_list = [Res1]*(MPC//2) + [Res2]*(MPC//2) + +# Defining bonds in the hydrogel for all different pairs +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]]) + +# Provide MPC and BOND_LENGTH to Diamond Lattice +diamond_lattice = DiamondLattice(MPC,generic_bond_length) +# Get the box length from the object of DiamondLattice +# and use it to create espresso system +espresso_system = espressomd.System(box_l = [diamond_lattice.BOXL]*3) +pmb.add_bonds_to_espresso(espresso_system = espresso_system) +# Currently pass the pyMBE object into lattice builder +# so you have only one instance of the pyMBE by setting +# self.pmb = pmb in constructor of the LatticeBuilder +lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) # Dont need pmb when integrated to pyMBE.py + +######---creating a hydrogel---########### +# 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()} +# reverse_node_labels looks like {0:"[0 0 0]",1:"[1 1 1]", ..} +connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) for i, j in connectivity} +chain_topology = {} + +def parse_node(node_str): + return list(map(int, node_str.strip("[]").split())) + +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 + +pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology) +node_positions = pmb.create_hydrogel("my_hydrogel", espresso_system) + +from mpl_toolkits.mplot3d import Axes3D +fig = plt.figure() +ax = fig.add_subplot(111,projection="3d") +lattice_builder.draw_lattice(ax) +lattice_builder.draw_simulation_box(ax) +plt.legend(fontsize=12) +plt.show() + + + + diff --git a/testsuite/define_and_create_molecules_unit_tests.py b/testsuite/define_and_create_molecules_unit_tests.py index 8a2c0df..89e7a44 100644 --- a/testsuite/define_and_create_molecules_unit_tests.py +++ b/testsuite/define_and_create_molecules_unit_tests.py @@ -100,6 +100,7 @@ "M2":{"name": "M2", "residue_list": ["R1","R2","R3"]}} + for parameter_set in molecule_parameters.values(): pmb.define_molecule(**parameter_set) @@ -178,6 +179,7 @@ bond = {'r_0' : 0.4*pmb.units.nm, 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + pmb.define_default_bond(bond_type = bond_type, bond_parameters = bond) diff --git a/testsuite/read-write-df_test.py b/testsuite/read-write-df_test.py index d493816..4bccdb7 100644 --- a/testsuite/read-write-df_test.py +++ b/testsuite/read-write-df_test.py @@ -117,6 +117,12 @@ # Read the same pyMBE df from a csv a load it in pyMBE read_df = pmb.read_pmb_df(filename = df_filename) +stored_df['node_map'] = stored_df['node_map'].astype(object) +stored_df['chain_map'] = stored_df['chain_map'].astype(object) + +read_df['node_map'] = read_df['node_map'].astype(object) +read_df['chain_map'] = read_df['chain_map'].astype(object) + # Preprocess data for the Unit Test # The espresso bond object must be converted to a dict in order to compare them using assert_frame_equal stored_df['bond_object'] = stored_df['bond_object'].apply(lambda x: ast.literal_eval(re.subn('HarmonicBond', '', str(x))[0]) if pd.notnull(x) else x) From 20f953b6f9a8c535c6ded6a8e5afbc88b70ee4b5 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Mon, 18 Nov 2024 13:58:27 +0100 Subject: [PATCH 02/12] hydrogels-feature refactored --- lib/lattice.py | 4 -- pyMBE.py | 91 +++------------------------------------------ samples/hydrogel.py | 26 ------------- 3 files changed, 5 insertions(+), 116 deletions(-) diff --git a/lib/lattice.py b/lib/lattice.py index 6fba8be..109e693 100644 --- a/lib/lattice.py +++ b/lib/lattice.py @@ -18,7 +18,6 @@ import itertools import numpy as np -import espressomd #pmb = pyMBE.pymbe_library(seed=42) @@ -52,7 +51,6 @@ 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.node_labels_arr = {str(indices):i for i, indices in enumerate(self.lattice.indices)} self.nodes = {label: "default_linker" for label in self.node_labels} self.colormap = {} self.chains = {} @@ -60,10 +58,8 @@ def __init__(self, lattice, strict=True): self.kwargs_monomers = {} self.kwargs_bonds = {} self.kwargs_box = {} - #self.pmb = pmb self.MPC = lattice.MPC self.BOXL = lattice.BOXL - #self.BOND_LENGTH = BOND_LENGTH.to("reduced_length") if BOND_LENGTH is not None else self.pmb.units.Quantity(1, 'reduced_length') def _get_node_by_label(self, node): diff --git a/pyMBE.py b/pyMBE.py index aca55bb..d75fe2c 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -54,7 +54,7 @@ def default(self, obj): return obj.item() return super().default(obj) - def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, Kw=None, lattice_builder=None): + def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, Kw=None): """ Initializes the pymbe_library by setting up the reduced unit system with `temperature` and `reduced_length` and sets up the `pmb.df` for bookkeeping. @@ -897,11 +897,11 @@ def format_node(node_list): node_positions = {} node_topology = self.df[self.df["name"]==name]["node_map"].iloc[0] - for index, node_info in node_topology.items(): + for node_info in node_topology.values(): node_index = node_info["lattice_index"] node_name = node_info["particle_name"] - node_id,node_pos = self.set_node(format_node(node_index), node_name, espresso_system) + node_pos = self.set_node(format_node(node_index), node_name, espresso_system) node_label = self.lattice_builder.node_labels[format_node(node_index)] node_positions[node_label]=node_pos @@ -909,49 +909,7 @@ def format_node(node_list): # Looping over all the 16 chains chain_topology = self.df[self.df["name"]==name]["chain_map"].iloc[0] - for chain_id,chain_info in chain_topology.items(): - 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; - - - def create_hydrogel(self, name, espresso_system): - """ - creates the hyrogel `name` in espresso_system - - Args: - name(`str`): Label of the hydrogel to be created. `name` must be defined in the `pmb.df` - espresso_system(`espressomd.system.System`): Instance of a system object from the espressomd library. - - Returns: - - """ - - def format_node(node_list): - return "[" + " ".join(map(str, node_list)) + "]" - - self.check_if_name_is_defined_in_df(name=name, - pmb_type_to_be_defined='hydrogel') - - # placing nodes - - node_positions = {} - node_topology = self.df[self.df["name"]==name]["node_map"].iloc[0] - for index, node_info in node_topology.items(): - node_index = node_info["lattice_index"] - node_name = node_info["particle_name"] - - node_id,node_pos = self.set_node(format_node(node_index), node_name, espresso_system) - node_label = self.lattice_builder.node_labels[format_node(node_index)] - 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_id,chain_info in chain_topology.items(): + 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) @@ -1583,45 +1541,6 @@ def define_hydrogel(self, name, node_map, chain_map): return; - - def define_hydrogel(self, name, node_map, chain_map): - """ - Defines a pyMBE object of type `hydrogel` in `pymbe.df`. - - Args: - name(`str`): Unique label that identifies the `hydrogel`. - node_map(`dict`): {"node_label": {"particle_name": , "lattice_index": }, ... } - chain_map(`dict`): {"chain_id": {"node_start": , "node_end": , "residue_list": , ... } - - """ - - - if self.check_if_name_is_defined_in_df(name=name,pmb_type_to_be_defined='hydrogel'): - return - - index = len(self.df) - self.df.at [index, "name"] = name - self.df.at [index, "pmb_type"] = "hydrogel" - self.add_value_to_df(index = index, - key = ('node_map',''), - new_value = node_map, - non_standard_value=True) - self.add_value_to_df(index = index, - key = ('chain_map',''), - new_value = chain_map, - non_standard_value=True) - - for chain_id in chain_map: - node_start = chain_map[chain_id]["node_start"] - node_end = chain_map[chain_id]["node_end"] - residue_list = chain_map[chain_id]['residue_list'] - # Molecule name - molecule_name = "chain_"+node_start+"_"+node_end - self.define_molecule(name=molecule_name, residue_list=residue_list) - - return; - - def define_molecule(self, name, residue_list): """ Defines a pyMBE object of type `molecule` in `pymbe.df`. @@ -2887,7 +2806,7 @@ def set_node(self, node, residue, espresso_system): key = self.lattice_builder._get_node_by_label(node) self.lattice_builder.nodes[key] = residue - return node_particle_id,node_position.tolist() + return node_position.tolist() def set_particle_acidity(self, name, acidity=None, default_charge_number=0, pka=None, verbose=True, overwrite=True): """ diff --git a/samples/hydrogel.py b/samples/hydrogel.py index 50394f3..f99d628 100644 --- a/samples/hydrogel.py +++ b/samples/hydrogel.py @@ -1,20 +1,9 @@ import pyMBE import espressomd -import numpy as np -from tqdm import tqdm -import re -from pathlib import Path import matplotlib.pyplot as plt -import json -import pandas as pd -from lib.handy_functions import minimize_espresso_system_energy -from lib.handy_functions import setup_langevin_dynamics -from lib.analysis import block_analyze from lib.lattice import DiamondLattice -# Create an instance of pyMBE library pmb = pyMBE.pymbe_library(seed=42) -# Creating instance from diamond lattice reduced_unit_set = pmb.get_reduced_units() # Monomers per chain MPC = 8 @@ -58,16 +47,10 @@ # Provide MPC and BOND_LENGTH to Diamond Lattice diamond_lattice = DiamondLattice(MPC,generic_bond_length) -# Get the box length from the object of DiamondLattice -# and use it to create espresso system espresso_system = espressomd.System(box_l = [diamond_lattice.BOXL]*3) pmb.add_bonds_to_espresso(espresso_system = espresso_system) -# Currently pass the pyMBE object into lattice builder -# so you have only one instance of the pyMBE by setting -# self.pmb = pmb in constructor of the LatticeBuilder lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) # Dont need pmb when integrated to pyMBE.py -######---creating a hydrogel---########### # Setting up node topology indices = diamond_lattice.indices node_topology = {} @@ -80,13 +63,9 @@ connectivity = diamond_lattice.connectivity node_labels = lattice_builder.node_labels reverse_node_labels = {v: k for k, v in node_labels.items()} -# reverse_node_labels looks like {0:"[0 0 0]",1:"[1 1 1]", ..} connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) for i, j in connectivity} chain_topology = {} -def parse_node(node_str): - return list(map(int, node_str.strip("[]").split())) - chain_id = 0 for node_s, node_e in connectivity_with_labels: chain_topology[chain_id]={'node_start':node_s, @@ -97,14 +76,9 @@ def parse_node(node_str): pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology) node_positions = pmb.create_hydrogel("my_hydrogel", espresso_system) -from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111,projection="3d") lattice_builder.draw_lattice(ax) lattice_builder.draw_simulation_box(ax) plt.legend(fontsize=12) plt.show() - - - - From abe38610ca7776006eb9a0d8f536db109db38d8e Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Mon, 18 Nov 2024 16:42:56 +0100 Subject: [PATCH 03/12] removed unused variables --- pyMBE.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyMBE.py b/pyMBE.py index d75fe2c..a2edcf6 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -2799,7 +2799,7 @@ 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 - node_particle_id = self.create_particle(name = residue, + self.create_particle(name = residue, espresso_system=espresso_system, number_of_particles=1, position = [node_position]) From 76c384f9fa00ec181805d880a4aebacc8fdb4342 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Wed, 20 Nov 2024 15:31:48 +0100 Subject: [PATCH 04/12] lattice_builder.py - new testsuite WIP --- pyMBE.py | 9 +- testsuite/lattice_builder.py | 214 +++++++++++++++++++++++------------ 2 files changed, 149 insertions(+), 74 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index f8732c8..db127e1 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -2704,7 +2704,11 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system): raise ValueError("LatticeBuilder is not initialized. Use `initialize_lattice_builder` first.") molecule_name = "chain_"+node_start+"_"+node_end - + if molecule_name not in self.df['name'].values: + raise KeyError( + f"Molecule '{molecule_name}' is not defined in the pyMBE DataFrame. " + "Please define it with the correct residue list before setting the chain." + ) sequence = self.df[self.df['name']==molecule_name].residue_list.values [0] assert len(sequence) != 0 and not isinstance(sequence, str) @@ -2720,6 +2724,9 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system): node1 = self.lattice_builder.node_labels[node_start] node2 = self.lattice_builder.node_labels[node_end] + + if node_positions[node1] is None or node_positions[node2] is None: + raise ValueError("Set node position before placing a chain between them") # Finding a backbone vector between node_start and node_end vec_between_nodes = np.array(node_positions[node2]) - np.array(node_positions[node1]) diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index 1481d70..806cbee 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -16,105 +16,173 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # - +import numpy as np import lib.lattice import unittest as ut import matplotlib import matplotlib.pyplot as plt +import pyMBE +import espressomd + matplotlib.use("Agg") # use a non-graphic backend +pmb = pyMBE.pymbe_library(seed=42) +MPC = 4 +BOND_LENGTH = 0.355 * pmb.units.nm + +# Define node particle +NodeType1 = "node_type1" +NodeType2 = "node_type2" + + +# define monomers +BeadType1 = "carbon" +BeadType2 = "oxygen" +BeadType3 = "nitrogen" + +Res1 = "res_1" +Res2 = "res_2" +Res3 = "res_3" + +# Defining bonds in the hydrogel for all different pairs +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} class Test(ut.TestCase): colormap = { - "default_linker": "C0", - "default_monomer": "C1", - "silicon": "C2", - "carbon": "C3", - "oxygen": "C4", - "nitrogen": "C5", - } - - def test_builder(self): - diamond = lib.lattice.DiamondLattice - lattice = lib.lattice.LatticeBuilder(diamond, strict=False) - sequence = ["nitrogen", "carbon", "oxygen", "carbon"] + "default_linker":"green", + "default_monomer":"blue", + Res3: "red", + NodeType2: "orange", + NodeType1: "cyan", + Res1: "yellow", + Res2: "magenta" + } + + @classmethod + def setUpClass(cls): + pmb.define_particle(name=NodeType1, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=NodeType2, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=BeadType1, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=BeadType2, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=BeadType3, sigma=0.355*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_residue( + name=Res1, + central_bead=BeadType1, + side_chains=[] + ) + pmb.define_residue( + name=Res2, + central_bead=BeadType2, + side_chains=[] + ) + pmb.define_residue( + name=Res3, + central_bead=BeadType3, + side_chains=[] + ) + pmb.define_bond(bond_type = 'harmonic', + bond_parameters = HARMONIC_parameters, particle_pairs = [[BeadType1, BeadType1], + [BeadType1, BeadType2], + [BeadType1, BeadType3], + [BeadType2, BeadType2], + [BeadType2, BeadType3], + [BeadType3, BeadType3], + [BeadType1, NodeType1], + [BeadType1, NodeType2], + [BeadType2, NodeType1], + [BeadType2, NodeType2], + [BeadType3, NodeType1], + [BeadType3, NodeType2]]) + + def test_lattice_setup(self): + + diamond = lib.lattice.DiamondLattice(MPC, BOND_LENGTH) + espresso_system = espressomd.System(box_l = [diamond.BOXL]*3) + pmb.add_bonds_to_espresso(espresso_system = espresso_system) + lattice = pmb.initialize_lattice_builder(diamond) + sequence = [Res3, Res1, Res2, Res1] # build default structure - self.assertEqual(len(lattice.nodes), len(diamond.indices)) - self.assertEqual(len(lattice.chains), 0) + assert len(lattice.nodes) == len(diamond.indices) + assert len(lattice.chains) == 0 + + # this function need some work lattice.add_default_chains(mpc=2) - self.assertEqual(len(lattice.chains), len(diamond.connectivity)) + assert len(lattice.chains) == len(diamond.connectivity) + # define custom nodes - self.assertEqual(lattice.get_node("[1 1 1]"), "default_linker") - lattice.set_node(node="[1 1 1]", residue="silicon") - self.assertEqual(lattice.get_node("[1 1 1]"), "silicon") - # define custom chain in forward direction - lattice.set_chain("[0 0 0]", "[1 1 1]", sequence=sequence) - self.assertEqual(lattice.get_chain("[0 0 0]", "[1 1 1]"), sequence) - self.assertEqual(lattice.get_chain("[1 1 1]", "[0 0 0]"), sequence[::-1]) + assert lattice.get_node("[1 1 1]") == "default_linker" + assert lattice.get_node("[0 0 0]") == "default_linker" + + pos_node1 = pmb.set_node("[1 1 1]", NodeType1, espresso_system=espresso_system) + np.testing.assert_equal(actual = lattice.get_node("[1 1 1]"), desired = NodeType1, verbose=True) + pos_node2 = pmb.set_node("[0 0 0]", NodeType2, espresso_system=espresso_system) + np.testing.assert_equal(actual = lattice.get_node("[0 0 0]"), desired = NodeType2, verbose=True) + pos_node3 = pmb.set_node("[2 2 0]", NodeType2, espresso_system=espresso_system) + np.testing.assert_equal(actual = lattice.get_node("[2 2 0]"), desired = NodeType2, verbose=True) + + node_positions={} + node1_label = lattice.node_labels["[1 1 1]"] + node_positions[node1_label]=pos_node1 + node2_label = lattice.node_labels["[0 0 0]"] + node_positions[node2_label]=pos_node2 + node3_label = lattice.node_labels["[2 2 0]"] + node_positions[node3_label]=pos_node3 + + # define molecule in forward direction + molecule_name = "chain_[1 1 1]_[0 0 0]" + pmb.define_molecule(name=molecule_name, residue_list=sequence) + pmb.set_chain("[1 1 1]", "[0 0 0]", node_positions, espresso_system=espresso_system) + np.testing.assert_equal(actual = lattice.get_chain("[1 1 1]", "[0 0 0]"), desired = sequence, verbose=True) + np.testing.assert_equal(actual = lattice.get_chain("[0 0 0]", "[1 1 1]"), desired = sequence[::-1], verbose=True) + # define custom chain in reverse direction - lattice.set_chain("[1 1 1]", "[0 0 0]", sequence=sequence) - self.assertEqual(lattice.get_chain("[0 0 0]", "[1 1 1]"), sequence[::-1]) - self.assertEqual(lattice.get_chain("[1 1 1]", "[0 0 0]"), sequence) + molecule_name = "chain_[0 0 0]_[1 1 1]" + pmb.define_molecule(name=molecule_name, residue_list=sequence) + pmb.set_chain("[0 0 0]", "[1 1 1]", node_positions, espresso_system=espresso_system) + np.testing.assert_equal(lattice.get_chain("[1 1 1]", "[0 0 0]"), sequence[::-1]) + np.testing.assert_equal(lattice.get_chain("[0 0 0]", "[1 1 1]"), sequence) + + ####---Raise Exceptions---#### # define custom chain between normally unconnected nodes - lattice.set_chain("[0 0 0]", "[2 2 0]", sequence=sequence) - self.assertEqual(lattice.get_chain("[0 0 0]", "[2 2 0]"), sequence) - self.assertEqual(lattice.get_chain("[2 2 0]", "[0 0 0]"), sequence[::-1]) + molecule_name = "chain_[0 0 0]_[2 2 0]" + pmb.define_molecule(name=molecule_name, residue_list=sequence) + np.testing.assert_raises(AssertionError, + pmb.set_chain, + "[0 0 0]", "[2 2 0]", node_positions, espresso_system=espresso_system) + # define custom chain that loops - lattice.set_chain("[0 0 0]", "[0 0 0]", sequence=["carbon"]) - self.assertEqual(lattice.get_chain("[0 0 0]", "[0 0 0]"), ["carbon"]) - # colors - lattice.set_colormap(self.colormap) - for index, (label, color) in enumerate(self.colormap.items()): - self.assertEqual(lattice.get_monomer_color(label), color) - self.assertEqual(lattice.get_monomer_color_index(label), index) + molecule_name = "chain_[0 0 0]_[0 0 0]" + pmb.define_molecule(name=molecule_name, residue_list=sequence) + np.testing.assert_raises(AssertionError, + pmb.set_chain, + "[0 0 0]", "[0 0 0]", node_positions, espresso_system=espresso_system) + + lattice.set_colormap(lattice.colormap) + for index, (label, color) in enumerate(lattice.colormap.items()): + np.testing.assert_equal(actual = lattice.get_monomer_color(label),desired = color, verbose=True) + np.testing.assert_equal(actual = lattice.get_monomer_color_index(label),desired = index, verbose=True) - def test_exceptions(self): - diamond = lib.lattice.DiamondLattice - lattice = lib.lattice.LatticeBuilder(diamond, strict=True) + + # Test invalid operations with self.assertRaisesRegex(RuntimeError, "monomer 'unknown' has no associated color in the colormap"): lattice.get_monomer_color("unknown") with self.assertRaises(AssertionError): lattice.set_colormap("red") - with self.assertRaises(AssertionError): - lattice.set_colormap([("oxygen", "red")]) + + # Test node operations with self.assertRaisesRegex(AssertionError, r"node '\[0 5 13\]' doesn't exist in a diamond lattice"): lattice.get_node("[0 5 13]") - with self.assertRaisesRegex(AssertionError, r"node '\[-1 0 0\]' doesn't exist in a diamond lattice"): - lattice.set_node("[-1 0 0]", "carbon") - with self.assertRaisesRegex(RuntimeError, r"no chain has been defined between '\[0 0 0\]' and '\[1 1 1\]' yet"): - lattice.get_chain("[0 0 0]", "[1 1 1]") - with self.assertRaisesRegex(AssertionError, r"node '\[0 5 13\]' doesn't exist in a diamond lattice"): - lattice.get_chain("[0 5 13]", "[1 1 1]") - with self.assertRaisesRegex(AssertionError, r"node '\[0 5 13\]' doesn't exist in a diamond lattice"): - lattice.get_chain("[0 0 0]", "[0 5 13]") - with self.assertRaisesRegex(AssertionError, r"there is no chain between '\[0 0 0\]' and '\[2 2 0\]' in a diamond lattice \(strict mode is enabled\)"): - lattice.get_chain("[0 0 0]", "[2 2 0]") - with self.assertRaisesRegex(AssertionError, r"there is no chain between '\[0 0 0\]' and '\[0 0 0\]' in a diamond lattice \(strict mode is enabled\)"): - lattice.get_chain("[0 0 0]", "[0 0 0]") - with self.assertRaises(AssertionError): - lattice.set_chain("[0 0 0]", "[1 1 1]", sequence=[]) - with self.assertRaises(AssertionError): - lattice.set_chain("[0 0 0]", "[1 1 1]", sequence="oxygen") - with self.assertRaisesRegex(AssertionError, r"chain cannot be defined between '\[0 0 0\]' and '\[0 0 0\]' since it would form a loop with a non-symmetric sequence \(under-defined stereocenter\)"): - lattice.strict = False - lattice.set_chain("[0 0 0]", "[0 0 0]", sequence=["carbon", "oxygen"]) - - def test_plot(self): - # smoke test: check for runtime errors when drawing on a 3D canvas - sequence = ["nitrogen", "carbon", "oxygen", "carbon"] - diamond = lib.lattice.DiamondLattice - lattice = lib.lattice.LatticeBuilder(diamond) - lattice.add_default_chains(mpc=2) - lattice.set_node(node="[1 1 1]", residue="silicon") - lattice.set_chain("[0 0 0]", "[1 1 1]", sequence=sequence) - # with default colormap + + # Test plot fig = plt.figure(figsize=(12, 12)) ax = fig.add_subplot(projection="3d", computed_zorder=False) lattice.draw_lattice(ax) lattice.draw_simulation_box(ax) - ax.legend() plt.close(fig) - # with custom colormap + fig = plt.figure(figsize=(12, 12)) ax = fig.add_subplot(projection="3d", computed_zorder=False) lattice.set_colormap(self.colormap) @@ -123,6 +191,6 @@ def test_plot(self): ax.legend() plt.close(fig) - if __name__ == "__main__": ut.main() + From 3e84660080b128a2e9fa810b7468b8fdc3e4a908 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 1 Dec 2024 20:35:35 +0100 Subject: [PATCH 05/12] lattice_builder modified --- testsuite/lattice_builder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index 806cbee..8dbde8b 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -160,8 +160,8 @@ def test_lattice_setup(self): pmb.set_chain, "[0 0 0]", "[0 0 0]", node_positions, espresso_system=espresso_system) - lattice.set_colormap(lattice.colormap) - for index, (label, color) in enumerate(lattice.colormap.items()): + lattice.set_colormap(self.colormap) + for index, (label, color) in enumerate(self.colormap.items()): np.testing.assert_equal(actual = lattice.get_monomer_color(label),desired = color, verbose=True) np.testing.assert_equal(actual = lattice.get_monomer_color_index(label),desired = index, verbose=True) From d04628e9c3292fc4027b9ce432d0ff38f3fe1b2e Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Fri, 20 Dec 2024 13:59:42 +0100 Subject: [PATCH 06/12] new testuites for hydrogel builder --- lib/lattice.py | 2 + pyMBE.py | 42 +++++-- samples/hydrogel.py | 2 + testsuite/hydrogel_builder.py | 225 ++++++++++++++++++++++++++++++++++ testsuite/lattice_builder.py | 6 +- 5 files changed, 264 insertions(+), 13 deletions(-) create mode 100644 testsuite/hydrogel_builder.py diff --git a/lib/lattice.py b/lib/lattice.py index 109e693..7e97690 100644 --- a/lib/lattice.py +++ b/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) diff --git a/pyMBE.py b/pyMBE.py index db127e1..dd90fc4 100644 --- a/pyMBE.py +++ b/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): """ diff --git a/samples/hydrogel.py b/samples/hydrogel.py index f99d628..32fe669 100644 --- a/samples/hydrogel.py +++ b/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) diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py new file mode 100644 index 0000000..3160aaa --- /dev/null +++ b/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 ***") diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index 8dbde8b..81aa28c 100644 --- a/testsuite/lattice_builder.py +++ b/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]" From 97f932acfb18fb168374894d44cf577ebb2cfe0b Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Fri, 20 Dec 2024 15:05:45 +0100 Subject: [PATCH 07/12] Fix R installation issue in testsuite.yml --- .github/workflows/testsuite.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/testsuite.yml b/.github/workflows/testsuite.yml index 3c408d5..1705acb 100644 --- a/.github/workflows/testsuite.yml +++ b/.github/workflows/testsuite.yml @@ -19,6 +19,8 @@ jobs: uses: actions/checkout@main - name: Validate CITATION.cff uses: dieghernan/cff-validator@v3 + with: + install-r: true # Ensure R is explicitly installed - name: Setup EESSI uses: eessi/github-action-eessi@v3 with: From 17188f72396f3e4c559396acf3e0621acef842e1 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Fri, 20 Dec 2024 15:15:46 +0100 Subject: [PATCH 08/12] removed unused variables in hydrogel_builder.py --- testsuite/hydrogel_builder.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py index 3160aaa..b50b0a7 100644 --- a/testsuite/hydrogel_builder.py +++ b/testsuite/hydrogel_builder.py @@ -1,9 +1,6 @@ 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 @@ -201,7 +198,7 @@ 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") + raise ValueError("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] From 23cc6907927912ec36c2cb43f4ac034c18d136d3 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Thu, 26 Dec 2024 18:51:22 +0100 Subject: [PATCH 09/12] refactoring hydrogel builder testsuite for code coverage --- pyMBE.py | 35 ++-- testsuite/hydrogel_builder.py | 325 ++++++++++++++++++++++++++++------ 2 files changed, 295 insertions(+), 65 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index dd90fc4..3950a8e 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -885,12 +885,9 @@ def create_hydrogel(self, name, espresso_system): Returns: """ - def format_node(node_list): - return "[" + " ".join(map(str, node_list)) + "]" - - self.check_if_name_is_defined_in_df(name=name, - pmb_type_to_be_defined='hydrogel') - + if not self.check_if_name_is_defined_in_df(name=name, pmb_type_to_be_defined='hydrogel'): + raise ValueError(f"Hydrogel with name '{name}' is not defined in the DataFrame.") + hydrogel_info={"name":name, "chains":{}, "nodes":{}} # placing nodes node_positions = {} @@ -898,9 +895,9 @@ def format_node(node_list): for node_info in node_topology.values(): node_index = node_info["lattice_index"] node_name = node_info["particle_name"] - 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_pos, node_id = self.set_node(self.format_node(node_index), node_name, espresso_system) + node_label = self.lattice_builder.node_labels[self.format_node(node_index)] + hydrogel_info["nodes"][self.format_node(node_index)]=node_id node_positions[node_label]=node_pos # Placing chains between nodes @@ -1550,15 +1547,14 @@ def define_hydrogel(self, name, node_map, chain_map): key = ('chain_map',''), new_value = chain_map, non_standard_value=True) - for chain_id in chain_map: node_start = chain_map[chain_id]["node_start"] node_end = chain_map[chain_id]["node_end"] residue_list = chain_map[chain_id]['residue_list'] # Molecule name molecule_name = "chain_"+node_start+"_"+node_end + self.delete_molecule_entry(molecule_name) self.define_molecule(name=molecule_name, residue_list=residue_list) - return; def define_molecule(self, name, residue_list): @@ -1770,7 +1766,17 @@ def define_residue(self, name, central_bead, side_chains): self.df.at [index,'pmb_type'] = 'residue' self.df.at [index,'central_bead'] = central_bead self.df.at [index,('side_chains','')] = side_chains - return + return + + def delete_molecule_entry(self, molecule_name): + """ + Deletes a molecule entry from the DataFrame if it exists. + + Args: + molecule_name (`str`): The name of the molecule to delete. + """ + if molecule_name in self.df["name"].values: + self.df = self.df[self.df["name"] != molecule_name].reset_index(drop=True) def destroy_pmb_object_in_system(self, name, espresso_system): """ @@ -1991,6 +1997,10 @@ def find_value_from_es_type(self, es_type, column_name): return column_name_value return None + def format_node(self, node_list): + return "[" + " ".join(map(str, node_list)) + "]" + + def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_samples): """ Generates coordinates outside a sphere centered at `center`. @@ -2731,7 +2741,6 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system): "Please define it with the correct residue list before setting the chain." ) sequence = self.df[self.df['name']==molecule_name].residue_list.values [0] - assert len(sequence) != 0 and not isinstance(sequence, str) assert len(sequence) == self.lattice_builder.MPC diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py index b50b0a7..d035e19 100644 --- a/testsuite/hydrogel_builder.py +++ b/testsuite/hydrogel_builder.py @@ -4,6 +4,8 @@ import pyMBE from lib.lattice import DiamondLattice import espressomd +from unittest.mock import patch +import io pmb = pyMBE.pymbe_library(seed=42) MPC = 0 @@ -53,6 +55,12 @@ with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): diamond_lattice = DiamondLattice(MPC, generic_bond_length) +with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): + diamond_lattice = DiamondLattice("invalid", generic_bond_length) # MPC is a string + +with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): + diamond_lattice = DiamondLattice(-5, generic_bond_length) # MPC is negative + print("*** Unit Test passed ***") MPC=8 @@ -83,51 +91,159 @@ '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))] +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))] +################################################ +incomplete_node_map = { + 0: {"particle_name": NodeType, "lattice_index": [0, 0, 0]}, + 1: {"particle_name": NodeType, "lattice_index": [1, 1, 1]}, + } +incomplete_chain_map = {0: {"node_start": "[0 0 0]", "node_end":"[1 1 1]" , "residue_list": residue_list}} + +def test_incomplete_node_map(): + #node_map = {0: {"particle_name": "node0", "lattice_index": [0, 0, 0]}} + #chain_map = {i: {"node_start": str(i), "node_end": str(i+1), "residue_list": ["A", "B", "C"]} for i in range(16)} + + with patch('sys.stdout', new=io.StringIO()) as mock_stdout: + pmb.define_hydrogel(name="test_hydrogel", node_map = incomplete_node_map, chain_map = chain_topology) + output = mock_stdout.getvalue().strip() + + expected_output = ( + "Incomplete hydrogel: node_map must contain exactly 8 node_labels (0 through 7). " + "Found: [0, 1]" + ) + assert output == expected_output, f"Unexpected print output: {output}" + + +def test_incomplete_chain_map(): + #node_map = {i: {"particle_name": f"node{i}", "lattice_index": [i, i, i]} for i in range(8)} + #chain_map = {i: {"node_start": str(i), "node_end": str(i+1), "residue_list": ["A", "B", "C"]} for i in range(10)} + + with patch('sys.stdout', new=io.StringIO()) as mock_stdout: + pmb.define_hydrogel(name="test_hydrogel", node_map = node_topology, chain_map = incomplete_chain_map) + output = mock_stdout.getvalue().strip() + + expected_output = ( + "Incomplete hydrogel: chain_map must contain exactly 16 chain_ids (0 through 15). " + "Found: [0]" + ) + assert output == expected_output, f"Unexpected print output: {output}" +test_incomplete_node_map() +test_incomplete_chain_map() +################################################################## +#pmb.define_hydrogel("test_hydrogel", incomplete_node_map, chain_map) +#pmb.create_hydrogel("test_hydrogel", espresso_system) +def compare_node_maps(map1, map2): + # Ensure keys are comparable + map1 = {int(k): v for k, v in map1.items()} + map2 = {int(k): v for k, v in map2.items()} + + # Compare key sets + if set(map1.keys()) != set(map2.keys()): + return False + + # Compare each node's details + for key in map1: + if map1[key]["particle_name"] != map2[key]["particle_name"]: + return False + if not np.array_equal(np.array(map1[key]["lattice_index"]), np.array(map2[key]["lattice_index"])): + return False + + return True + +def parse_string_to_array(string): + """ + Convert a string representation of a list (e.g., '[3 1 3]') into a numpy array. + """ + string = string.strip("[]") # Remove brackets + elements = map(int, string.split()) # Split by spaces and convert to integers + return np.array(list(elements)) + +def compare_chain_maps(map1, map2): + # Normalize keys to integers + map1 = {int(k): v for k, v in map1.items()} + map2 = {int(k): v for k, v in map2.items()} + + # Compare key sets + if set(map1.keys()) != set(map2.keys()): + return False + + # Compare each chain's details + for key in map1: + # Compare node_start and node_end + node_start_1 = parse_string_to_array(map1[key]["node_start"]) + node_start_2 = parse_string_to_array(map2[key]["node_start"]) + if not np.array_equal(node_start_1, node_start_2): + return False + + node_end_1 = parse_string_to_array(map1[key]["node_end"]) + node_end_2 = parse_string_to_array(map2[key]["node_end"]) + if not np.array_equal(node_end_1, node_end_2): + return False + + # Compare residue_list + if map1[key]["residue_list"] != map2[key]["residue_list"]: + return False + + return True +####################################################### pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology) +#def test_hydrogel_already_defined(): +existing_hydrogel_name = "my_hydrogel" +pmb.define_hydrogel(existing_hydrogel_name,node_topology, chain_topology) +hydrogel_count = len(pmb.df[pmb.df["name"] == existing_hydrogel_name]) +assert hydrogel_count == 1, f"Hydrogel '{existing_hydrogel_name}' should not be redefined." +#test_hydrogel_already_defined() +assert existing_hydrogel_name in pmb.df["name"].values +assert pmb.df.loc[pmb.df["name"] == existing_hydrogel_name, "pmb_type"].values[0] == "hydrogel" + +# Verify node_map and chain_map are correctly added +assert compare_node_maps(pmb.df.loc[pmb.df["name"] == existing_hydrogel_name, "node_map"].values[0], node_topology) +assert compare_chain_maps(pmb.df.loc[pmb.df["name"] == existing_hydrogel_name, "chain_map"].values[0], chain_topology) +for chain_id in chain_topology: + molecule_name = f"chain_{chain_topology[chain_id]['node_start']}_{chain_topology[chain_id]['node_end']}" + assert molecule_name in pmb.df["name"].values 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 +#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'] +#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) +#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 +#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] +#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) -] +#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") +#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'] +#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)) +#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 ***") @@ -135,30 +251,135 @@ 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 ***") +############################ +def test_format_node(): + assert pmb.format_node([1, 2, 3]) == "[1 2 3]" + assert pmb.format_node([4, 5, 6]) == "[4 5 6]" +test_format_node() +############################ + +def test_hydrogel_info(): + assert hydrogel_info["name"] == "my_hydrogel" + +test_hydrogel_info() +########################### +def test_node_positions(): + for node_index, node_id in hydrogel_info["nodes"].items(): + # Ensure that the node is added to the espresso system and has a position + node_pos = espresso_system.part.by_id(int(node_id[0])).pos + node_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "particle") & (pmb.df["particle_id"] == node_id[0])]["name"].values[0] + node_data = node_topology[node_id[0]] + node_name = node_data["particle_name"] + + # Assert node's name and position are correctly set + ut.TestCase().assertEqual(node_name_in_espresso, node_name) + ut.TestCase().assertTrue(np.allclose(node_pos, np.array(node_data["lattice_index"]) * 0.25 * diamond_lattice.BOXL, atol=1e-7)) +test_node_positions() +############################# +def test_chain_placement_and_connectivity(): + for molecule_id, molecule_data in hydrogel_info["chains"].items(): + # Ensure that chain's node_start and node_end are correctly set + node_start = molecule_data["node_start"] + node_end = molecule_data["node_end"] + chain_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "molecule") & (pmb.df["molecule_id"] == molecule_id)]["name"].values[0] + + # Assert chain's node_start and node_end + ut.TestCase().assertEqual(chain_name_in_espresso, f"chain_{node_start}_{node_end}") + + # Check if chain is connected in the espresso system (e.g., check bond or distance between node_start and node_end) + node_start_id = hydrogel_info["nodes"][node_start][0] + node_end_id = hydrogel_info["nodes"][node_end][0] + start_pos = espresso_system.part.by_id(int(node_start_id)).pos + end_pos = espresso_system.part.by_id(int(node_end_id)).pos + vec_between_nodes = end_pos - start_pos + # Ensure that the chain is connected (check distance, should be within acceptable bond length range) + vec_between_nodes = vec_between_nodes - diamond_lattice.BOXL * np.round(vec_between_nodes / diamond_lattice.BOXL) + distance_between_nodes = np.linalg.norm(vec_between_nodes) + ut.TestCase().assertAlmostEqual(distance_between_nodes, (diamond_lattice.MPC+1)*generic_bond_length.magnitude, delta=0.0000001) +test_chain_placement_and_connectivity() +######################################## +def test_all_residue_placement(): + for chain_id, chain_data in hydrogel_info["chains"].items(): + node_start = chain_data["node_start"] + node_end = chain_data["node_end"] + + # Get the list of residues (keys in chain_data, excluding node_start/node_end) + residues_in_chain = {k: v for k, v in chain_data.items() if isinstance(k, int)} + + # Loop through each residue in the chain + for (res_id, res_data) in residues_in_chain.items(): + central_bead_id = res_data["central_bead_id"] + + # Get the position of the central bead from the espresso system + central_bead_pos = espresso_system.part.by_id(central_bead_id).pos + + # Calculate the expected position of the residue's central bead + residue_index = list(residues_in_chain.keys()).index(res_id) + + vec_between_nodes = (np.array([float(x) for x in node_end.strip('[]').split()]) - + np.array([float(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 = vec_between_nodes / (diamond_lattice.MPC + 1) + expected_position = (np.array([float(x) for x in node_start.strip('[]').split()]) * 0.25 * diamond_lattice.BOXL + + (residue_index + 1) * backbone_vector) + + # Validate that the central bead's position matches the expected position + np.allclose(central_bead_pos, expected_position, atol=1e-2) + expected_res_name = chain_topology[chain_id]["residue_list"][residue_index] + expected_node_start = chain_topology[chain_id]["node_start"] + expected_node_end = chain_topology[chain_id]["node_end"] + residue_name = pmb.df[(pmb.df["pmb_type"]=="residue") & (pmb.df["residue_id"]==res_id)]["name"].values[0] + ut.TestCase().assertEqual(node_start, expected_node_start) + ut.TestCase().assertEqual(node_end, expected_node_end) + ut.TestCase().assertEqual(residue_name, expected_res_name) +test_all_residue_placement() +############################################################################################################################### + +#####-- Invalid hydrogel name --##### +# Test if create_hydrogel raises an exception when provided with invalid data +print("*** Unit Test: Check invalid inputs for create_hydrogel ***") + +with ut.TestCase().assertRaises(ValueError): + pmb.create_hydrogel("invalid_hydrogel", espresso_system) # No such hydrogel in df + +print("*** Invalid Input Test Passed ***") +####################################################################### + +# Check if the molecules (chains) are correctly stored in the hydrogel data +for ((molecule_id, molecule_data),chain_id_in_input_dict) in zip(hydrogel_info["chains"].items(),chain_topology.items()): + molecule_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "molecule") & (pmb.df["molecule_id"] == molecule_id)]["name"].values[0] + ut.TestCase().assertEqual(molecule_name_in_espresso, f"chain_{molecule_data['node_start']}_{molecule_data['node_end']}") + Residue_list_in_espresso = pmb.df[(pmb.df["pmb_type"]=="molecule") & (pmb.df["molecule_id"]==molecule_id)]["residue_list"].values[0] + residue_list = chain_id_in_input_dict[1]["residue_list"] + ut.TestCase().assertEqual(Residue_list_in_espresso, residue_list) + +######################################################################## + + +#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 ***") - +random_chain_id = random.choice(list(chain_topology.keys())) # Choose a random chain of 0-15 molecule_random = hydrogel_info["chains"][random_chain_id] numeric_keys = {key: value for key, value in molecule_random.items() if isinstance(key, int)} From 660451a67d85f981ceb8e5450ed6983a972caac0 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Thu, 26 Dec 2024 19:13:34 +0100 Subject: [PATCH 10/12] refactoring hydrogel builder testsuite for code coverage --- testsuite/hydrogel_builder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py index d035e19..4a3421e 100644 --- a/testsuite/hydrogel_builder.py +++ b/testsuite/hydrogel_builder.py @@ -264,7 +264,7 @@ def test_hydrogel_info(): test_hydrogel_info() ########################### def test_node_positions(): - for node_index, node_id in hydrogel_info["nodes"].items(): + for _, node_id in hydrogel_info["nodes"].items(): # Ensure that the node is added to the espresso system and has a position node_pos = espresso_system.part.by_id(int(node_id[0])).pos node_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "particle") & (pmb.df["particle_id"] == node_id[0])]["name"].values[0] From d9883d85340def11bcf251d66df355da89b7961c Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Mon, 6 Jan 2025 19:56:09 +0100 Subject: [PATCH 11/12] refactoringhydrogel builder testsuite for code coverage --- pyMBE.py | 17 +- testsuite/hydrogel_builder.py | 302 +++++++++++----------------------- testsuite/lattice_builder.py | 8 + 3 files changed, 108 insertions(+), 219 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index 3950a8e..8971fc3 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -1519,19 +1519,17 @@ def define_hydrogel(self, name, node_map, chain_map): 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( + raise ValueError( f"Incomplete hydrogel: node_map must contain exactly 8 node_labels (0 through 7). " - f"Found: {sorted(actual_node_labels)}" - ) + 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( + raise ValueError( f"Incomplete hydrogel: chain_map must contain exactly 16 chain_ids (0 through 15). " - f"Found: {sorted(actual_chain_ids)}" - ) + f"Found: {sorted(actual_chain_ids)}") if self.check_if_name_is_defined_in_df(name=name,pmb_type_to_be_defined='hydrogel'): return @@ -2735,11 +2733,6 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system): raise ValueError("LatticeBuilder is not initialized. Use `initialize_lattice_builder` first.") molecule_name = "chain_"+node_start+"_"+node_end - if molecule_name not in self.df['name'].values: - raise KeyError( - f"Molecule '{molecule_name}' is not defined in the pyMBE DataFrame. " - "Please define it with the correct residue list before setting the chain." - ) sequence = self.df[self.df['name']==molecule_name].residue_list.values [0] assert len(sequence) != 0 and not isinstance(sequence, str) assert len(sequence) == self.lattice_builder.MPC @@ -2755,7 +2748,7 @@ def set_chain(self, node_start, node_end, node_positions, espresso_system): node1 = self.lattice_builder.node_labels[node_start] node2 = self.lattice_builder.node_labels[node_end] - if node_positions[node1] is None or node_positions[node2] is None: + if not node1 in node_positions or not node2 in node_positions: raise ValueError("Set node position before placing a chain between them") # Finding a backbone vector between node_start and node_end diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py index 4a3421e..e244760 100644 --- a/testsuite/hydrogel_builder.py +++ b/testsuite/hydrogel_builder.py @@ -4,19 +4,13 @@ import pyMBE from lib.lattice import DiamondLattice import espressomd -from unittest.mock import patch -import io 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')) @@ -36,7 +30,7 @@ central_bead=BeadType2, # Define the central bead name side_chains=[] # Assuming no side chains for the monomer ) - +# define bond parameters 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, @@ -51,17 +45,18 @@ [NodeType, BeadType2]]) print("*** Unit Test: check that only non-negative values of monomers per chain are allowed ***") +MPC=0 +np.testing.assert_raises(ValueError, DiamondLattice, MPC, generic_bond_length) +MPC="invalid" +np.testing.assert_raises(ValueError, DiamondLattice, MPC, generic_bond_length) +MPC=-5 +np.testing.assert_raises(ValueError, DiamondLattice, MPC, generic_bond_length) +print("*** Unit Test passed ***") -with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): - diamond_lattice = DiamondLattice(MPC, generic_bond_length) - -with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): - diamond_lattice = DiamondLattice("invalid", generic_bond_length) # MPC is a string - -with ut.TestCase().assertRaisesRegex(ValueError, "MPC must be a non-zero positive integer."): - diamond_lattice = DiamondLattice(-5, generic_bond_length) # MPC is negative +print("*** Unit test: check that any objects are other than DiamondLattice passed to initialize_lattice_builder raises a TypeError ***") +np.testing.assert_raises(TypeError, pmb.initialize_lattice_builder, None) +print("*** Unit test passed ***") -print("*** Unit Test passed ***") MPC=8 diamond_lattice = DiamondLattice(MPC, generic_bond_length) @@ -103,39 +98,10 @@ } incomplete_chain_map = {0: {"node_start": "[0 0 0]", "node_end":"[1 1 1]" , "residue_list": residue_list}} -def test_incomplete_node_map(): - #node_map = {0: {"particle_name": "node0", "lattice_index": [0, 0, 0]}} - #chain_map = {i: {"node_start": str(i), "node_end": str(i+1), "residue_list": ["A", "B", "C"]} for i in range(16)} - - with patch('sys.stdout', new=io.StringIO()) as mock_stdout: - pmb.define_hydrogel(name="test_hydrogel", node_map = incomplete_node_map, chain_map = chain_topology) - output = mock_stdout.getvalue().strip() - - expected_output = ( - "Incomplete hydrogel: node_map must contain exactly 8 node_labels (0 through 7). " - "Found: [0, 1]" - ) - assert output == expected_output, f"Unexpected print output: {output}" - - -def test_incomplete_chain_map(): - #node_map = {i: {"particle_name": f"node{i}", "lattice_index": [i, i, i]} for i in range(8)} - #chain_map = {i: {"node_start": str(i), "node_end": str(i+1), "residue_list": ["A", "B", "C"]} for i in range(10)} - - with patch('sys.stdout', new=io.StringIO()) as mock_stdout: - pmb.define_hydrogel(name="test_hydrogel", node_map = node_topology, chain_map = incomplete_chain_map) - output = mock_stdout.getvalue().strip() - - expected_output = ( - "Incomplete hydrogel: chain_map must contain exactly 16 chain_ids (0 through 15). " - "Found: [0]" - ) - assert output == expected_output, f"Unexpected print output: {output}" -test_incomplete_node_map() -test_incomplete_chain_map() -################################################################## -#pmb.define_hydrogel("test_hydrogel", incomplete_node_map, chain_map) -#pmb.create_hydrogel("test_hydrogel", espresso_system) +np.testing.assert_raises(ValueError, pmb.define_hydrogel, "test_hydrogel", incomplete_node_map, chain_topology) + +np.testing.assert_raises(ValueError, pmb.define_hydrogel, "test_hydrogel", node_topology, incomplete_chain_map) + def compare_node_maps(map1, map2): # Ensure keys are comparable map1 = {int(k): v for k, v in map1.items()} @@ -151,7 +117,6 @@ def compare_node_maps(map1, map2): return False if not np.array_equal(np.array(map1[key]["lattice_index"]), np.array(map2[key]["lattice_index"])): return False - return True def parse_string_to_array(string): @@ -187,16 +152,14 @@ def compare_chain_maps(map1, map2): # Compare residue_list if map1[key]["residue_list"] != map2[key]["residue_list"]: return False - return True + ####################################################### pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology) -#def test_hydrogel_already_defined(): existing_hydrogel_name = "my_hydrogel" pmb.define_hydrogel(existing_hydrogel_name,node_topology, chain_topology) hydrogel_count = len(pmb.df[pmb.df["name"] == existing_hydrogel_name]) assert hydrogel_count == 1, f"Hydrogel '{existing_hydrogel_name}' should not be redefined." -#test_hydrogel_already_defined() assert existing_hydrogel_name in pmb.df["name"].values assert pmb.df.loc[pmb.df["name"] == existing_hydrogel_name, "pmb_type"].values[0] == "hydrogel" @@ -207,176 +170,101 @@ def compare_chain_maps(map1, map2): molecule_name = f"chain_{chain_topology[chain_id]['node_start']}_{chain_topology[chain_id]['node_end']}" assert molecule_name in pmb.df["name"].values -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 ***") -############################ -def test_format_node(): - assert pmb.format_node([1, 2, 3]) == "[1 2 3]" - assert pmb.format_node([4, 5, 6]) == "[4 5 6]" -test_format_node() ############################ -def test_hydrogel_info(): - assert hydrogel_info["name"] == "my_hydrogel" - -test_hydrogel_info() -########################### -def test_node_positions(): - for _, node_id in hydrogel_info["nodes"].items(): - # Ensure that the node is added to the espresso system and has a position - node_pos = espresso_system.part.by_id(int(node_id[0])).pos - node_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "particle") & (pmb.df["particle_id"] == node_id[0])]["name"].values[0] - node_data = node_topology[node_id[0]] - node_name = node_data["particle_name"] - - # Assert node's name and position are correctly set - ut.TestCase().assertEqual(node_name_in_espresso, node_name) - ut.TestCase().assertTrue(np.allclose(node_pos, np.array(node_data["lattice_index"]) * 0.25 * diamond_lattice.BOXL, atol=1e-7)) -test_node_positions() -############################# -def test_chain_placement_and_connectivity(): - for molecule_id, molecule_data in hydrogel_info["chains"].items(): - # Ensure that chain's node_start and node_end are correctly set - node_start = molecule_data["node_start"] - node_end = molecule_data["node_end"] - chain_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "molecule") & (pmb.df["molecule_id"] == molecule_id)]["name"].values[0] - - # Assert chain's node_start and node_end - ut.TestCase().assertEqual(chain_name_in_espresso, f"chain_{node_start}_{node_end}") - - # Check if chain is connected in the espresso system (e.g., check bond or distance between node_start and node_end) - node_start_id = hydrogel_info["nodes"][node_start][0] - node_end_id = hydrogel_info["nodes"][node_end][0] - start_pos = espresso_system.part.by_id(int(node_start_id)).pos - end_pos = espresso_system.part.by_id(int(node_end_id)).pos - vec_between_nodes = end_pos - start_pos - # Ensure that the chain is connected (check distance, should be within acceptable bond length range) - vec_between_nodes = vec_between_nodes - diamond_lattice.BOXL * np.round(vec_between_nodes / diamond_lattice.BOXL) - distance_between_nodes = np.linalg.norm(vec_between_nodes) - ut.TestCase().assertAlmostEqual(distance_between_nodes, (diamond_lattice.MPC+1)*generic_bond_length.magnitude, delta=0.0000001) -test_chain_placement_and_connectivity() -######################################## -def test_all_residue_placement(): - for chain_id, chain_data in hydrogel_info["chains"].items(): - node_start = chain_data["node_start"] - node_end = chain_data["node_end"] - - # Get the list of residues (keys in chain_data, excluding node_start/node_end) - residues_in_chain = {k: v for k, v in chain_data.items() if isinstance(k, int)} - - # Loop through each residue in the chain - for (res_id, res_data) in residues_in_chain.items(): - central_bead_id = res_data["central_bead_id"] - - # Get the position of the central bead from the espresso system - central_bead_pos = espresso_system.part.by_id(central_bead_id).pos - - # Calculate the expected position of the residue's central bead - residue_index = list(residues_in_chain.keys()).index(res_id) - - vec_between_nodes = (np.array([float(x) for x in node_end.strip('[]').split()]) - - np.array([float(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 = vec_between_nodes / (diamond_lattice.MPC + 1) - expected_position = (np.array([float(x) for x in node_start.strip('[]').split()]) * 0.25 * diamond_lattice.BOXL + - (residue_index + 1) * backbone_vector) - - # Validate that the central bead's position matches the expected position - np.allclose(central_bead_pos, expected_position, atol=1e-2) - expected_res_name = chain_topology[chain_id]["residue_list"][residue_index] - expected_node_start = chain_topology[chain_id]["node_start"] - expected_node_end = chain_topology[chain_id]["node_end"] - residue_name = pmb.df[(pmb.df["pmb_type"]=="residue") & (pmb.df["residue_id"]==res_id)]["name"].values[0] - ut.TestCase().assertEqual(node_start, expected_node_start) - ut.TestCase().assertEqual(node_end, expected_node_end) - ut.TestCase().assertEqual(residue_name, expected_res_name) -test_all_residue_placement() -############################################################################################################################### +class Test(ut.TestCase): + + def test_format_node(self): + assert pmb.format_node([1, 2, 3]) == "[1 2 3]" + assert pmb.format_node([4, 5, 6]) == "[4 5 6]" + + def test_hydrogel_info(self): + assert hydrogel_info["name"] == "my_hydrogel" + + def test_node_positions(self): + for _, node_id in hydrogel_info["nodes"].items(): + node_pos = espresso_system.part.by_id(int(node_id[0])).pos + node_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "particle") & (pmb.df["particle_id"] == node_id[0])]["name"].values[0] + node_data = node_topology[node_id[0]] + node_name = node_data["particle_name"] + # Assert node's name and position are correctly set + np.testing.assert_equal(node_name_in_espresso, node_name) + np.testing.assert_allclose(node_pos, np.array(node_data["lattice_index"]) * 0.25 * diamond_lattice.BOXL, atol=1e-7) + + def test_chain_placement_and_connectivity(self): + for molecule_id, molecule_data in hydrogel_info["chains"].items(): + # Ensure that chain's node_start and node_end are correctly set + node_start = molecule_data["node_start"] + node_end = molecule_data["node_end"] + chain_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "molecule") & (pmb.df["molecule_id"] == molecule_id)]["name"].values[0] + # Assert chain's node_start and node_end + np.testing.assert_equal(chain_name_in_espresso, f"chain_{node_start}_{node_end}") + # Check if chain is connected in the espresso system (e.g., check bond or distance between node_start and node_end) + node_start_id = hydrogel_info["nodes"][node_start][0] + node_end_id = hydrogel_info["nodes"][node_end][0] + start_pos = espresso_system.part.by_id(int(node_start_id)).pos + end_pos = espresso_system.part.by_id(int(node_end_id)).pos + vec_between_nodes = end_pos - start_pos + # Ensure that the chain is connected (check distance, should be within acceptable bond length range) + vec_between_nodes = vec_between_nodes - diamond_lattice.BOXL * np.round(vec_between_nodes / diamond_lattice.BOXL) + distance_between_nodes = np.linalg.norm(vec_between_nodes) + np.testing.assert_allclose(distance_between_nodes, (diamond_lattice.MPC+1)*generic_bond_length.magnitude, atol=0.0000001) + + def test_all_residue_placement(self): + for chain_id, chain_data in hydrogel_info["chains"].items(): + node_start = chain_data["node_start"] + node_end = chain_data["node_end"] + + # Get the list of residues (keys in chain_data, excluding node_start/node_end) + residues_in_chain = {k: v for k, v in chain_data.items() if isinstance(k, int)} + + # Loop through each residue in the chain + for (res_id, res_data) in residues_in_chain.items(): + central_bead_id = res_data["central_bead_id"] + + # Get the position of the central bead from the espresso system + central_bead_pos = espresso_system.part.by_id(central_bead_id).pos + + # Calculate the expected position of the residue's central bead + residue_index = list(residues_in_chain.keys()).index(res_id) + + vec_between_nodes = (np.array([float(x) for x in node_end.strip('[]').split()]) - + np.array([float(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 = vec_between_nodes / (diamond_lattice.MPC + 1) + expected_position = (np.array([float(x) for x in node_start.strip('[]').split()]) * 0.25 * diamond_lattice.BOXL + + (residue_index + 1) * backbone_vector) + + # Validate that the central bead's position matches the expected position + np.testing.assert_allclose(central_bead_pos, expected_position, atol=0.05) + expected_res_name = chain_topology[chain_id]["residue_list"][residue_index] + expected_node_start = chain_topology[chain_id]["node_start"] + expected_node_end = chain_topology[chain_id]["node_end"] + residue_name = pmb.df[(pmb.df["pmb_type"]=="residue") & (pmb.df["residue_id"]==res_id)]["name"].values[0] + np.testing.assert_equal(node_start, expected_node_start) + np.testing.assert_equal(node_end, expected_node_end) + np.testing.assert_equal(residue_name, expected_res_name) + +if __name__ == "__main__": + ut.main() #####-- Invalid hydrogel name --##### # Test if create_hydrogel raises an exception when provided with invalid data print("*** Unit Test: Check invalid inputs for create_hydrogel ***") - -with ut.TestCase().assertRaises(ValueError): - pmb.create_hydrogel("invalid_hydrogel", espresso_system) # No such hydrogel in df - +np.testing.assert_raises(ValueError, pmb.create_hydrogel, "invalid_hydrogel", espresso_system) print("*** Invalid Input Test Passed ***") -####################################################################### # Check if the molecules (chains) are correctly stored in the hydrogel data for ((molecule_id, molecule_data),chain_id_in_input_dict) in zip(hydrogel_info["chains"].items(),chain_topology.items()): molecule_name_in_espresso = pmb.df[(pmb.df["pmb_type"] == "molecule") & (pmb.df["molecule_id"] == molecule_id)]["name"].values[0] - ut.TestCase().assertEqual(molecule_name_in_espresso, f"chain_{molecule_data['node_start']}_{molecule_data['node_end']}") + np.testing.assert_equal(molecule_name_in_espresso, f"chain_{molecule_data['node_start']}_{molecule_data['node_end']}") Residue_list_in_espresso = pmb.df[(pmb.df["pmb_type"]=="molecule") & (pmb.df["molecule_id"]==molecule_id)]["residue_list"].values[0] residue_list = chain_id_in_input_dict[1]["residue_list"] - ut.TestCase().assertEqual(Residue_list_in_espresso, residue_list) - -######################################################################## - - -#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 ***") + np.testing.assert_equal(Residue_list_in_espresso, residue_list) print("*** Checking if the ends of the randomly chosen chain is connected to node_start and node_end ***") random_chain_id = random.choice(list(chain_topology.keys())) # Choose a random chain of 0-15 diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index 81aa28c..ee24727 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -102,6 +102,8 @@ def test_lattice_setup(self): diamond = lib.lattice.DiamondLattice(MPC, BOND_LENGTH) espresso_system = espressomd.System(box_l = [diamond.BOXL]*3) pmb.add_bonds_to_espresso(espresso_system = espresso_system) + np.testing.assert_raises(ValueError, pmb.set_node, "[1 1 1]", NodeType1, espresso_system) + np.testing.assert_raises(ValueError, pmb.set_chain, "[0 0 0]", "[1 1 1]", {0:[0,0,0],1:diamond.BOXL/4.0*np.ones(3)},espresso_system) lattice = pmb.initialize_lattice_builder(diamond) sequence = [Res3, Res1, Res2, Res1] # build default structure @@ -122,6 +124,8 @@ def test_lattice_setup(self): np.testing.assert_equal(actual = lattice.get_node("[0 0 0]"), desired = NodeType2, verbose=True) pos_node3 = pmb.set_node("[2 2 0]", NodeType2, espresso_system=espresso_system) np.testing.assert_equal(actual = lattice.get_node("[2 2 0]"), desired = NodeType2, verbose=True) + pos_node4,_ = pmb.set_node("[3 1 3]", NodeType1, espresso_system=espresso_system) + np.testing.assert_equal(actual = lattice.get_node("[3 1 3]"), desired = NodeType1, verbose=True) node_positions={} node1_label = lattice.node_labels["[1 1 1]"] @@ -137,6 +141,10 @@ def test_lattice_setup(self): pmb.set_chain("[1 1 1]", "[0 0 0]", node_positions, espresso_system=espresso_system) np.testing.assert_equal(actual = lattice.get_chain("[1 1 1]", "[0 0 0]"), desired = sequence, verbose=True) np.testing.assert_equal(actual = lattice.get_chain("[0 0 0]", "[1 1 1]"), desired = sequence[::-1], verbose=True) + # set chain before set node + molecule_name = "chain_[3 1 3]_[0 0 0]" + pmb.define_molecule(name=molecule_name, residue_list=sequence) + np.testing.assert_raises(ValueError, pmb.set_chain, "[3 1 3]", "[0 0 0]", node_positions, espresso_system) # define custom chain in reverse direction molecule_name = "chain_[0 0 0]_[1 1 1]" From bf635a4c4442df398df535e51f99e5c0e1efbcd9 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Mon, 6 Jan 2025 20:07:21 +0100 Subject: [PATCH 12/12] refactoringhydrogel builder testsuite for code coverage --- testsuite/lattice_builder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index ee24727..7b9e5e3 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -124,7 +124,7 @@ def test_lattice_setup(self): np.testing.assert_equal(actual = lattice.get_node("[0 0 0]"), desired = NodeType2, verbose=True) pos_node3 = pmb.set_node("[2 2 0]", NodeType2, espresso_system=espresso_system) np.testing.assert_equal(actual = lattice.get_node("[2 2 0]"), desired = NodeType2, verbose=True) - pos_node4,_ = pmb.set_node("[3 1 3]", NodeType1, espresso_system=espresso_system) + _,_ = pmb.set_node("[3 1 3]", NodeType1, espresso_system=espresso_system) np.testing.assert_equal(actual = lattice.get_node("[3 1 3]"), desired = NodeType1, verbose=True) node_positions={}