From bb2d1d35fde949932ef4eb04fc3c24ffe97d9936 Mon Sep 17 00:00:00 2001 From: Marius Aarsten Date: Wed, 13 Mar 2024 14:17:47 +0100 Subject: [PATCH 01/10] implemented perp. function --- pyMBE.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index 8e17db2..b347ab0 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -1528,7 +1528,30 @@ def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_sample coord_list.append (coord) counter += 1 return coord_list + def generate_trial_perpendicular_vector(self,vector,center,radius): + """ + Generates an orthogonal vector with the same size as the input vector. + Args: + vector: Input vector for which an orthogonal vector needs to be generated. + + Returns: + perpendicular_vector: Orthogonal vector with the same magnitude as the input vector. + """ + np_vec = self.np.array(vector) + if self.np.all(np_vec == 0): + raise ValueError('Zero vector') + + # Generate a random vector with the same size as the input vector + random_vector = self.generate_trialvectors(center=center, radius=radius, 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) / self.np.dot(np_vec, 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) / self.np.linalg.norm(np_vec) + return perpendicular_vector + ''' def generate_trial_perpendicular_vector(self, vector, center, radius): """ Generates a random vector perpendicular to `vector`. @@ -1545,8 +1568,8 @@ def generate_trial_perpendicular_vector(self, vector, center, radius): 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 - + return center+norm_perp_vec*radius + ''' def generate_trialvectors(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 From 9876d54926491dc6807fc2c285e8fc96da5f4bcf Mon Sep 17 00:00:00 2001 From: blancoapa Date: Wed, 13 Mar 2024 15:55:36 +0100 Subject: [PATCH 02/10] fix bug in create_molecule(), fix bux in generate_trial_perpendicular(), improve docs --- pyMBE.py | 77 +++++++++++++++++++++++++++----------------------------- 1 file changed, 37 insertions(+), 40 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index b347ab0..f429b5d 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -678,11 +678,11 @@ def create_molecule(self, name, number_of_molecules, espresso_system, first_resi 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',''), @@ -706,10 +706,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','')): @@ -929,6 +930,7 @@ 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, radius=l0, @@ -938,7 +940,7 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead bead_position=self.generate_trial_perpendicular_vector(vector=backbone_vector, center=central_bead_position, radius=l0) - + side_bead_id = self.create_particle(name=side_chain_element, espresso_system=espresso_system, position=[bead_position], @@ -1528,48 +1530,43 @@ def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_sample coord_list.append (coord) counter += 1 return coord_list + def generate_trial_perpendicular_vector(self,vector,center,radius): """ - Generates an orthogonal vector with the same size as the input vector. + Generates an orthogonal vector to the input `vector`. + + Args: + vector(`lst`): arbitrary vector. + center(`lst`): origin of the orthogonal vector. + radius(`float`): magnitude of the orthogonal vector. + + Returns: + perpendicular_vector: Orthogonal vector with the same magnitude as the input vector. + """ + + # FIXME: + # We should consider renaming the input arguments for clarity + # center --> origin + # radius --> magnitude - Args: - vector: Input vector for which an orthogonal vector needs to be generated. - - Returns: - perpendicular_vector: Orthogonal vector with the same magnitude as the input vector. - """ np_vec = self.np.array(vector) + if self.np.all(np_vec == 0): raise ValueError('Zero vector') # Generate a random vector with the same size as the input vector - random_vector = self.generate_trialvectors(center=center, radius=radius, n_samples=1, on_surface=True)[0] + random_vector = self.generate_trialvectors(center=center, + radius=radius, + 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) / self.np.dot(np_vec, 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) / self.np.linalg.norm(np_vec) - return perpendicular_vector - ''' - 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 - ''' + perpendicular_vector /= self.np.linalg.norm(perpendicular_vector) + return center+perpendicular_vector*radius + def generate_trialvectors(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 @@ -1579,7 +1576,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: From 647772653cd8318b9cebe7d66e1643a419d6c80a Mon Sep 17 00:00:00 2001 From: Marius Aarsten Date: Thu, 14 Mar 2024 15:24:08 +0100 Subject: [PATCH 03/10] removed excessive normalization from projection in perp_function --- pyMBE.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index f429b5d..cc34314 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -940,7 +940,6 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead bead_position=self.generate_trial_perpendicular_vector(vector=backbone_vector, center=central_bead_position, radius=l0) - side_bead_id = self.create_particle(name=side_chain_element, espresso_system=espresso_system, position=[bead_position], @@ -1561,11 +1560,12 @@ def generate_trial_perpendicular_vector(self,vector,center,radius): on_surface=True)[0] # Project the random vector onto the input vector and subtract the projection - projection = self.np.dot(random_vector, np_vec) / self.np.dot(np_vec, np_vec) * np_vec + 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 center+perpendicular_vector*radius + def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface=False): """ From 382e81c872d0bf99bce5e8b5564d41ee473798d3 Mon Sep 17 00:00:00 2001 From: blancoapa Date: Thu, 14 Mar 2024 15:45:20 +0100 Subject: [PATCH 04/10] solve merge conflict --- pyMBE.py | 40 +++++++------------ samples/branched_polyampholyte.py | 4 +- samples/peptide.py | 14 +++---- samples/peptide_mixture_grxmc_ideal.py | 8 ++-- .../peptide_mixture_grxmc_unified_ideal.py | 8 ++-- 5 files changed, 32 insertions(+), 42 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index cc34314..da0ce19 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -1,4 +1,5 @@ class pymbe_library(): + """ The library for the Molecular Builder for ESPResSo (pyMBE) @@ -938,8 +939,9 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead on_surface=True)[0] else: bead_position=self.generate_trial_perpendicular_vector(vector=backbone_vector, - center=central_bead_position, - radius=l0) + origin=central_bead_position, + magnitude=l0) + side_bead_id = self.create_particle(name=side_chain_element, espresso_system=espresso_system, position=[bead_position], @@ -971,8 +973,8 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead on_surface=True)[0] else: residue_position=self.generate_trial_perpendicular_vector(vector=backbone_vector, - center=central_bead_position, - radius=l0) + origin=central_bead_position, + 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] @@ -1492,8 +1494,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)]] @@ -1530,42 +1531,32 @@ def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_sample counter += 1 return coord_list - def generate_trial_perpendicular_vector(self,vector,center,radius): + def generate_trial_perpendicular_vector(self,vector,origin,magnitude): """ Generates an orthogonal vector to the input `vector`. Args: vector(`lst`): arbitrary vector. - center(`lst`): origin of the orthogonal vector. - radius(`float`): magnitude of the orthogonal vector. + origin(`lst`): origin of the orthogonal vector. + magnitude(`float`): magnitude of the orthogonal vector. Returns: - perpendicular_vector: Orthogonal vector with the same magnitude as the input vector. + (`lst`): Orthogonal vector with the same magnitude as the input vector. """ - - # FIXME: - # We should consider renaming the input arguments for clarity - # center --> origin - # radius --> magnitude - - np_vec = self.np.array(vector) - + np_vec = self.np.array(vector) if self.np.all(np_vec == 0): raise ValueError('Zero vector') - # Generate a random vector with the same size as the input vector - random_vector = self.generate_trialvectors(center=center, - radius=radius, + random_vector = self.generate_trialvectors(center=[0,0,0], + radius=1, 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 center+perpendicular_vector*radius - + return origin+perpendicular_vector*magnitude def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface=False): """ @@ -1594,7 +1585,6 @@ 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] diff --git a/samples/branched_polyampholyte.py b/samples/branched_polyampholyte.py index f75d6a3..46dba8c 100644 --- a/samples/branched_polyampholyte.py +++ b/samples/branched_polyampholyte.py @@ -16,8 +16,8 @@ pmb = pyMBE.pymbe_library() # Load some functions from the handy_scripts library for convinience -from handy_scripts.handy_functions import setup_langevin_dynamics -from handy_scripts.handy_functions import block_analyze +from lib.handy_functions import setup_langevin_dynamics +from lib.analysis import block_analyze # The trajectories of the simulations will be stored using espresso built-up functions in separed files in the folder 'frames' if not os.path.exists('./frames'): diff --git a/samples/peptide.py b/samples/peptide.py index 8bf7434..815cc67 100644 --- a/samples/peptide.py +++ b/samples/peptide.py @@ -16,10 +16,10 @@ pmb = pyMBE.pymbe_library() # Load some functions from the handy_scripts library for convinience -from handy_scripts.handy_functions import setup_electrostatic_interactions -from handy_scripts.handy_functions import minimize_espresso_system_energy -from handy_scripts.handy_functions import setup_langevin_dynamics -from handy_scripts.handy_functions import block_analyze +from lib.handy_functions import setup_electrostatic_interactions +from lib.handy_functions import minimize_espresso_system_energy +from lib.handy_functions import setup_langevin_dynamics +from lib.analysis import block_analyze # The trajectories of the simulations will be stored using espresso built-up functions in separed files in the folder 'frames' if not os.path.exists('./frames'): @@ -39,15 +39,15 @@ # 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 # Load peptide parametrization from Lunkad, R. et al. Molecular Systems Design & Engineering (2021), 6(2), 122-131. -path_to_interactions=pmb.get_resource("reference_parameters/interaction_parameters/Lunkad2021.txt") -path_to_pka=pmb.get_resource("reference_parameters/pka_sets/Hass2015.txt") +path_to_interactions=pmb.get_resource("parameters/peptides/Lunkad2021.txt") +path_to_pka=pmb.get_resource("parameters/pka_sets/Hass2015.txt") pmb.load_interaction_parameters (filename=path_to_interactions) pmb.load_pka_set (path_to_pka) diff --git a/samples/peptide_mixture_grxmc_ideal.py b/samples/peptide_mixture_grxmc_ideal.py index 4a6ce4c..611a1d4 100644 --- a/samples/peptide_mixture_grxmc_ideal.py +++ b/samples/peptide_mixture_grxmc_ideal.py @@ -20,8 +20,8 @@ os.makedirs('./frames') #Import functions from handy_functions script -from handy_scripts.handy_functions import minimize_espresso_system_energy -from handy_scripts.handy_functions import block_analyze +from lib.handy_functions import minimize_espresso_system_energy +from lib.analysis import block_analyze # Simulation parameters pmb.set_reduced_units(unit_length=0.4*pmb.units.nm) @@ -46,8 +46,8 @@ N_peptide2_chains = 10 # Load peptide parametrization from Lunkad, R. et al. Molecular Systems Design & Engineering (2021), 6(2), 122-131. -path_to_interactions=pmb.get_resource("reference_parameters/interaction_parameters/Lunkad2021.txt") -path_to_pka=pmb.get_resource("reference_parameters/pka_sets/Hass2015.txt") +path_to_interactions=pmb.get_resource("parameters/peptides/Lunkad2021.txt") +path_to_pka=pmb.get_resource("parameters/pka_sets/Hass2015.txt") pmb.load_interaction_parameters (filename=path_to_interactions) pmb.load_pka_set (path_to_pka) diff --git a/samples/peptide_mixture_grxmc_unified_ideal.py b/samples/peptide_mixture_grxmc_unified_ideal.py index 570ce0c..9c6ac1d 100644 --- a/samples/peptide_mixture_grxmc_unified_ideal.py +++ b/samples/peptide_mixture_grxmc_unified_ideal.py @@ -20,8 +20,8 @@ os.makedirs('./frames') #Import functions from handy_functions script -from handy_scripts.handy_functions import minimize_espresso_system_energy -from handy_scripts.handy_functions import block_analyze +from lib.handy_functions import minimize_espresso_system_energy +from lib.analysis import block_analyze # Simulation parameters pmb.set_reduced_units(unit_length=0.4*pmb.units.nm) @@ -47,8 +47,8 @@ # Load peptide parametrization from Lunkad, R. et al. Molecular Systems Design & Engineering (2021), 6(2), 122-131. -path_to_interactions=pmb.get_resource("reference_parameters/interaction_parameters/Lunkad2021.txt") -path_to_pka=pmb.get_resource("reference_parameters/pka_sets/Hass2015.txt") +path_to_interactions=pmb.get_resource("parameters/peptides/Lunkad2021.txt") +path_to_pka=pmb.get_resource("parameters/pka_sets/Hass2015.txt") pmb.load_interaction_parameters (filename=path_to_interactions) pmb.load_pka_set (path_to_pka) From 383c01297ca16ac9c80bbe3b9a01de4b84149db4 Mon Sep 17 00:00:00 2001 From: blancoapa Date: Thu, 14 Mar 2024 15:48:45 +0100 Subject: [PATCH 05/10] update contributors list --- AUTHORS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/AUTHORS.md b/AUTHORS.md index 79d098f..b71cf3c 100644 --- a/AUTHORS.md +++ b/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) From 2727ec18506f1a09137d3f61fed24263bb475b0c Mon Sep 17 00:00:00 2001 From: Marius Aarsten Date: Tue, 2 Apr 2024 15:32:34 +0200 Subject: [PATCH 06/10] test script for gen_trial_perp --- testsuite/perpendicular_test.py | 50 +++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 testsuite/perpendicular_test.py diff --git a/testsuite/perpendicular_test.py b/testsuite/perpendicular_test.py new file mode 100644 index 0000000..93c0c84 --- /dev/null +++ b/testsuite/perpendicular_test.py @@ -0,0 +1,50 @@ +import sys +import os +import inspect +import numpy as np +import random +current_dir= os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) +path_end_index=current_dir.find("dendrimer_dna") + +code_path=current_dir[0:path_end_index]+"dendrimer_dna" +pyMBE_path = code_path + '/pyMBE' +sys.path.insert(0, pyMBE_path) +import pyMBE +pmb = pyMBE.pymbe_library() + + +#Defining variables for test +center = [0, 0, 0] +radius = 1 +n = 10 +tolerance = 1e-3 +original_vector = pmb.generate_trialvectors(center=center,radius=radius, n_samples=1, seed=None, on_surface=True)[0] +perpendicular_vectors = [] +for _ in range(n): + perpendicular_vector = pmb.generate_trial_perpendicular_vector(original_vector, center, radius) + perpendicular_vectors.append(perpendicular_vector) + +print(f"***Perpendicular vector unit test***") +print(f"***Unit test: Checks that the generated vector is perpendicular to the original vector {n} times.***") +for vector in perpendicular_vectors: + np.testing.assert_almost_equal(actual = np.dot(original_vector, vector), + desired = 0, + decimal = 5, + verbose = True) +print(f"***Unit test passed***") + +print(f"***Unit test: Checks that {n} generated perpendicular vectors are unique.***") +unique_set = set() +for vector in perpendicular_vectors: + vector_tuple = tuple(round(coord, 2) for coord in vector) + unique_set.add(vector_tuple) +np.testing.assert_equal(actual = len(unique_set), desired = n, verbose = True) +print(f"***Unit test passed***") + + + +print(f"***Unit test: Checks that {n} generated perpendicular vectors have the same magnitude as the input vector {original_vector} generated by pyMBE function.***") +for vector in perpendicular_vectors: + np.testing.assert_almost_equal(actual = np.linalg.norm(vector), desired = np.linalg.norm(original_vector), decimal = 5, verbose = True) +print(f"***Unit test passed") +print(f"*** All unit tests passed ***") \ No newline at end of file From fdde72dc713f5fee56209d12740399720f130321 Mon Sep 17 00:00:00 2001 From: "Pablo M. Blanco" Date: Tue, 2 Apr 2024 17:26:37 +0200 Subject: [PATCH 07/10] Add perpendicular test to the CI --- Makefile | 1 + ...pendicular_test.py => generate_perpendicular_vectors_test.py} | 0 2 files changed, 1 insertion(+) rename testsuite/{perpendicular_test.py => generate_perpendicular_vectors_test.py} (100%) diff --git a/Makefile b/Makefile index d1728df..1ef55a9 100644 --- a/Makefile +++ b/Makefile @@ -10,6 +10,7 @@ docs: tests: python3 testsuite/lj_tests.py python3 testsuite/peptide_tests.py + python3 testsuite/generate_perpendicular_vectors_test.py sample: python3 sample_scripts/peptide_simulation_example.py diff --git a/testsuite/perpendicular_test.py b/testsuite/generate_perpendicular_vectors_test.py similarity index 100% rename from testsuite/perpendicular_test.py rename to testsuite/generate_perpendicular_vectors_test.py From 8b7377c075e161e73024e4ef5d5ff462c26db89e Mon Sep 17 00:00:00 2001 From: "Pablo M. Blanco" Date: Tue, 2 Apr 2024 18:06:33 +0200 Subject: [PATCH 08/10] clean up test --- .../generate_perpendicular_vectors_test.py | 46 +++++++++---------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/testsuite/generate_perpendicular_vectors_test.py b/testsuite/generate_perpendicular_vectors_test.py index 93c0c84..3533fea 100644 --- a/testsuite/generate_perpendicular_vectors_test.py +++ b/testsuite/generate_perpendicular_vectors_test.py @@ -1,50 +1,48 @@ -import sys -import os -import inspect import numpy as np -import random -current_dir= os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) -path_end_index=current_dir.find("dendrimer_dna") - -code_path=current_dir[0:path_end_index]+"dendrimer_dna" -pyMBE_path = code_path + '/pyMBE' -sys.path.insert(0, pyMBE_path) import pyMBE pmb = pyMBE.pymbe_library() - #Defining variables for test -center = [0, 0, 0] -radius = 1 +origin = [1, 2, 3] +magnitude = 1 n = 10 tolerance = 1e-3 -original_vector = pmb.generate_trialvectors(center=center,radius=radius, n_samples=1, seed=None, on_surface=True)[0] +original_vector = pmb.generate_trialvectors(center=origin, + radius=magnitude, + n_samples=1, + seed=None, + on_surface=True)[0] perpendicular_vectors = [] for _ in range(n): - perpendicular_vector = pmb.generate_trial_perpendicular_vector(original_vector, center, radius) + perpendicular_vector = pmb.generate_trial_perpendicular_vector(vector=original_vector, + origin=origin, + magnitude=magnitude) perpendicular_vectors.append(perpendicular_vector) -print(f"***Perpendicular vector unit test***") -print(f"***Unit test: Checks that the generated vector is perpendicular to the original vector {n} times.***") +print(f"*** generate_trial_perpendicular_vector unit tests ***") +print(f"*** Unit test: Check that the function creates {n} perpendicular vectors to an input vector. ***") for vector in perpendicular_vectors: np.testing.assert_almost_equal(actual = np.dot(original_vector, vector), desired = 0, decimal = 5, verbose = True) -print(f"***Unit test passed***") -print(f"***Unit test: Checks that {n} generated perpendicular vectors are unique.***") +print(f"***Unit test passed***") +print(f"***Unit test: Check that the {n} perpendicular vectors are different.***") unique_set = set() for vector in perpendicular_vectors: vector_tuple = tuple(round(coord, 2) for coord in vector) unique_set.add(vector_tuple) -np.testing.assert_equal(actual = len(unique_set), desired = n, verbose = True) +np.testing.assert_equal(actual = len(unique_set), + desired = n, + verbose = True) print(f"***Unit test passed***") - - -print(f"***Unit test: Checks that {n} generated perpendicular vectors have the same magnitude as the input vector {original_vector} generated by pyMBE function.***") +print(f"***Unit test: Check that the perpendicular vectors have the same magnitude as the input magnitude.***") for vector in perpendicular_vectors: - np.testing.assert_almost_equal(actual = np.linalg.norm(vector), desired = np.linalg.norm(original_vector), decimal = 5, verbose = True) + np.testing.assert_almost_equal(actual = np.linalg.norm(vector), + desired = magnitude, + decimal = 5, + verbose = True) print(f"***Unit test passed") print(f"*** All unit tests passed ***") \ No newline at end of file From 367b226a393a863e98fabf965002134c0afd442e Mon Sep 17 00:00:00 2001 From: blancoapa Date: Tue, 2 Apr 2024 20:29:25 +0200 Subject: [PATCH 09/10] fix issue with origin, add additional tests with different magnitudes and origins --- Makefile | 2 +- pyMBE.py | 92 ++++++++-------- .../generate_perpendicular_vectors_test.py | 100 +++++++++++------- 3 files changed, 103 insertions(+), 91 deletions(-) diff --git a/Makefile b/Makefile index 1ef55a9..2972bd0 100644 --- a/Makefile +++ b/Makefile @@ -9,8 +9,8 @@ docs: tests: python3 testsuite/lj_tests.py - python3 testsuite/peptide_tests.py python3 testsuite/generate_perpendicular_vectors_test.py + python3 testsuite/peptide_tests.py sample: python3 sample_scripts/peptide_simulation_example.py diff --git a/pyMBE.py b/pyMBE.py index 5faaf71..d956745 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -695,7 +695,8 @@ 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] @@ -954,14 +955,13 @@ def create_residue(self, name, espresso_system, number_of_residues, central_bead 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, - origin=central_bead_position, - magnitude=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, @@ -988,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, - origin=central_bead_position, - magnitude=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] @@ -1226,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: @@ -1557,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`. @@ -1577,40 +1575,15 @@ def generate_coordinates_outside_sphere(self, center, radius, max_dist, n_sample coord_list = [] counter = 0 while counter=radius: coord_list.append (coord) counter += 1 return coord_list - def generate_trial_perpendicular_vector(self,vector,origin,magnitude): - """ - Generates an orthogonal vector to the input `vector`. - - Args: - vector(`lst`): arbitrary vector. - origin(`lst`): origin of the orthogonal 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) - if self.np.all(np_vec == 0): - raise ValueError('Zero vector') - # Generate a random vector with the same size as the input vector - random_vector = self.generate_trialvectors(center=[0,0,0], - radius=1, - 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 origin+perpendicular_vector*magnitude - - 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. @@ -1642,11 +1615,37 @@ def generate_trialvectors(self, center, radius, n_samples, seed=None, on_surface 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. @@ -2679,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',''), @@ -2707,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): @@ -2802,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): diff --git a/testsuite/generate_perpendicular_vectors_test.py b/testsuite/generate_perpendicular_vectors_test.py index 3533fea..f8e02b7 100644 --- a/testsuite/generate_perpendicular_vectors_test.py +++ b/testsuite/generate_perpendicular_vectors_test.py @@ -1,48 +1,66 @@ import numpy as np import pyMBE +from itertools import combinations pmb = pyMBE.pymbe_library() -#Defining variables for test -origin = [1, 2, 3] -magnitude = 1 -n = 10 -tolerance = 1e-3 -original_vector = pmb.generate_trialvectors(center=origin, - radius=magnitude, - n_samples=1, - seed=None, - on_surface=True)[0] -perpendicular_vectors = [] -for _ in range(n): - perpendicular_vector = pmb.generate_trial_perpendicular_vector(vector=original_vector, - origin=origin, - magnitude=magnitude) - perpendicular_vectors.append(perpendicular_vector) +def check_if_different_perpendicular_vectors_are_generated(vector,magnitude,n=50): + """ + Checks if pmb.generate_trial_perpendicular_vector generates perpendicular vectors to `vector`. -print(f"*** generate_trial_perpendicular_vector unit tests ***") -print(f"*** Unit test: Check that the function creates {n} perpendicular vectors to an input vector. ***") -for vector in perpendicular_vectors: - np.testing.assert_almost_equal(actual = np.dot(original_vector, vector), - desired = 0, - decimal = 5, - verbose = True) + Args: + vector(`array`,`float`): vector to which new perpendicular vectors will be generated. + magnitude(`float`): magnitude of the perpendicular vectors. + n(`int`): number of perpendicular vectors to be generated. + """ + perpendicular_vectors = [] + for _ in range(n): + perpendicular_vector = pmb.generate_trial_perpendicular_vector(vector=vector, + magnitude=magnitude) + perpendicular_vectors.append(perpendicular_vector) + # Check that the function creates {n} perpendicular vectors to an input vector + for pvector in perpendicular_vectors: + np.testing.assert_almost_equal(actual = np.dot(pvector, vector), + desired = 0, + decimal = 5, + verbose = True) + # Check that the {n} perpendicular vectors are different + for vector_pair in combinations(perpendicular_vectors, 2): + if np.array_equal(vector_pair[0],vector_pair[1]): + raise Exception(f"Error: pmb.generate_trial_perpendicular_vector two equal perpendicular vectors v1 = {vector_pair[0]} v2 = {vector_pair[1]}") + # Check that the perpendicular vectors have the same magnitude as the input magnitude + for pvector in perpendicular_vectors: + np.testing.assert_almost_equal(actual = np.linalg.norm(pvector), + desired = magnitude, + decimal = 5, + verbose = True) + return -print(f"***Unit test passed***") -print(f"***Unit test: Check that the {n} perpendicular vectors are different.***") -unique_set = set() -for vector in perpendicular_vectors: - vector_tuple = tuple(round(coord, 2) for coord in vector) - unique_set.add(vector_tuple) -np.testing.assert_equal(actual = len(unique_set), - desired = n, - verbose = True) -print(f"***Unit test passed***") - -print(f"***Unit test: Check that the perpendicular vectors have the same magnitude as the input magnitude.***") -for vector in perpendicular_vectors: - np.testing.assert_almost_equal(actual = np.linalg.norm(vector), - desired = magnitude, - decimal = 5, - verbose = True) -print(f"***Unit test passed") +print(f"*** generate_trial_perpendicular_vector unit tests ***") +print(f"*** Unit test: Check that the function creates perpendicular vectors to an arbitrary vector of the same magnitude ***") +vector = pmb.generate_random_points_in_a_sphere(center=[0,0,0], + radius=1, + n_samples=1, + seed=None, + on_surface=True)[0] +check_if_different_perpendicular_vectors_are_generated(vector=vector, + magnitude=1) +print(f"*** Unit test passed ***") +print(f"*** Unit test: Check that the function creates perpendicular vectors to a vector with an arbitrary origin ***") +vector = pmb.generate_random_points_in_a_sphere(center=[1,2,3], + radius=1, + n_samples=1, + seed=None, + on_surface=True)[0] +check_if_different_perpendicular_vectors_are_generated(vector=vector, + magnitude=1) +print(f"*** Unit test passed ***") +print(f"*** Unit test: Check that the function creates perpendicular vectors with a different magnitude ***") +vector = pmb.generate_random_points_in_a_sphere(center=[1,2,3], + radius=2, + n_samples=1, + seed=None, + on_surface=True)[0] +check_if_different_perpendicular_vectors_are_generated(vector=vector, + magnitude=2) +print(f"*** Unit test passed ***") print(f"*** All unit tests passed ***") \ No newline at end of file From cc9dbcfbef730723c87355295229cdfd477f822b Mon Sep 17 00:00:00 2001 From: blancoapa Date: Tue, 2 Apr 2024 20:40:54 +0200 Subject: [PATCH 10/10] fix test to have different magnitude for the input and perpendicular vectors --- testsuite/generate_perpendicular_vectors_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/generate_perpendicular_vectors_test.py b/testsuite/generate_perpendicular_vectors_test.py index f8e02b7..1ab0a0e 100644 --- a/testsuite/generate_perpendicular_vectors_test.py +++ b/testsuite/generate_perpendicular_vectors_test.py @@ -61,6 +61,6 @@ def check_if_different_perpendicular_vectors_are_generated(vector,magnitude,n=50 seed=None, on_surface=True)[0] check_if_different_perpendicular_vectors_are_generated(vector=vector, - magnitude=2) + magnitude=3) print(f"*** Unit test passed ***") print(f"*** All unit tests passed ***") \ No newline at end of file