From 367b226a393a863e98fabf965002134c0afd442e Mon Sep 17 00:00:00 2001 From: blancoapa Date: Tue, 2 Apr 2024 20:29:25 +0200 Subject: [PATCH] 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