added single distance calculations to lib.distances

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)

package/MDAnalysis/analysis/hbonds/hbond_analysis.py

@@ -328,7 +328,7 @@ 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
 
 
 logger = logging.getLogger('MDAnalysis.analysis.hbonds')
@@ -806,7 +806,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
@@ -992,9 +992,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))
@@ -1017,9 +1018,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))
@@ -1034,21 +1036,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."""
-        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])
-
-    @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]
-
     @property
     def timeseries(self):
         """Time series of hydrogen bonds.

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"\

package/MDAnalysis/core/topologyobjects.py

@@ -191,19 +191,17 @@ 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
         """
-        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 +218,7 @@ class Angle(TopologyObject):
     """
     btype = 'angle'
 
-    def angle(self):
+    def angle(self, pbc=True):
         """Returns the angle in degrees of this Angle.
 
         Angle between atoms 0 and 2 with apex at 1::
@@ -238,10 +236,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, default True
         """
-        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 +266,7 @@ class Dihedral(TopologyObject):
     # http://cbio.bmt.tue.nl/pumma/uploads/Theory/dihedral.png
     btype = 'dihedral'
 
-    def dihedral(self):
+    def dihedral(self, pbc=True):
         """Calculate the dihedral angle in degrees.
 
         Dihedral angle around axis connecting atoms 1 and 2 (i.e. the angle
@@ -284,12 +285,14 @@ def dihedral(self):
         4 decimals (and is only tested to 3 decimals).
 
         .. versionadded:: 0.9.0
+        .. versionchanged:: 0.18.1
+           Added pbc keyword, default True
         """
+        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
 

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
@@ -879,3 +880,53 @@ 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
+        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])
+
+
+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
+
+
+    .. versionadded:: 0.18.1
+    """
+    return np.rad2deg(
+        calc_dihedrals(a[None, :], b[None, :], c[None, :], d[None, :], box)[0])

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

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):