Skip to content

Commit

Permalink
Added single distance calculations to lib.distances
Browse files Browse the repository at this point in the history
calc_distance, calc_angle and calc_dihedral

Fixes Issue #1262 #1938

redid 180 degree angle test for new function
  • Loading branch information
richardjgowers committed Jun 29, 2018
1 parent 339c510 commit b5c2665
Show file tree
Hide file tree
Showing 5 changed files with 75 additions and 28 deletions.
2 changes: 2 additions & 0 deletions package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 4 additions & 8 deletions package/MDAnalysis/analysis/hbonds/hbond_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand Down Expand Up @@ -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):
Expand Down
33 changes: 17 additions & 16 deletions package/MDAnalysis/core/topologyobjects.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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::
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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

Expand Down
48 changes: 48 additions & 0 deletions package/MDAnalysis/lib/distances.py
Original file line number Diff line number Diff line change
Expand Up @@ -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])
8 changes: 4 additions & 4 deletions testsuite/MDAnalysisTests/core/test_topologyobjects.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down

0 comments on commit b5c2665

Please sign in to comment.