Merge pull request #26 from mariusaarsten/fixing_perpendicular_vector

Fixing perpendicular vector

AUTHORS.md

@@ -10,12 +10,13 @@ The following people have contributed to pyMBE:
 - Paola B. Torres (Universidad Tecnologica Nacional, Argentina)
 
 ## Contributors
+- Marius Aarsten (Norwegian University of Science and Tecnology, Norway)
 - Claudio F. Narambuena (Universidad Tecnologica Nacional, Argentina)
 
 ## Early testers of the library
 - Corinna Dannert (Norwegian University of Science and Technology, Norway)
 - Rita S. Dias (Norwegian University of Science and Technology, Norway)
--  Sergio Madurga (University of Barcelona, Spain)
+- Sergio Madurga (University of Barcelona, Spain)
 - Albert Martinez-Serra (Royal College of Surgeons, Ireland)
 - Francesc Mas (University of Barcelona, Spain)
 - Magdaléna Nejedlá (Charles University, Czech Republic)

Makefile

@@ -9,6 +9,7 @@ docs:
 
 tests:
 	python3 testsuite/lj_tests.py
+	python3 testsuite/generate_perpendicular_vectors_test.py
 	python3 testsuite/peptide_tests.py
 
 sample:

pyMBE.py

@@ -1,4 +1,5 @@
 class pymbe_library():
+
     """
     The library for the Molecular Builder for ESPResSo (pyMBE)
 
@@ -694,16 +695,17 @@ def create_molecule(self, name, number_of_molecules, espresso_system, first_resi
             for residue in residue_list:
                 if first_residue:
                     residue_position = first_residue_position
-                    backbone_vector = self.generate_trialvectors(center=[0,0,0], 
+                    # Generate an arbitrary random unit vector
+                    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,
-                                                                    number_of_residues=1, 
-                                                                    espresso_system=espresso_system, 
-                                                                    central_bead_position=first_residue_position,  
-                                                                    use_default_bond= use_default_bond, 
-                                                                    backbone_vector=backbone_vector)
+                                                        number_of_residues=1, 
+                                                        espresso_system=espresso_system, 
+                                                        central_bead_position=first_residue_position,  
+                                                        use_default_bond= use_default_bond, 
+                                                        backbone_vector=backbone_vector)
                     residue_id = next(iter(residues_info))
                     for index in self.df[self.df['residue_id']==residue_id].index:
                         self.add_value_to_df(key=('molecule_id',''),
@@ -727,10 +729,11 @@ def create_molecule(self, name, number_of_molecules, espresso_system, first_resi
                                               use_default_bond=use_default_bond)                
                     residue_position = residue_position+backbone_vector*l0
                     residues_info = self.create_residue(name=residue, 
-                                                                    number_of_residues=1, 
-                                                                    espresso_system=espresso_system, 
-                                                                    central_bead_position=[residue_position],
-                                                                    use_default_bond= use_default_bond)
+                                                        number_of_residues=1, 
+                                                        espresso_system=espresso_system, 
+                                                        central_bead_position=[residue_position],
+                                                        use_default_bond= use_default_bond, 
+                                                        backbone_vector=backbone_vector)
                     residue_id = next(iter(residues_info))      
                     for index in self.df[self.df['residue_id']==residue_id].index:
                         if not self.check_if_df_cell_has_a_value(index=index,key=('molecule_id','')):
@@ -950,16 +953,16 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead
                                               particle_name2=side_chain_element, 
                                               hard_check=True, 
                                               use_default_bond=use_default_bond)
+
                     if backbone_vector is None:
-                        bead_position=self.generate_trialvectors(center=central_bead_position, 
+                        bead_position=self.generate_random_points_in_a_sphere(center=central_bead_position, 
                                                                  radius=l0, 
                                                                  n_samples=1,
                                                                  on_surface=True)[0]
                     else:
-                        bead_position=self.generate_trial_perpendicular_vector(vector=backbone_vector,
-                                                                            center=central_bead_position, 
-                                                                            radius=l0)
-
+                        bead_position=central_bead_position+self.generate_trial_perpendicular_vector(vector=backbone_vector,
+                                                                                                    magnitude=l0)
+
                     side_bead_id = self.create_particle(name=side_chain_element, 
                                                                     espresso_system=espresso_system,
                                                                     position=[bead_position], 
@@ -985,14 +988,13 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead
                                               hard_check=True, 
                                               use_default_bond=use_default_bond)
                     if backbone_vector is None:
-                        residue_position=self.generate_trialvectors(center=central_bead_position, 
+                        residue_position=self.generate_random_points_in_a_sphere(center=central_bead_position, 
                                                                     radius=l0, 
                                                                     n_samples=1,
                                                                     on_surface=True)[0]
                     else:
-                        residue_position=self.generate_trial_perpendicular_vector(vector=backbone_vector,
-                                                                                center=central_bead_position, 
-                                                                                radius=l0)
+                        residue_position=central_bead_position+self.generate_trial_perpendicular_vector(vector=backbone_vector,
+                                                                                                        magnitude=l0)
                     lateral_residue_info = self.create_residue(name=side_chain_element, espresso_system=espresso_system,
                         number_of_residues=1, central_bead_position=[residue_position],use_default_bond=use_default_bond)
                     lateral_residue_dict=list(lateral_residue_info.values())[0]
@@ -1223,9 +1225,9 @@ def define_particle(self, name, q=0, acidity='inert', pka=None, sigma=None, epsi
             q (`int`, optional): Permanent charge of this particle type. Defaults to 0.
             acidity (`str`, optional): Identifies whether if the particle is `acidic` or `basic`, used to setup constant pH simulations. Defaults to `inert`.
             pka (`float`, optional): If `particle` is an acid or a base, it defines its  pka-value. Defaults to None.
-            sigma (`pint.Quantity`, optional): Sigma parameters used to set up Lennard-Jones interactions for this particle type. Defaults to None.
-            cutoff (`pint.Quantity`, optional): Sigma parameters used to set up Lennard-Jones interactions for this particle type. Defaults to None.
-            offset (`pint.Quantity`, optional): Sigma parameters used to set up Lennard-Jones interactions for this particle type. Defaults to None.
+            sigma (`pint.Quantity`, optional): Sigma parameter used to set up Lennard-Jones interactions for this particle type. Defaults to None.
+            cutoff (`pint.Quantity`, optional): cutoff parameter used to set up Lennard-Jones interactions for this particle type. Defaults to None.
+            offset (`pint.Quantity`, optional): offset parameter used to set up Lennard-Jones interactions for this particle type. Defaults to None.
             epsilon (`pint.Quantity`, optional): Epsilon parameter used to setup Lennard-Jones interactions for this particle tipe. Defaults to None.
 
         Note:
@@ -1543,8 +1545,7 @@ def find_value_from_es_type(self, es_type, column_name):
         """
         idx = self.pd.IndexSlice
         for state in ['state_one', 'state_two']:
-            index = self.np.where(self.df[(state, 'es_type')] == es_type)[0]
-
+            index = self.np.where(self.df[(state, 'es_type')] == es_type)[0]      
             if len(index) > 0:
                 if column_name == 'label':
                     label = self.df.loc[idx[index[0]], idx[(state,column_name)]]
@@ -1555,7 +1556,6 @@ def find_value_from_es_type(self, es_type, column_name):
         return None
 
     def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_samples):
-
         """
         Generates coordinates outside a sphere centered at `center`.
 
@@ -1575,31 +1575,15 @@ def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_sample
         coord_list = []
         counter = 0
         while counter<n_samples:
-            coord = self.generate_trialvectors(center=center, radius=max_dist,n_samples=1)[0]
+            coord = self.generate_random_points_in_a_sphere(center=center, 
+                                            radius=max_dist,
+                                            n_samples=1)[0]
             if self.np.linalg.norm(coord-self.np.asarray(center))>=radius:
                 coord_list.append (coord)
                 counter += 1
         return coord_list
-
-    def generate_trial_perpendicular_vector(self, vector, center, radius):
-        """
-        Generates a random vector perpendicular to `vector`.
-        
-        Args:
-            vector(`lst`): Coordinates of the vector.
-            magnitude(`float`): magnitude of the vector perpendicular to `vector`.
-
-        Returns:
-            perpendicular_vector(`np.array`): random vector perpendicular to `vector`.
-        """
-        np_vec=self.np.array(vector)
-        if np_vec[1] == 0 and np_vec[2] == 0 and np_vec[0] == 1:
-            raise ValueError('zero vector')
-        perp_vec = self.np.cross(np_vec, self.generate_trialvectors(center=center, radius=radius, n_samples=1, on_surface=True)[0])
-        norm_perp_vec = perp_vec/self.np.linalg.norm(perp_vec)
-        return center+norm_perp_vec*radius  
 
-    def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface=False):
+    def generate_random_points_in_a_sphere(self, center, radius, n_samples, seed=None, on_surface=False):
         """
         Uniformly samples points from a hypersphere. If on_surface is set to True, the points are
         uniformly sampled from the surface of the hypersphere.
@@ -1608,7 +1592,7 @@ def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface
             center(`lst`): Array with the coordinates of the center of the spheres.
             radius(`float`): Radius of the sphere.
             n_samples(`int`): Number of sample points to generate inside the sphere.
-            seed (`int`, optional): Seed for the random number generator
+            seed (`int`, optional): Seed for the random number generator.
             on_surface (`bool`, optional): If set to True, points will be uniformly sampled on the surface of the hypersphere.
 
         Returns:
@@ -1626,17 +1610,42 @@ def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface
         # make the samples lie on the surface of the unit hypersphere
         normalize_radii = self.np.linalg.norm(samples, axis=1)[:, self.np.newaxis]
         samples /= normalize_radii
-
         if not on_surface:
             # make the samples lie inside the hypersphere with the correct density
             uniform_points = rng.uniform(size=n_samples)[:, self.np.newaxis]
             new_radii = self.np.power(uniform_points, 1/d)
             samples *= new_radii
-
         # scale the points to have the correct radius and center
         samples = samples * radius + center
         return samples 
 
+    def generate_trial_perpendicular_vector(self,vector,magnitude):
+        """
+        Generates an orthogonal vector to the input `vector`.
+
+        Args:
+            vector(`lst`): arbitrary vector.
+            magnitude(`float`): magnitude of the orthogonal vector.
+            
+        Returns:
+            (`lst`): Orthogonal vector with the same magnitude as the input vector.
+        """ 
+        np_vec = self.np.array(vector) 
+        np_vec /= self.np.linalg.norm(np_vec) 
+        if self.np.all(np_vec == 0):
+            raise ValueError('Zero vector')
+        # Generate a random vector 
+        random_vector = self.generate_random_points_in_a_sphere(radius=1, 
+                                                                center=[0,0,0],
+                                                                n_samples=1, 
+                                                                on_surface=True)[0]
+        # Project the random vector onto the input vector and subtract the projection
+        projection = self.np.dot(random_vector, np_vec) * np_vec
+        perpendicular_vector = random_vector - projection
+        # Normalize the perpendicular vector to have the same magnitude as the input vector
+        perpendicular_vector /= self.np.linalg.norm(perpendicular_vector) 
+        return perpendicular_vector*magnitude
+
     def get_bond_length(self, particle_name1, particle_name2, hard_check=False, use_default_bond=False) :
         """
         Searches for bonds between the particle types given by `particle_name1` and `particle_name2` in `pymbe.df` and returns the initial bond length.
@@ -2669,7 +2678,6 @@ def setup_df (self):
         Returns:
             columns_names(`obj`): pandas multiindex object with the column names of the pyMBE's dataframe
         """
-
         columns_names = self.pd.MultiIndex.from_tuples ([
                         ('name',''),
                         ('pmb_type',''),
@@ -2697,9 +2705,7 @@ def setup_df (self):
                         ('bond_object',''),
                         ('parameters_of_the_potential','')
                         ])
-
         self.df = self.pd.DataFrame (columns = columns_names)
-
         return columns_names
 
     def setup_lj_interactions(self, espresso_system, shift_potential=True, combining_rule='Lorentz-Berthelot', warnings=True):
@@ -2792,11 +2798,9 @@ def setup_particle_sigma(self, topology_dict):
             residue_sigma  = topology_dict[residue]['sigma']
             sigma = residue_sigma*self.units.nm
             index = self.df[(self.df['residue_id']==residue_number) & (self.df['name']==residue_name) ].index.values[0]
-
             self.add_value_to_df(key= ('sigma',''),
                         index=int (index),
-                        new_value=sigma)
-
+                        new_value=sigma)           
         return 
 
     def write_output_vtf_file(self, espresso_system, filename):

samples/peptide.py

@@ -39,7 +39,7 @@
 
 # Peptide parameters
 
-sequence = 'nGEGGHc'
+sequence = 'nEEEEEc'
 model = '2beadAA'  # Model with 2 beads per each aminoacid
 pep_concentration = 5.56e-4 *pmb.units.mol/pmb.units.L
 N_peptide_chains = 4