Skip to content

Commit

Permalink
fix issue with origin, add additional tests with different magnitudes…
Browse files Browse the repository at this point in the history
… and origins
  • Loading branch information
pm-blanco committed Apr 2, 2024
1 parent 8b7377c commit 367b226
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 91 deletions.
2 changes: 1 addition & 1 deletion Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
92 changes: 43 additions & 49 deletions pyMBE.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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,
Expand All @@ -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]
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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`.
Expand All @@ -1577,40 +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,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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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',''),
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down
100 changes: 59 additions & 41 deletions testsuite/generate_perpendicular_vectors_test.py
Original file line number Diff line number Diff line change
@@ -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 ***")

0 comments on commit 367b226

Please sign in to comment.