include backbone vector in create_molecule (#99)

* include backbone vector in create_molecule

* include backbone vector in create molecule

* Test for backbone vector in create_molecule()

* backbone_vector added in the docstring and node_map/chain_map removed from setup_df

* Add somesh to the list of authors, update changelog, include backbone_vector to pmb.create_pmb_objetct

---------

Co-authored-by: pm-blanco <pablb@ntnu.no>

Author: 1234somesh

AUTHORS.md

@@ -12,6 +12,7 @@ The following people have contributed to pyMBE (Github username, Name, current a
 - pinedaps, Sebastian P. Pineda (Charles University, Czech Republic)
 
 ## Tier-1 Contributors
+- 1234somesh, Somesh Kurahatti (University of Stuttgart, Germany)
 - amartise, Albert Martinez-Serra (Royal College of Surgeons, Ireland) 
 - mariusaarsten, Marius Aarsten (Norwegian University of Science and Tecnology, Norway)
 - Zitzeronion, Stefan Zitz (Research Center Pharmaceutical Engineering, Austria)

CHANGELOG.md

@@ -17,6 +17,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Switched to Ctest for testing, allowing to run the tests on paralel (#87)
 
 ### Added
+- New optional argument `backbone_vector` enabling to build molecules along an input vector using `pmb.create_molecule` and `pmb.create_pmb_object` (#99)
+- Unit testing for reaction methods (#86)
 - New boolean flag `--ideal` as argparse argument of `samples/globular_protein.py` enabling to run the script without setting up interactions.
 - Unit tests for `pmb.create_protein`, `pmb.enable_motion_of_rigid_object`, `pmb.protein_sequence_parser`, `pmb.define_protein`, `pmb.read_protein_vtf_in_df` (#101)
 - Library `lattice.py`, a general builder for crystaline lattices. This library is part of on-going project to support hydrogels in pyMBE. (#93)

pyMBE.py

@@ -873,15 +873,16 @@ def create_counterions(self, object_name, cation_name, anion_name, espresso_syst
                 print(f'Ion type: {name} created number: {counterion_number[name]}')
         return counterion_number
 
-    def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, use_default_bond=False):
+    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`.
 
         Args:
             name(`str`): Label of the molecule type to be created. `name` must be defined in `pmb.df`
             espresso_system(`espressomd.system.System`): Instance of a system object from espressomd library.
             number_of_molecules(`int`): Number of molecules of type `name` to be created.
-            list_of_first_residue_positions(`list`, optional): List of coordinates where the central bead of the first_residue_position will be created, random by default
+            list_of_first_residue_positions(`list`, optional): List of coordinates where the central bead of the first_residue_position will be created, random by default.
+            backbone_vector(`list` of `float`): Backbone vector of the molecule, random by default. Central beads of the residues in the `residue_list` are placed along this vector. 
             use_default_bond(`bool`, optional): Controls if a bond of type `default` is used to bond particle with undefined bonds in `pymbe.df`
 
         Returns:
@@ -926,7 +927,8 @@ 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], 
+                    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]
@@ -1039,7 +1041,7 @@ def create_particle(self, name, espresso_system, number_of_particles, position=N
             self.add_value_to_df(key=('particle_id',''),index=df_index,new_value=bead_id, verbose=False)                  
         return created_pid_list
 
-    def create_pmb_object(self, name, number_of_objects, espresso_system, position=None, use_default_bond=False):
+    def create_pmb_object(self, name, number_of_objects, espresso_system, position=None, use_default_bond=False, backbone_vector=None):
         """
         Creates all `particle`s associated to `pmb object` into  `espresso` a number of times equal to `number_of_objects`.
         
@@ -1049,6 +1051,7 @@ def create_pmb_object(self, name, number_of_objects, espresso_system, position=N
             espresso_system(`espressomd.system.System`): Instance of an espresso system object from espressomd library.
             position(`list`): Coordinates where the particles should be created.
             use_default_bond(`bool`,optional): Controls if a `default` bond is used to bond particles with undefined bonds in `pmb.df`. Defaults to `False`.
+            backbone_vector(`list` of `float`): Backbone vector of the molecule, random by default. Central beads of the residues in the `residue_list` are placed along this vector. 
 
         Note:
             - If no `position` is given, particles will be created in random positions. For bonded particles, they will be created at a distance equal to the bond length. 
@@ -1066,13 +1069,15 @@ def create_pmb_object(self, name, number_of_objects, espresso_system, position=N
             self.create_residue(name=name,  
                                 espresso_system=espresso_system, 
                                 central_bead_position=position,
-                                use_default_bond=use_default_bond)
+                                use_default_bond=use_default_bond,
+                                backbone_vector=backbone_vector)
         elif pmb_type == 'molecule':
             self.create_molecule(name=name, 
                                 number_of_molecules=number_of_objects, 
                                 espresso_system=espresso_system, 
                                 use_default_bond=use_default_bond, 
-                                list_of_first_residue_positions=position)
+                                list_of_first_residue_positions=position,
+                                backbone_vector=backbone_vector)
         return
 
     def create_protein(self, name, number_of_proteins, espresso_system, topology_dict):
@@ -1855,7 +1860,7 @@ def find_bond_key(self, particle_name1, particle_name2, use_default_bond=False):
         bond_keys = [particle_name1 +'-'+ particle_name2, particle_name2 +'-'+ particle_name1 ]
         bond_defined=False
         for bond_key in bond_keys:
-            if bond_key in self.df.values:
+            if bond_key in self.df["name"].values:
                 bond_defined=True
                 correct_key=bond_key
                 break
@@ -3221,7 +3226,7 @@ def setup_df (self):
                 '': object},
             'l0': {
                 '': float},
-        }
+            }
 
         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()]))
 

testsuite/define_and_create_molecules_unit_tests.py

@@ -325,13 +325,18 @@
 # Additional unit tests for define_molecule are in create_molecule_position_test
 print("*** Unit test: check that create_molecule() creates a simple molecule into the espresso_system with the properties defined in pmb.df  ***")
 
-pmb.create_molecule(name="M2",
+backbone_vector = np.array([1,3,-4])
+magnitude = np.linalg.norm(backbone_vector)
+backbone_vector = backbone_vector/magnitude
+molecule_info_M2 = pmb.create_molecule(name="M2",
                     number_of_molecules=2,
                     espresso_system=espresso_system,
+                    backbone_vector = backbone_vector,
                     use_default_bond=True)
 
 particle_ids=pmb.get_particle_id_map(object_name="M2")["all"]
 
+
 residue_ids={9: 2, 10: 3, 11: 3, 12: 3, 13: 4, 14: 4, 15: 4, 16: 4, # First molecule
             17: 5, 18: 6, 19: 6, 20: 6, 21: 7, 22: 7, 23: 7, 24: 7} # Second molecule
 
@@ -412,6 +417,41 @@
                             verbose=True)
 
 print("*** Unit test passed ***")
+
+print("*** Unit test: check the backbone vector of the molecule in espresso and the given input backbone vector are same ***")
+
+
+central_bead_positions = []
+
+for residue_name in molecule_parameters["M2"]["residue_list"]:
+    
+    mol_id = pmb.df[pmb.df["name"]=="M2"]["molecule_id"].values[0]
+    res_id = pmb.df[(pmb.df["molecule_id"]==mol_id) & (pmb.df['name']==residue_name)]["residue_id"].values[0]
+    central_bead_id = molecule_info_M2[mol_id][res_id]['central_bead_id']
+    central_bead_pos = espresso_system.part.by_id(central_bead_id).pos
+    central_bead_positions.append(central_bead_pos)
+
+if len(central_bead_positions) == len(molecule_parameters["M2"]["residue_list"]):
+    
+    backbone_direction_1 = central_bead_positions[1] - central_bead_positions[0]
+    backbone_direction_2 = central_bead_positions[2] - central_bead_positions[1]
+    backbone_direction_1 /= np.linalg.norm(backbone_direction_1)
+    backbone_direction_2 /= np.linalg.norm(backbone_direction_2)
+    np.testing.assert_almost_equal(
+        actual = backbone_direction_1,
+        desired = backbone_vector,
+        verbose = True)
+    np.testing.assert_almost_equal(
+        actual = backbone_direction_2,
+        desired = backbone_vector,
+        verbose = True)
+
+else:
+
+    raise ValueError("Expected 3 central bead positions for residues R1, R2, and R3")        
+
+print("*** Unit test passed ***")
+
 print("*** Unit test: check that create_molecule() does not create any molcule for number_of_molecules <= 0  ***")
 
 starting_number_of_particles=len(espresso_system.part.all())