From b5c2665157106a8b88600759df5fb4d3b15ea4a2 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Thu, 14 Jun 2018 07:18:44 -0500 Subject: [PATCH 01/13] Added single distance calculations to lib.distances calc_distance, calc_angle and calc_dihedral Fixes Issue #1262 #1938 redid 180 degree angle test for new function --- package/CHANGELOG | 2 + .../analysis/hbonds/hbond_analysis.py | 12 ++--- package/MDAnalysis/core/topologyobjects.py | 33 ++++++------- package/MDAnalysis/lib/distances.py | 48 +++++++++++++++++++ .../core/test_topologyobjects.py | 8 ++-- 5 files changed, 75 insertions(+), 28 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index d54a313609a..ae564c3d795 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -51,6 +51,8 @@ Fixes * PDBWriter now properly sets start value Changes + * Added calc_distance, calc_angle and calc_dihedral to lib.distances + (Issue #1262 #1938) * TopologyAttrs are now statically typed (Issue #1876) * updated meta data for new PyPi (#1837) * AtomGroup.atoms, ResidueGroup.residues, and SegmentGroup.segments now return diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py index bf6b376479b..f7edc9d031e 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py @@ -328,7 +328,8 @@ class HydrogenBondAnalysis_OtherFF(HydrogenBondAnalysis): from MDAnalysis import MissingDataWarning, NoDataError, SelectionError, SelectionWarning from MDAnalysis.lib.log import ProgressMeter, _set_verbose from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch -from MDAnalysis.lib.distances import calc_bonds, calc_angles +from MDAnalysis.lib import distances +from MDAnalysis.lib.distances import calc_bonds logger = logging.getLogger('MDAnalysis.analysis.hbonds') @@ -1037,17 +1038,12 @@ def _get_timestep(): @staticmethod def calc_angle(d, h, a, box=None): """Calculate the angle (in degrees) between two atoms with H at apex.""" - v1= d.position - v2= h.position - v3= a.position - angle= calc_angles(v1[None, :], v2[None, :], v3[None, :], box=box) - - return np.rad2deg(angle[0]) + return distances.calc_angle(d.position, h.position, a.position, box) @staticmethod def calc_eucl_distance(a1, a2, box=None): """Calculate the Euclidean distance between two atoms. """ - return calc_bonds(a1.position[None, :], a2.position[None, :], box=box)[0] + return distances.calc_distance(a1.position, a2.position, box) @property def timeseries(self): diff --git a/package/MDAnalysis/core/topologyobjects.py b/package/MDAnalysis/core/topologyobjects.py index bd586ec0556..2ca6b111c2d 100644 --- a/package/MDAnalysis/core/topologyobjects.py +++ b/package/MDAnalysis/core/topologyobjects.py @@ -197,13 +197,9 @@ def length(self, pbc=False): .. versionchanged:: 0.11.0 Added pbc keyword """ - if pbc: - box = self.universe.dimensions - return distances.self_distance_array( - np.array([self[0].position, self[1].position]), - box=box)[0] - else: - return mdamath.norm(self[0].position - self[1].position) + box = self.universe.dimensions if pbc else None + + return distances.calc_distance(self[0].position, self[1].position, box) value = length @@ -220,7 +216,7 @@ class Angle(TopologyObject): """ btype = 'angle' - def angle(self): + def angle(self, pbc=False): """Returns the angle in degrees of this Angle. Angle between atoms 0 and 2 with apex at 1:: @@ -238,10 +234,13 @@ def angle(self): .. versionadded:: 0.9.0 .. versionchanged:: 0.17.0 Fixed angles close to 180 giving NaN + .. versionchanged:: 0.18.1 + Added pbc keyword """ - a = self[0].position - self[1].position - b = self[2].position - self[1].position - return np.rad2deg(mdamath.angle(a, b)) + box = self.universe.dimensions if pbc else None + + return distances.calc_angle( + self[0].position, self[1].position, self[2].position, box) value = angle @@ -265,7 +264,7 @@ class Dihedral(TopologyObject): # http://cbio.bmt.tue.nl/pumma/uploads/Theory/dihedral.png btype = 'dihedral' - def dihedral(self): + def dihedral(self, pbc=False): """Calculate the dihedral angle in degrees. Dihedral angle around axis connecting atoms 1 and 2 (i.e. the angle @@ -284,12 +283,14 @@ def dihedral(self): 4 decimals (and is only tested to 3 decimals). .. versionadded:: 0.9.0 + .. versionchanged:: 0.18.1 + Added pbc keyword """ + box = self.universe.dimensions if pbc else None A, B, C, D = self.atoms - ab = A.position - B.position - bc = B.position - C.position - cd = C.position - D.position - return np.rad2deg(mdamath.dihedral(ab, bc, cd)) + + return distances.calc_dihedral( + A.position, B.position, C.position, D.position, box) value = dihedral diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 01bbeb868d0..4890fbf0283 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -879,3 +879,51 @@ def apply_PBC(incoords, box, backend="serial"): backend=backend) return coords + + +def calc_distance(a, b, box=None): + """Distance between a and b + + Parameters + ---------- + a, b : numpy.ndarray + single coordinate vectors + box : numpy.ndarray, optional + simulation box, if given periodic boundary conditions will be applied + + .. versionadded:: 0.18.1 + """ + return calc_bonds(a[None, :], b[None, :], box=box)[0] + + +def calc_angle(a, b, c, box=None): + """Angle (in degrees) between a, b and c, where b is apex of angle + + Parameters + ---------- + a, b, c : numpy.ndarray + single coordinate vectors + box : numpy.ndarray, optional + simulation box, if given periodic boundary conditions will be applied + to the vectors a->b and b->c. + + .. versionadded:: 0.18.1 + """ + return np.rad2deg(calc_angles(a[None, :], b[None, :], c[None, :], box=box)[0]) + + +def calc_dihedral(a, b, c, d, box=None): + """Dihedral angle (in degrees) between planes (a, b, c) and (b, c, d) + + Parameters + ---------- + a, b, c, d : numpy.ndarray + single coordinate vectors + box : numpy.ndarray, optional + simulation box, if given periodic boundary conditions will be applied + to the vectors a->b and b->c. + + .. versionadded:: 0.18.1 + """ + return np.rad2deg( + calc_dihedrals(a[None, :], b[None, :], c[None, :], d[None, :], box)[0]) diff --git a/testsuite/MDAnalysisTests/core/test_topologyobjects.py b/testsuite/MDAnalysisTests/core/test_topologyobjects.py index 6a1f9adb708..244771bcfd0 100644 --- a/testsuite/MDAnalysisTests/core/test_topologyobjects.py +++ b/testsuite/MDAnalysisTests/core/test_topologyobjects.py @@ -163,14 +163,14 @@ def test_angle_180(self): # we edit the coordinates, so make our own universe u = mda.Universe(PSF, DCD) angle = u.atoms[210].angles[0] - coords = np.array([[ 83.37363434, 65.01402283, 35.03564835], - [ 82.28535461, 59.99816513, 34.94399261], - [ 81.19387817, 54.97631073, 34.85218048]], + coords = np.array([[1, 1, 1], + [2, 1, 1], + [3, 1, 1]], dtype=np.float32) angle.atoms.positions = coords - assert_almost_equal(angle.value(), -180.0, self.precision) + assert_almost_equal(angle.value(), 180.0, self.precision) # Dihedral class check def test_dihedral(self, PSFDCD): From a40b33ec2a8f5fceb2330c17caf3e71ba551a7bf Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 16 Jun 2018 18:29:26 -0500 Subject: [PATCH 02/13] tests for new functions in lib.distances --- .../MDAnalysisTests/lib/test_distances.py | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 9a6c5aff6c6..571f4893338 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -37,3 +37,92 @@ def test_transform_StoR_pass(coord_dtype): test_r = mda.lib.distances.transform_StoR(s, box) assert_equal(original_r, test_r) + + +# different boxlengths to shift a coordinate +shifts = [ + ((0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0)), # no shifting + ((1, 0, 0), (0, 1, 1), (0, 0, 1), (1, 1, 0)), # single box lengths + ((-1, 0, 1), (0, -1, 0), (1, 0, 1), (-1, -1, -1)), # negative single + ((4, 3, -2), (-2, 2, 2), (-5, 2, 2), (0, 2, 2)), # multiple boxlengths +] + + +@pytest.mark.parametrize('shift', shifts) +@pytest.mark.parametrize('periodic', [True, False]) +def test_calc_distance(shift, periodic): + box = np.array([10, 10, 10, 90., 90., 90.], dtype=np.float32) + + coords = np.array([[1, 1, 1], [3, 1, 1]], dtype=np.float32) + + shift1, shift2, _, _ = shift + + coords[0] += shift1 * box[:3] + coords[1] += shift2 * box[:3] + + box = box if periodic else None + result = mda.lib.distances.calc_distance(coords[0], coords[1], box) + + reference = 2.0 if periodic else np.linalg.norm(coords[0] - coords[1]) + + assert result == pytest.approx(reference) + + +@pytest.mark.parametrize('case', [ + # 90 degree angle + (np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1]], dtype=np.float32), 90.), + # straight line / 180. + (np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]], dtype=np.float32), 180.), +]) +@pytest.mark.parametrize('shift', shifts) +@pytest.mark.parametrize('periodic', [True, False]) +def test_calc_angle(case, shift, periodic): + def manual_angle(x, y, z): + return np.rad2deg(mda.lib.mdamath.angle(y - x, y - z)) + + box = np.array([10, 20, 30, 90., 90., 90.], dtype=np.float32) + (a, b, c), ref = case + + shift1, shift2, shift3, _ = shift + + a += shift1 * box[:3] + b += shift2 * box[:3] + c += shift3 * box[:3] + + box = box if periodic else None + result = mda.lib.distances.calc_angle(a, b, c, box) + + reference = ref if periodic else manual_angle(a, b, c) + + assert result == pytest.approx(reference, abs=1e-3) + + +@pytest.mark.parametrize('case', [ + # 0 degree angle (cis) + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 2, 1]], dtype=np.float32), 0.), + # 180 degree (trans) + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 0, 1]], dtype=np.float32), 180.), +]) +@pytest.mark.parametrize('shift', shifts) +@pytest.mark.parametrize('periodic', [True, False]) +def test_calc_dihedral(case, shift, periodic): + def manual_dihedral(a, b, c, d): + return np.rad2deg(mda.lib.mdamath.dihedral(b - a, c - b, d - c)) + + box = np.array([10., 10., 10., 90., 90., 90.], dtype=np.float32) + + (a, b, c, d), ref = case + + shift1, shift2, shift3, shift4 = shift + + a += shift1 * box[:3] + b += shift2 * box[:3] + c += shift3 * box[:3] + d += shift4 * box[:3] + + box = box if periodic else None + result = mda.lib.distances.calc_dihedral(a, b, c, d, box) + + reference = ref if periodic else manual_dihedral(a, b, c, d) + + assert abs(result) == pytest.approx(abs(reference), abs=1e-3) From 0da4f4fdc42dc6499b4d6a0b8fc9c3886983f788 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 16 Jun 2018 18:32:09 -0500 Subject: [PATCH 03/13] simplified docstrings --- package/MDAnalysis/lib/distances.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 4890fbf0283..e82455935f9 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -905,7 +905,6 @@ def calc_angle(a, b, c, box=None): single coordinate vectors box : numpy.ndarray, optional simulation box, if given periodic boundary conditions will be applied - to the vectors a->b and b->c. .. versionadded:: 0.18.1 """ @@ -921,7 +920,6 @@ def calc_dihedral(a, b, c, d, box=None): single coordinate vectors box : numpy.ndarray, optional simulation box, if given periodic boundary conditions will be applied - to the vectors a->b and b->c. .. versionadded:: 0.18.1 """ From 0931383da977b4b851169811c939433d7ee4ae09 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 17 Jun 2018 13:59:51 -0500 Subject: [PATCH 04/13] added more tests for single distance calculations --- testsuite/MDAnalysisTests/lib/test_distances.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 571f4893338..611b31c60d1 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -51,7 +51,7 @@ def test_transform_StoR_pass(coord_dtype): @pytest.mark.parametrize('shift', shifts) @pytest.mark.parametrize('periodic', [True, False]) def test_calc_distance(shift, periodic): - box = np.array([10, 10, 10, 90., 90., 90.], dtype=np.float32) + box = np.array([10, 20, 30, 90., 90., 90.], dtype=np.float32) coords = np.array([[1, 1, 1], [3, 1, 1]], dtype=np.float32) @@ -73,6 +73,8 @@ def test_calc_distance(shift, periodic): (np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1]], dtype=np.float32), 90.), # straight line / 180. (np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]], dtype=np.float32), 180.), + # 45 + (np.array([[1, 1, 1], [1, 2, 1], [2, 1, 1]], dtype=np.float32), 45.), ]) @pytest.mark.parametrize('shift', shifts) @pytest.mark.parametrize('periodic', [True, False]) @@ -102,6 +104,14 @@ def manual_angle(x, y, z): (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 2, 1]], dtype=np.float32), 0.), # 180 degree (trans) (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 0, 1]], dtype=np.float32), 180.), + # 90 degree + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 1, 2]], dtype=np.float32), 90.), + # other 90 degree + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 1, 0]], dtype=np.float32), 90.), + # 45 degree + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 2, 2]], dtype=np.float32), 45.), + # 135 + (np.array([[1, 2, 1], [1, 1, 1], [2, 1, 1], [2, 0, 2]], dtype=np.float32), 135.), ]) @pytest.mark.parametrize('shift', shifts) @pytest.mark.parametrize('periodic', [True, False]) From 7a11becb293eaa12127c702ea85e797c209634c2 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 17 Jun 2018 14:02:39 -0500 Subject: [PATCH 05/13] finalised changelog defaulted Bond.value pbc to True to match other objects --- package/CHANGELOG | 6 ++++-- package/MDAnalysis/core/topologyobjects.py | 4 +++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index ae564c3d795..b104e006aad 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -35,6 +35,10 @@ Enhancements * Calculations in *Group.center() are performed in double precision (#PR1936) * Functions in lib.distances accept coordinate arrays of arbitrary dtype (PR #1936) + * Added calc_distance, calc_angle and calc_dihedral to lib.distances + (Issue #1262 #1938) + * Added pbc kwarg to Bond/Angle/Dihedral/Improper object value method, + default True. (Issue #1938) Fixes * rewind in the SingleFrameReader now reads the frame from the file (Issue #1929) @@ -51,8 +55,6 @@ Fixes * PDBWriter now properly sets start value Changes - * Added calc_distance, calc_angle and calc_dihedral to lib.distances - (Issue #1262 #1938) * TopologyAttrs are now statically typed (Issue #1876) * updated meta data for new PyPi (#1837) * AtomGroup.atoms, ResidueGroup.residues, and SegmentGroup.segments now return diff --git a/package/MDAnalysis/core/topologyobjects.py b/package/MDAnalysis/core/topologyobjects.py index 2ca6b111c2d..f26936c4074 100644 --- a/package/MDAnalysis/core/topologyobjects.py +++ b/package/MDAnalysis/core/topologyobjects.py @@ -191,11 +191,13 @@ def partner(self, atom): else: raise ValueError("Unrecognised Atom") - def length(self, pbc=False): + def length(self, pbc=True): """Length of the bond. .. versionchanged:: 0.11.0 Added pbc keyword + .. versionchanged:: 0.18.1 + Changed default of pbc to True """ box = self.universe.dimensions if pbc else None From 0f005f2a0d0d56a966ede42b8d44febcb252911d Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 17 Jun 2018 14:04:06 -0500 Subject: [PATCH 06/13] changed default to pbc to True in TopologyObjects --- package/MDAnalysis/core/topologyobjects.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyobjects.py b/package/MDAnalysis/core/topologyobjects.py index f26936c4074..6abfb02eb33 100644 --- a/package/MDAnalysis/core/topologyobjects.py +++ b/package/MDAnalysis/core/topologyobjects.py @@ -218,7 +218,7 @@ class Angle(TopologyObject): """ btype = 'angle' - def angle(self, pbc=False): + def angle(self, pbc=True): """Returns the angle in degrees of this Angle. Angle between atoms 0 and 2 with apex at 1:: @@ -266,7 +266,7 @@ class Dihedral(TopologyObject): # http://cbio.bmt.tue.nl/pumma/uploads/Theory/dihedral.png btype = 'dihedral' - def dihedral(self, pbc=False): + def dihedral(self, pbc=True): """Calculate the dihedral angle in degrees. Dihedral angle around axis connecting atoms 1 and 2 (i.e. the angle From 5353f95a4ed918058e9084714db3d9936a4752cc Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Tue, 19 Jun 2018 11:47:42 -0500 Subject: [PATCH 07/13] fixed make_whole now Bond.length() defaults to pbc=True --- package/MDAnalysis/core/topologyobjects.py | 4 ++-- package/MDAnalysis/lib/mdamath.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/core/topologyobjects.py b/package/MDAnalysis/core/topologyobjects.py index 6abfb02eb33..0844efa987e 100644 --- a/package/MDAnalysis/core/topologyobjects.py +++ b/package/MDAnalysis/core/topologyobjects.py @@ -237,7 +237,7 @@ def angle(self, pbc=True): .. versionchanged:: 0.17.0 Fixed angles close to 180 giving NaN .. versionchanged:: 0.18.1 - Added pbc keyword + Added pbc keyword, default True """ box = self.universe.dimensions if pbc else None @@ -286,7 +286,7 @@ def dihedral(self, pbc=True): .. versionadded:: 0.9.0 .. versionchanged:: 0.18.1 - Added pbc keyword + Added pbc keyword, default True """ box = self.universe.dimensions if pbc else None A, B, C, D = self.atoms diff --git a/package/MDAnalysis/lib/mdamath.py b/package/MDAnalysis/lib/mdamath.py index eccc43fdf55..dd58ea12bab 100644 --- a/package/MDAnalysis/lib/mdamath.py +++ b/package/MDAnalysis/lib/mdamath.py @@ -416,7 +416,7 @@ def make_whole(atomgroup, reference_atom=None): continue if other in ref_points: continue - if b.length() > box_length: + if b.length(pbc=False) > box_length: # Vector from ref atom to other vec = other.position - ref.position # Apply pbc to this vector From 61157394f803a73370bf01b4292a6182cb0ec500 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Tue, 19 Jun 2018 11:54:13 -0500 Subject: [PATCH 08/13] completely removed HBondAnalysis' distance methods --- .../analysis/hbonds/hbond_analysis.py | 22 ++++++---------- .../analysis/hbonds/wbridge_analysis.py | 25 +++++++++++++------ 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py index f7edc9d031e..42b3d9f37f4 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py @@ -807,7 +807,7 @@ def _get_bonded_hydrogens_list(self, atom, **kwargs): hydrogens = [ a for a in self.u.atoms[atom.index + 1:atom.index + 4] if (a.name.startswith(('H', '1H', '2H', '3H')) and - self.calc_eucl_distance(atom, a, box) < self.r_cov[atom.name[0]]) + distances.calc_distance(atom.position, a.position, box) < self.r_cov[atom.name[0]]) ] except IndexError: hydrogens = [] # weird corner case that atom is the last one in universe @@ -993,9 +993,10 @@ def _get_timestep(): for h in donor_h_set: res = ns_acceptors.search(h, self.distance) for a in res: - angle = self.calc_angle(d, h, a, box=box) + angle = distances.calc_angle(d.position, h.position, + a.position, box=box) donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a, box) + dist = distances.calc_distance(donor_atom.position, a.position, box) if angle >= self.angle and dist <= self.distance: self.logger_debug( "S1-D: {0!s} <-> S2-A: {1!s} {2:f} A, {3:f} DEG".format(h.index, a.index, dist, angle)) @@ -1018,9 +1019,10 @@ def _get_timestep(): (h.index, a.index) in already_found or (a.index, h.index) in already_found): continue - angle = self.calc_angle(d, h, a, box=box) + angle = distances.calc_angle(d.position, h.position, + a.position, box=box) donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a, box) + dist = distances.calc_distance(donor_atom.position, a.position, box) if angle >= self.angle and dist <= self.distance: self.logger_debug( "S1-A: {0!s} <-> S2-D: {1!s} {2:f} A, {3:f} DEG".format(a.index, h.index, dist, angle)) @@ -1035,16 +1037,6 @@ def _get_timestep(): logger.info("HBond analysis: complete; timeseries %s.timeseries", self.__class__.__name__) - @staticmethod - def calc_angle(d, h, a, box=None): - """Calculate the angle (in degrees) between two atoms with H at apex.""" - return distances.calc_angle(d.position, h.position, a.position, box) - - @staticmethod - def calc_eucl_distance(a1, a2, box=None): - """Calculate the Euclidean distance between two atoms. """ - return distances.calc_distance(a1.position, a2.position, box) - @property def timeseries(self): """Time series of hydrogen bonds. diff --git a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py index 308f7d077a4..e55272c7a2a 100644 --- a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py @@ -292,6 +292,7 @@ class WaterBridgeAnalysis_OtherFF(WaterBridgeAnalysis): from .hbond_analysis import HydrogenBondAnalysis from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch from MDAnalysis.lib.log import ProgressMeter, _set_verbose +from MDAnalysis.lib import distances from MDAnalysis import SelectionWarning logger = logging.getLogger('MDAnalysis.analysis.wbridges') @@ -634,9 +635,11 @@ def run(self, **kwargs): res = ns_acceptors.search(h, self.distance) for a in res: donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a) + dist = distances.calc_distance(donor_atom.position, + a.position) if dist <= self.distance: - angle = self.calc_angle(d, h, a) + angle = distances.calc_angle(d.position, h.position, + a.position) if angle >= self.angle: self.logger_debug( "S1-D: {0!s} <-> W-A: {1!s} {2:f} A, {3:f} DEG"\ @@ -658,9 +661,11 @@ def run(self, **kwargs): res = ns_acceptors.search(h, self.distance) for a in res: donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a) + dist = distances.calc_distance(donor_atom.position, + a.position) if dist <= self.distance: - angle = self.calc_angle(d, h, a) + angle = distances.calc_angle(d.position, h.position, + a.position) if angle >= self.angle: self.logger_debug( "S1-A: {0!s} <-> W-D: {1!s} {2:f} A, {3:f} DEG"\ @@ -706,9 +711,11 @@ def run(self, **kwargs): res = ns_acceptors.search(h, self.distance) for a in res: donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a) + dist = distances.calc_distance(donor_atom.position, + a.position) if dist <= self.distance: - angle = self.calc_angle(d, h, a) + angle = distances.calc_angle(d.position, h.position, + a.position) if angle >= self.angle: self.logger_debug( "WB-D: {0!s} <-> S2-A: {1!s} {2:f} A, {3:f} DEG"\ @@ -728,9 +735,11 @@ def run(self, **kwargs): res = ns_acceptors.search(h, self.distance) for a in res: donor_atom = h if self.distance_type != 'heavy' else d - dist = self.calc_eucl_distance(donor_atom, a) + dist = distances.calc_distance(donor_atom.position, + a.position) if dist <= self.distance: - angle = self.calc_angle(d, h, a) + angle = distances.calc_angle(d.position, h.position, + a.position) if angle >= self.angle: self.logger_debug( "WB-A: {0!s} <-> S2-D: {1!s} {2:f} A, {3:f} DEG"\ From be1ce8cfaa96dfa457efd556d5a4c7b6a5c11b45 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Wed, 20 Jun 2018 09:32:28 -0500 Subject: [PATCH 09/13] removed unused import --- package/MDAnalysis/analysis/hbonds/hbond_analysis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py index 42b3d9f37f4..ac38eceb8a9 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py @@ -329,7 +329,6 @@ class HydrogenBondAnalysis_OtherFF(HydrogenBondAnalysis): from MDAnalysis.lib.log import ProgressMeter, _set_verbose from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch from MDAnalysis.lib import distances -from MDAnalysis.lib.distances import calc_bonds logger = logging.getLogger('MDAnalysis.analysis.hbonds') From c4fc8ec6d423668d493f79e61ea0421ec6c577e0 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 23 Jun 2018 19:11:04 -0500 Subject: [PATCH 10/13] maybe fixed the docs? --- package/MDAnalysis/lib/distances.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index e82455935f9..d25f5cd8aac 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -887,9 +887,9 @@ def calc_distance(a, b, box=None): Parameters ---------- a, b : numpy.ndarray - single coordinate vectors + single coordinate vectors box : numpy.ndarray, optional - simulation box, if given periodic boundary conditions will be applied + simulation box, if given periodic boundary conditions will be applied .. versionadded:: 0.18.1 """ @@ -902,9 +902,10 @@ def calc_angle(a, b, c, box=None): Parameters ---------- a, b, c : numpy.ndarray - single coordinate vectors - box : numpy.ndarray, optional - simulation box, if given periodic boundary conditions will be applied + single coordinate vectors + box : numpy.ndarray + simulation box if given periodic boundary conditions will be applied to + the vectors between atoms .. versionadded:: 0.18.1 """ @@ -917,9 +918,9 @@ def calc_dihedral(a, b, c, d, box=None): Parameters ---------- a, b, c, d : numpy.ndarray - single coordinate vectors + single coordinate vectors box : numpy.ndarray, optional - simulation box, if given periodic boundary conditions will be applied + simulation box, if given periodic boundary conditions will be applied .. versionadded:: 0.18.1 """ From 068934372d42048bc4e92d903fa434cfc89ee86b Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 24 Jun 2018 10:51:59 -0500 Subject: [PATCH 11/13] really really fixed the docs --- package/MDAnalysis/lib/distances.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index d25f5cd8aac..0411a6d4ab0 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -891,6 +891,7 @@ def calc_distance(a, b, box=None): box : numpy.ndarray, optional simulation box, if given periodic boundary conditions will be applied + .. versionadded:: 0.18.1 """ return calc_bonds(a[None, :], b[None, :], box=box)[0] @@ -907,6 +908,7 @@ def calc_angle(a, b, c, box=None): simulation box if given periodic boundary conditions will be applied to the vectors between atoms + .. versionadded:: 0.18.1 """ return np.rad2deg(calc_angles(a[None, :], b[None, :], c[None, :], box=box)[0]) @@ -922,6 +924,7 @@ def calc_dihedral(a, b, c, d, box=None): box : numpy.ndarray, optional simulation box, if given periodic boundary conditions will be applied + .. versionadded:: 0.18.1 """ return np.rad2deg( From f3cbffdc8bee42a321c6142a5680554268c6541c Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 24 Jun 2018 11:03:23 -0500 Subject: [PATCH 12/13] more doc fixes in lib.distances --- package/MDAnalysis/lib/distances.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 0411a6d4ab0..773eed09336 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -413,7 +413,7 @@ def transform_RtoS(inputcoords, box, backend="serial"): The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - backend : str + backend : str, optional Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). @@ -422,6 +422,7 @@ def transform_RtoS(inputcoords, box, backend="serial"): outcoords : array A n x 3 array of fractional coordiantes. + .. versionchanged:: 0.13.0 Added *backend* keyword. .. versionchanged:: 0.19.0 @@ -475,14 +476,14 @@ def transform_StoR(inputcoords, box, backend="serial"): The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - backend : str + backend : str, optional Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). Returns ------- outcoords : array - A n x 3 array of fracional coordiantes. + A n x 3 array of fractional coordiantes. .. versionchanged:: 0.13.0 From 68446bd398ab83a1c3348cc94e2e039f63465ab8 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Fri, 29 Jun 2018 10:12:51 -0500 Subject: [PATCH 13/13] use np.testing.assert_almost_equal not pytest.approx --- testsuite/MDAnalysisTests/lib/test_distances.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 611b31c60d1..61857fffe73 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -23,7 +23,7 @@ from __future__ import absolute_import import pytest import numpy as np -from numpy.testing import assert_equal +from numpy.testing import assert_equal, assert_almost_equal import MDAnalysis as mda @@ -65,7 +65,7 @@ def test_calc_distance(shift, periodic): reference = 2.0 if periodic else np.linalg.norm(coords[0] - coords[1]) - assert result == pytest.approx(reference) + assert_almost_equal(result, reference, decimal=3) @pytest.mark.parametrize('case', [ @@ -96,7 +96,7 @@ def manual_angle(x, y, z): reference = ref if periodic else manual_angle(a, b, c) - assert result == pytest.approx(reference, abs=1e-3) + assert_almost_equal(result, reference, decimal=3) @pytest.mark.parametrize('case', [ @@ -135,4 +135,4 @@ def manual_dihedral(a, b, c, d): reference = ref if periodic else manual_dihedral(a, b, c, d) - assert abs(result) == pytest.approx(abs(reference), abs=1e-3) + assert_almost_equal(abs(result), abs(reference), decimal=3)