Source code for mdtraj.utils.unitcell

##############################################################################
# MDTraj: A Python Library for Loading, Saving, and Manipulating
#         Molecular Dynamics Trajectories.
# Copyright 2012-2013 Stanford University and the Authors
#
# Authors: Robert McGibbon
# Contributors:
#
# MDTraj is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 2.1
# of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with MDTraj. If not, see <http://www.gnu.org/licenses/>.
##############################################################################


import warnings

import numpy as np

##############################################################################
# Functions
##############################################################################


[docs] def lengths_and_angles_to_box_vectors(a_length, b_length, c_length, alpha, beta, gamma): """Convert from the lengths/angles of the unit cell to the box vectors (Bravais vectors). The angles should be in degrees. Parameters ---------- a_length : scalar or np.ndarray length of Bravais unit vector **a** b_length : scalar or np.ndarray length of Bravais unit vector **b** c_length : scalar or np.ndarray length of Bravais unit vector **c** alpha : scalar or np.ndarray angle between vectors **b** and **c**, in degrees. beta : scalar or np.ndarray angle between vectors **c** and **a**, in degrees. gamma : scalar or np.ndarray angle between vectors **a** and **b**, in degrees. Returns ------- a : np.ndarray If the inputs are scalar, the vectors will one dimesninoal (length 3). If the inputs are one dimension, shape=(n_frames, ), then the output will be (n_frames, 3) b : np.ndarray If the inputs are scalar, the vectors will one dimesninoal (length 3). If the inputs are one dimension, shape=(n_frames, ), then the output will be (n_frames, 3) c : np.ndarray If the inputs are scalar, the vectors will one dimesninoal (length 3). If the inputs are one dimension, shape=(n_frames, ), then the output will be (n_frames, 3) Examples -------- >>> import numpy as np >>> result = lengths_and_angles_to_box_vectors(1, 1, 1, 90.0, 90.0, 90.0) Notes ----- This code is adapted from gyroid, which is licensed under the BSD http://pythonhosted.org/gyroid/_modules/gyroid/unitcell.html """ if np.all(alpha < 2 * np.pi) and np.all(beta < 2 * np.pi) and np.all(gamma < 2 * np.pi): warnings.warn("All your angles were less than 2*pi. Did you accidentally give me radians?") alpha = alpha * np.pi / 180 beta = beta * np.pi / 180 gamma = gamma * np.pi / 180 a = np.array([a_length, np.zeros_like(a_length), np.zeros_like(a_length)]) b = np.array([b_length * np.cos(gamma), b_length * np.sin(gamma), np.zeros_like(b_length)]) cx = c_length * np.cos(beta) cy = c_length * (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) / np.sin(gamma) cz = np.sqrt(c_length * c_length - cx * cx - cy * cy) c = np.array([cx, cy, cz]) if not a.shape == b.shape == c.shape: raise TypeError("Shape is messed up.") # Make sure that all vector components that are _almost_ 0 are set exactly # to 0 tol = 1e-6 a[np.logical_and(a > -tol, a < tol)] = 0.0 b[np.logical_and(b > -tol, b < tol)] = 0.0 c[np.logical_and(c > -tol, c < tol)] = 0.0 return a.T, b.T, c.T
[docs] def box_vectors_to_lengths_and_angles(a, b, c): """Convert box vectors into the lengths and angles defining the box. Parameters ---------- a : np.ndarray the vector defining the first edge of the periodic box (length 3), or an array of this vector in multiple frames, where a[i,:] gives the length 3 array of vector a in each frame of a simulation b : np.ndarray the vector defining the second edge of the periodic box (length 3), or an array of this vector in multiple frames, where b[i,:] gives the length 3 array of vector a in each frame of a simulation c : np.ndarray the vector defining the third edge of the periodic box (length 3), or an array of this vector in multiple frames, where c[i,:] gives the length 3 array of vector a in each frame of a simulation Examples -------- >>> a = np.array([2,0,0], dtype=float) >>> b = np.array([0,1,0], dtype=float) >>> c = np.array([0,1,1], dtype=float) >>> l1, l2, l3, alpha, beta, gamma = box_vectors_to_lengths_and_angles(a, b, c) >>> (l1 == 2.0) and (l2 == 1.0) and (l3 == np.sqrt(2)) True >>> np.abs(alpha - 45) < 1e-6 True >>> np.abs(beta - 90.0) < 1e-6 True >>> np.abs(gamma - 90.0) < 1e-6 True Returns ------- a_length : scalar or np.ndarray length of Bravais unit vector **a** b_length : scalar or np.ndarray length of Bravais unit vector **b** c_length : scalar or np.ndarray length of Bravais unit vector **c** alpha : scalar or np.ndarray angle between vectors **b** and **c**, in degrees. beta : scalar or np.ndarray angle between vectors **c** and **a**, in degrees. gamma : scalar or np.ndarray angle between vectors **a** and **b**, in degrees. """ if not a.shape == b.shape == c.shape: raise TypeError("Shape is messed up.") if not a.shape[-1] == 3: raise TypeError("The last dimension must be length 3") if a.ndim not in [1, 2]: raise ValueError( "vectors must be 1d or 2d (for a vectorized " "operation on multiple frames)", ) last_dim = a.ndim - 1 a_length = np.sqrt(np.sum(a * a, axis=last_dim)) b_length = np.sqrt(np.sum(b * b, axis=last_dim)) c_length = np.sqrt(np.sum(c * c, axis=last_dim)) # we allow 2d input, where the first dimension is the frame index # so we want to do the dot product only over the last dimension alpha = np.arccos(np.einsum("...i, ...i", b, c) / (b_length * c_length)) beta = np.arccos(np.einsum("...i, ...i", c, a) / (c_length * a_length)) gamma = np.arccos(np.einsum("...i, ...i", a, b) / (a_length * b_length)) # convert to degrees alpha = alpha * 180.0 / np.pi beta = beta * 180.0 / np.pi gamma = gamma * 180.0 / np.pi return a_length, b_length, c_length, alpha, beta, gamma
def lengths_and_angles_to_tilt_factors( a_length, b_length, c_length, alpha, beta, gamma, ): """ Parameters ---------- a_length : scalar or np.ndarray length of Bravais unit vector **a** b_length : scalar or np.ndarray length of Bravais unit vector **b** c_length : scalar or np.ndarray length of Bravais unit vector **c** alpha : scalar or np.ndarray angle between vectors **b** and **c**, in degrees. beta : scalar or np.ndarray angle between vectors **c** and **a**, in degrees. gamma : scalar or np.ndarray angle between vectors **a** and **b**, in degrees. Returns ------- lx : scalar Extent in x direction ly : scalar Extent in y direction lz : scalar Extent in z direction xy : scalar Unit vector **b** tilt with respect to **a** xz : scalar Unit vector of **c** tilt with respect to **a** yz : scalar Unit vector of **c** tilt with respect to **b** """ lx = a_length xy = b_length * np.cos(np.deg2rad(gamma)) xz = c_length * np.cos(np.deg2rad(beta)) ly = np.sqrt(b_length**2 - xy**2) yz = (b_length * c_length * np.cos(np.deg2rad(alpha)) - xy * xz) / ly lz = np.sqrt(c_length**2 - xz**2 - yz**2) return np.array([lx, ly, lz, xy, xz, yz])