diff --git a/package/CHANGELOG b/package/CHANGELOG index d54a313609a..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) diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py index bf6b376479b..ac38eceb8a9 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ b/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. 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"\ diff --git a/package/MDAnalysis/core/topologyobjects.py b/package/MDAnalysis/core/topologyobjects.py index bd586ec0556..0844efa987e 100644 --- a/package/MDAnalysis/core/topologyobjects.py +++ b/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 diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 01bbeb868d0..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 @@ -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]) 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 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): diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 9a6c5aff6c6..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 @@ -37,3 +37,102 @@ 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, 20, 30, 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_almost_equal(result, reference, decimal=3) + + +@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.), + # 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]) +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_almost_equal(result, reference, decimal=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.), + # 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]) +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_almost_equal(abs(result), abs(reference), decimal=3)