diff --git a/cmeutils/geometry.py b/cmeutils/geometry.py index 167e031..61949dc 100644 --- a/cmeutils/geometry.py +++ b/cmeutils/geometry.py @@ -60,8 +60,8 @@ def get_plane_normal(points): return ctr, normal -def angle_between_vectors(u, v, min_angle=True): - """Calculate the angle between two vectors in degrees. +def angle_between_vectors(u, v, min_angle=True, degrees=True): + """Calculate the angle between two vectors in radians or degrees. Parameters ---------- @@ -73,15 +73,21 @@ def angle_between_vectors(u, v, min_angle=True): Whether to return the supplement if the angle is greater than 90 degrees. Useful for calculating the minimum angle between the normal vectors of planes as direction doesn't matter. + degrees : bool, default True + If True, the angle is returned in degrees. + If False, the angle is returned in radians. Returns ------- angle: float - Angle between u and v in degrees + Angle between u and v + """ - angle = np.rad2deg( - np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v))) - ) - if angle > 90 and min_angle: - return 180 - angle + # Angle in radians + angle = np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v))) + if angle > np.pi/2 and min_angle: + angle = np.pi - angle + + if degrees: + return np.rad2deg(angle) return angle diff --git a/cmeutils/plotting.py b/cmeutils/plotting.py index e0232b7..1f9dc90 100644 --- a/cmeutils/plotting.py +++ b/cmeutils/plotting.py @@ -1,7 +1,8 @@ import matplotlib.pyplot as plt import numpy as np -def get_histogram(data, normalize=False, bins="auto"): + +def get_histogram(data, normalize=False, bins="auto", x_range=None): """Bins a 1-D array of data into a histogram using the numpy.histogram method. @@ -16,6 +17,9 @@ def get_histogram(data, normalize=False, bins="auto"): bins : float, int, or str, default="auto" Method used by numpy to determine bin borders. Check the numpy.histogram docs for more details. + x_range : (float, float), default = None + The lower and upper range of the histogram bins. + If set to None, then the min and max values of data are used. Returns ------- @@ -25,7 +29,7 @@ def get_histogram(data, normalize=False, bins="auto"): Array of the bin height values """ - bin_heights, bin_borders = np.histogram(data, bins=bins) + bin_heights, bin_borders = np.histogram(data, bins=bins, range=x_range) if normalize is True: bin_heights = bin_heights/sum(bin_heights) bin_widths = np.diff(bin_borders) diff --git a/cmeutils/structure.py b/cmeutils/structure.py index 9f9a3a9..fef9e12 100644 --- a/cmeutils/structure.py +++ b/cmeutils/structure.py @@ -1,3 +1,5 @@ +import warnings + import freud import gsd import gsd.hoomd @@ -16,7 +18,10 @@ def angle_distribution( C_name, start=0, stop=-1, + degrees=False, histogram=False, + theta_min=0.0, + theta_max=None, normalize=False, bins="auto" ): @@ -35,10 +40,19 @@ def angle_distribution( Negative numbers index from the end. (default 0) stop : int Final frame index for accumulating bond lengths. (default -1) + degrees : bool, default=False + If True, the angle values are returned in degrees. + if False, the angle values are returned in radians. histogram : bool, default=False If set to True, places the resulting angles into a histogram and retrums the histogram's bin centers and heights as opposed to the actual calcualted angles. + theta_min : float, default = 0.0 + Sets the minimum theta value to be included in the distribution + theta_max : float, default = None + Sets the maximum theta value to be included in the distribution + If left as None, then theta_max will be either pi radians or + 180 degrees depending on the value set for the degrees parameter normalize : bool, default=False If set to True, normalizes the angle distribution by the sum of the bin heights, so that the distribution adds up to 1. @@ -55,7 +69,11 @@ def angle_distribution( If histogram is True, returns a 2D array of bin centers and bin heights. """ - angles = [] + if not degrees and theta_max is None: + theta_max = np.pi + elif degrees and theta_max is None: + theta_max = 180 + trajectory = gsd.hoomd.open(gsd_file, mode="rb") name = "-".join([A_name, B_name, C_name]) name_rev = "-".join([C_name, B_name, A_name]) @@ -81,13 +99,23 @@ def angle_distribution( pos1_unwrap = pos1 + (img1 * snap.configuration.box[:3]) pos2_unwrap = pos2 + (img2 * snap.configuration.box[:3]) pos3_unwrap = pos3 + (img3 * snap.configuration.box[:3]) - u = pos2_unwrap - pos1_unwrap + u = pos1_unwrap - pos2_unwrap v = pos3_unwrap - pos2_unwrap - angles.append(np.round(angle_between_vectors(u, v, False), 3)) + angles.append( + np.round(angle_between_vectors(u, v, False, degrees), 3) + ) trajectory.close() if histogram: - bin_centers, bin_heights = get_histogram(np.array(angles), bins=bins) + if angles[0] < theta_min or angles[-1] > theta_max: + warnings.warn("There are bond angles that fall outside of " + "your set theta_min and theta_max range. " + "You may want to adjust this range to " + "include all bond angles." + ) + bin_centers, bin_heights = get_histogram( + np.array(angles), bins=bins, x_range=(theta_min, theta_max) + ) return np.stack((bin_centers, bin_heights)).T else: return np.array(angles) @@ -100,6 +128,8 @@ def bond_distribution( start=0, stop=-1, histogram=False, + l_min=0.0, + l_max=4.0, normalize=True, bins=100 ): @@ -121,6 +151,10 @@ def bond_distribution( If set to True, places the resulting bonds into a histogram and retrums the histogram's bin centers and heights as opposed to the actual calcualted bonds. + l_min : float, default = 0.0 + Sets the minimum bond length to be included in the distribution + l_max : float, default = 5.0 + Sets the maximum bond length value to be included in the distribution normalize : bool, default=False If set to True, normalizes the angle distribution by the sum of the bin heights, so that the distribution adds up to 1. @@ -162,7 +196,14 @@ def bond_distribution( trajectory.close() if histogram: - bin_centers, bin_heights = get_histogram(np.array(bonds), bins=bins) + if bonds[0] < l_min or bonds[-1] > l_max: + warnings.warn("There are bond lengths that fall outside of " + "your set l_min and l_max range. You may want to adjust " + "this range to include all bond lengths." + ) + bin_centers, bin_heights = get_histogram( + np.array(bonds), bins=bins, x_range=(l_min, l_max) + ) return np.stack((bin_centers, bin_heights)).T else: return np.array(bonds) diff --git a/cmeutils/tests/test_geometry.py b/cmeutils/tests/test_geometry.py index c0d3857..786a897 100644 --- a/cmeutils/tests/test_geometry.py +++ b/cmeutils/tests/test_geometry.py @@ -1,3 +1,5 @@ +import math + import numpy as np import pytest @@ -20,10 +22,25 @@ def test_get_plane_normal(self): with pytest.raises(AssertionError): get_plane_normal(points[:2]) - def test_angle_between_vectors(self): + def test_angle_between_vectors_deg(self): assert np.isclose( 90, angle_between_vectors(np.array([1,0,0]),np.array([0,1,0])) ) assert np.isclose( 0, angle_between_vectors(np.array([1,0,0]),np.array([-1,0,0])) ) + + def test_angle_between_vectors_rad(self): + assert np.isclose( + math.pi/2, + angle_between_vectors( + np.array([1,0,0]),np.array([0,1,0]), degrees=False + ) + ) + assert np.isclose( + 0, + angle_between_vectors( + np.array([1,0,0]),np.array([-1,0,0]), degrees=False + ) + ) + diff --git a/cmeutils/tests/test_structure.py b/cmeutils/tests/test_structure.py index 01d46de..6a7118e 100644 --- a/cmeutils/tests/test_structure.py +++ b/cmeutils/tests/test_structure.py @@ -18,11 +18,32 @@ ) class TestStructure(BaseTest): - def test_angle_distribution(self, p3ht_gsd): - angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1) + def test_angle_distribution_deg(self, p3ht_gsd): + angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, degrees=True) for ang in angles: assert 80 < ang < 100 + def test_angle_distribution_range(self, p3ht_gsd): + angles = angle_distribution( + p3ht_gsd, + "cc", + "ss", + "cc", + start=0, + stop=1, + histogram=True, + degrees=True, + theta_min=10, + theta_max=180 + ) + assert np.allclose(angles[0,0],10, atol=0.5) + assert np.allclose(angles[-1,0], 180, atol=0.5) + + def test_angle_distribution_rad(self, p3ht_gsd): + angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, degrees=False) + for ang in angles: + assert 1.40 < ang < 1.75 + def test_angle_distribution_order(self, p3ht_gsd): angles = angle_distribution(p3ht_gsd, "ss", "cc", "cd", start=0, stop=1) angles2 = angle_distribution(p3ht_gsd, "cd", "cc", "ss", start=0, stop=1) @@ -37,8 +58,23 @@ def test_angle_not_found(self, p3ht_gsd): def test_bond_distribution(self, p3ht_gsd): bonds = bond_distribution(p3ht_gsd, "cc", "ss", start=0, stop=1) for bond in bonds: + print(bond) assert 0.45 < bond < 0.52 + def test_bond_distribution_range(self, p3ht_gsd): + bonds = bond_distribution( + p3ht_gsd, + "cc", + "ss", + start=0, + stop=1, + l_min=0, + l_max=1, + histogram=True + ) + assert np.allclose(bonds[0,0],0, atol=0.5) + assert np.allclose(bonds[-1,0], 1, atol=0.5) + def test_bond_distribution_order(self, p3ht_gsd): bonds = bond_distribution(p3ht_gsd, "cc", "ss", start=0, stop=1) bonds2 = bond_distribution(p3ht_gsd, "ss", "cc", start=0, stop=1) @@ -59,6 +95,19 @@ def test_bond_histogram(self, p3ht_gsd): assert bonds_hist.ndim == 2 assert bonds_no_hist.ndim == 1 + def test_bond_range_outside(self, p3ht_gsd): + with pytest.warns(UserWarning): + bonds_hist = bond_distribution( + p3ht_gsd, + "cc", + "ss", + start=0, + stop=1, + histogram=True, + l_min=0.52, + l_max=0.60 + ) + def test_angle_histogram(self, p3ht_gsd): angles_hist = angle_distribution( p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, histogram=True @@ -68,6 +117,20 @@ def test_angle_histogram(self, p3ht_gsd): ) assert angles_hist.ndim == 2 assert angles_no_hist.ndim == 1 + + def test_angle_range_outside(self, p3ht_gsd): + with pytest.warns(UserWarning): + angles_hist = angle_distribution( + p3ht_gsd, + "cc", + "ss", + "cc", + start=0, + stop=1, + histogram=True, + theta_min = 120, + theta_max=180 + ) def test_gsd_rdf(self, gsdfile_bond): rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "B")