From b5c2665157106a8b88600759df5fb4d3b15ea4a2 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Thu, 14 Jun 2018 07:18:44 -0500 Subject: [PATCH] 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):