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: diff --git a/lib/lattice.py b/lib/lattice.py index 8b1e942..7e97690 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. @@ -20,17 +19,16 @@ import itertools import numpy as np +#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,6 +44,8 @@ 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 @@ -58,6 +58,9 @@ def __init__(self, lattice, strict=True): self.kwargs_monomers = {} self.kwargs_bonds = {} self.kwargs_box = {} + self.MPC = lattice.MPC + self.BOXL = lattice.BOXL + 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 +96,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 +106,6 @@ def add_default_chains(self, mpc): def draw_lattice(self, ax): """ Draw the lattice in an `Axes3D `_ canvas. - Args: ax: Axes. """ @@ -150,8 +151,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 +166,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 +184,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 +202,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 +214,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 +226,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 +244,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 +258,11 @@ class DiamondLattice: (2, 5), (3, 6), (4, 7), (5, 0), (5, 3), (5, 4), (6, 0), (6, 2), (6, 4), (7, 0), (7, 2), (7, 3)} + + def __init__(self,MPC,BOND_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 b1fcdcb..8971fc3 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, DiamondLattice class pymbe_library(): @@ -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,7 +874,46 @@ 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 + 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: + + """ + 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 = {} + 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, 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 + # 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"] + 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): """ Creates `number_of_molecules` molecule of type `name` into `espresso_system` and bookkeeps them into `pmb.df`. @@ -929,9 +970,10 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi # Generate an arbitrary random unit vector if backbone_vector is None: backbone_vector = self.generate_random_points_in_a_sphere(center=[0,0,0], - radius=1, + 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, @@ -1463,6 +1505,55 @@ 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": , ... } + + """ + 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: + raise ValueError( + 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: + raise ValueError( + 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 + + 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.delete_molecule_entry(molecule_name) + self.define_molecule(name=molecule_name, residue_list=residue_list) + return; def define_molecule(self, name, residue_list): """ @@ -1673,7 +1764,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): """ @@ -1894,6 +1995,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`. @@ -2215,6 +2320,16 @@ 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. + """ + 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 + def load_interaction_parameters(self, filename, verbose=False, overwrite=False): """ Loads the interaction parameters stored in `filename` into `pmb.df` @@ -2605,6 +2720,126 @@ 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] + + 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 + 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) + return chain_molecule_info + + 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 + 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(), p_id + 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. @@ -3226,6 +3461,10 @@ 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..32fe669 --- /dev/null +++ b/samples/hydrogel.py @@ -0,0 +1,86 @@ +import pyMBE +import espressomd +import matplotlib.pyplot as plt +from lib.lattice import DiamondLattice + +pmb = pyMBE.pymbe_library(seed=42) +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) +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) # Dont need pmb when integrated to pyMBE.py + +# 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 = {} + +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 + +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) + +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 435309b..78745c1 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/hydrogel_builder.py b/testsuite/hydrogel_builder.py new file mode 100644 index 0000000..e244760 --- /dev/null +++ b/testsuite/hydrogel_builder.py @@ -0,0 +1,331 @@ +import numpy as np +import random +import unittest as ut +import pyMBE +from lib.lattice import DiamondLattice +import espressomd + +pmb = pyMBE.pymbe_library(seed=42) + +# 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 +) +# 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, + '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 ***") +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 ***") + +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 ***") + + +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))] +################################################ +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}} + +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()} + 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) +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." +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 + +# Creating hydrogel +hydrogel_info = pmb.create_hydrogel("my_hydrogel", espresso_system) +print("*** Hydrogel created: Unit test to verify their name and positions ***") + +############################ + +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 ***") +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] + 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"] + 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 +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("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 1481d70..7b9e5e3 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -16,105 +16,181 @@ # 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) + 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 - 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) + _,_ = 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]"] + node_positions[node1_label]=pos_node1[0] + node2_label = lattice.node_labels["[0 0 0]"] + node_positions[node2_label]=pos_node2[0] + node3_label = lattice.node_labels["[2 2 0]"] + node_positions[node3_label]=pos_node3[0] + + # 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) + # 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 - 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 + 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(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) + 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 +199,6 @@ def test_plot(self): ax.legend() plt.close(fig) - if __name__ == "__main__": ut.main() + 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)