1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
3#
4# MDAnalysis --- https://www.mdanalysis.org
5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
6# (see the file AUTHORS for the full list of names)
7#
8# Released under the GNU Public Licence, v2 or any higher version
9#
10# Please cite your use of MDAnalysis in published work:
11#
12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
17#
18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
21#
22
23"""
24Mathematical helper functions --- :mod:`MDAnalysis.lib.mdamath`
25===============================================================
26
27Helper functions for common mathematical operations
28
29.. autofunction:: normal
30.. autofunction:: norm
31.. autofunction:: angle
32.. autofunction:: dihedral
33.. autofunction:: stp
34.. autofunction:: triclinic_box
35.. autofunction:: triclinic_vectors
36.. autofunction:: box_volume
37.. autofunction:: make_whole
38.. autofunction:: find_fragments
39
40.. versionadded:: 0.11.0
41"""
42from __future__ import division, absolute_import
43from six.moves import zip
44import numpy as np
45
46from ..exceptions import NoDataError
47from . import util
48from ._cutil import make_whole, find_fragments
49
50# geometric functions
51def norm(v):
52    r"""Calculate the norm of a vector v.
53
54    .. math:: v = \sqrt{\mathbf{v}\cdot\mathbf{v}}
55
56    This version is faster then numpy.linalg.norm because it only works for a
57    single vector and therefore can skip a lot of the additional fuss
58    linalg.norm does.
59
60    Parameters
61    ----------
62    v : array_like
63        1D array of shape (N) for a vector of length N
64
65    Returns
66    -------
67    float
68        norm of the vector
69
70    """
71    return np.sqrt(np.dot(v, v))
72
73
74def normal(vec1, vec2):
75    r"""Returns the unit vector normal to two vectors.
76
77    .. math::
78
79       \hat{\mathbf{n}} = \frac{\mathbf{v}_1 \times \mathbf{v}_2}{|\mathbf{v}_1 \times \mathbf{v}_2|}
80
81    If the two vectors are collinear, the vector :math:`\mathbf{0}` is returned.
82
83    .. versionchanged:: 0.11.0
84       Moved into lib.mdamath
85    """
86    normal = np.cross(vec1, vec2)
87    n = norm(normal)
88    if n == 0.0:
89        return normal  # returns [0,0,0] instead of [nan,nan,nan]
90    return normal / n  # ... could also use numpy.nan_to_num(normal/norm(normal))
91
92
93def angle(a, b):
94    """Returns the angle between two vectors in radians
95
96    .. versionchanged:: 0.11.0
97       Moved into lib.mdamath
98    """
99    x = np.dot(a, b) / (norm(a) * norm(b))
100    # catch roundoffs that lead to nan otherwise
101    if x > 1.0:
102        return 0.0
103    elif x < -1.0:
104        return -np.pi
105    return np.arccos(x)
106
107
108def stp(vec1, vec2, vec3):
109    r"""Takes the scalar triple product of three vectors.
110
111    Returns the volume *V* of the parallel epiped spanned by the three
112    vectors
113
114    .. math::
115
116        V = \mathbf{v}_3 \cdot (\mathbf{v}_1 \times \mathbf{v}_2)
117
118    .. versionchanged:: 0.11.0
119       Moved into lib.mdamath
120    """
121    return np.dot(vec3, np.cross(vec1, vec2))
122
123
124def dihedral(ab, bc, cd):
125    r"""Returns the dihedral angle in radians between vectors connecting A,B,C,D.
126
127    The dihedral measures the rotation around bc::
128
129         ab
130       A---->B
131              \ bc
132              _\'
133                C---->D
134                  cd
135
136    The dihedral angle is restricted to the range -π <= x <= π.
137
138    .. versionadded:: 0.8
139    .. versionchanged:: 0.11.0
140       Moved into lib.mdamath
141    """
142    x = angle(normal(ab, bc), normal(bc, cd))
143    return (x if stp(ab, bc, cd) <= 0.0 else -x)
144
145
146def _angle(a, b):
147    """Angle between two vectors *a* and *b* in degrees.
148
149    If one of the lengths is 0 then the angle is returned as 0
150    (instead of `nan`).
151    """
152    # This function has different limits than angle?
153
154    angle = np.arccos(np.dot(a, b) / (norm(a) * norm(b)))
155    if np.isnan(angle):
156        return 0.0
157    return np.rad2deg(angle)
158
159
160def triclinic_box(x, y, z):
161    """Convert the three triclinic box vectors to [A,B,C,alpha,beta,gamma].
162
163    Angles are in degrees.
164
165    * alpha  = angle(y,z)
166    * beta   = angle(x,z)
167    * gamma  = angle(x,y)
168
169    Note
170    ----
171    Definition of angles: http://en.wikipedia.org/wiki/Lattice_constant
172    """
173    A, B, C = [norm(v) for v in (x, y, z)]
174    alpha = _angle(y, z)
175    beta = _angle(x, z)
176    gamma = _angle(x, y)
177    return np.array([A, B, C, alpha, beta, gamma], dtype=np.float32)
178
179
180def triclinic_vectors(dimensions):
181    """Convert `[A,B,C,alpha,beta,gamma]` to a triclinic box representation.
182
183    Original `code by Tsjerk Wassenaar`_ posted on the Gromacs mailinglist.
184
185    If *dimensions* indicates a non-periodic system (i.e. all lengths
186    0) then null-vectors are returned.
187
188    .. _code by Tsjerk Wassenaar:
189       http://www.mail-archive.com/gmx-users@gromacs.org/msg28032.html
190
191    Parameters
192    ----------
193    dimensions : [A, B, C, alpha, beta, gamma]
194        list of box lengths and angles (in degrees) such as
195        ``[A,B,C,alpha,beta,gamma]``
196
197    Returns
198    -------
199    numpy.array
200        numpy 3x3 array B, with ``B[0]`` as first box vector,
201        ``B[1]``as second vector, ``B[2]`` as third box vector.
202
203    Notes
204    -----
205    The first vector is always pointing along the X-axis,
206    i.e., parallel to (1, 0, 0).
207
208
209    .. versionchanged:: 0.7.6
210       Null-vectors are returned for non-periodic (or missing) unit cell.
211
212    """
213    B = np.zeros((3, 3), dtype=np.float32)
214    x, y, z, a, b, c = dimensions[:6]
215
216    if np.all(dimensions[:3] == 0):
217        return B
218
219    B[0][0] = x
220    if a == 90. and b == 90. and c == 90.:
221        B[1][1] = y
222        B[2][2] = z
223    else:
224        a = np.deg2rad(a)
225        b = np.deg2rad(b)
226        c = np.deg2rad(c)
227        B[1][0] = y * np.cos(c)
228        B[1][1] = y * np.sin(c)
229        B[2][0] = z * np.cos(b)
230        B[2][1] = z * (np.cos(a) - np.cos(b) * np.cos(c)) / np.sin(c)
231        B[2][2] = np.sqrt(z * z - B[2][0] ** 2 - B[2][1] ** 2)
232    return B
233
234
235def box_volume(dimensions):
236    """Return the volume of the unitcell described by *dimensions*.
237
238    The volume is computed as `det(x1,x2,x2)` where the xi are the
239    triclinic box vectors from :func:`triclinic_vectors`.
240
241    Parameters
242    ----------
243    dimensions : [A, B, C, alpha, beta, gamma]
244        list of box lengths and angles (in degrees) such as
245        [A,B,C,alpha,beta,gamma]
246
247    Returns
248    -------
249    volume : float
250    """
251    return np.linalg.det(triclinic_vectors(dimensions))
252
253