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