check unitcell conversion

test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "from MDAnalysis.lib import mdamath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _angle(a, b):\n",
    "    \"\"\"Angle between two vectors *a* and *b* in degrees.\n",
    "\n",
    "    If one of the lengths is 0 then the angle is returned as 0\n",
    "    (instead of `nan`).\n",
    "    \"\"\"\n",
    "    # This function has different limits than angle?\n",
    "\n",
    "    angle = np.arccos(np.dot(a, b) / (norm(a) * norm(b)))\n",
    "    if np.isnan(angle):\n",
    "        return 0.0\n",
    "    return np.rad2deg(angle)\n",
    "\n",
    "\n",
    "def norm(v):\n",
    "    r\"\"\"Calculate the norm of a vector v.\n",
    "\n",
    "    .. math:: v = \\sqrt{\\mathbf{v}\\cdot\\mathbf{v}}\n",
    "\n",
    "    This version is faster then numpy.linalg.norm because it only works for a\n",
    "    single vector and therefore can skip a lot of the additional fuss\n",
    "    linalg.norm does.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    v : array_like\n",
    "        1D array of shape (N) for a vector of length N\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        norm of the vector\n",
    "\n",
    "    \"\"\"\n",
    "    return np.sqrt(np.dot(v, v))\n",
    "\n",
    "\n",
    "def triclinic_box(x, y, z):\n",
    "    \"\"\"Convert the three triclinic box vectors to [A,B,C,alpha,beta,gamma].\n",
    "\n",
    "    Angles are in degrees.\n",
    "\n",
    "    * alpha  = angle(y,z)\n",
    "    * beta   = angle(x,z)\n",
    "    * gamma  = angle(x,y)\n",
    "\n",
    "    Note\n",
    "    ----\n",
    "    Definition of angles: http://en.wikipedia.org/wiki/Lattice_constant\n",
    "    \"\"\"\n",
    "    A, B, C = [norm(v) for v in (x, y, z)]\n",
    "    alpha = _angle(y, z)\n",
    "    beta = _angle(x, z)\n",
    "    gamma = _angle(x, y)\n",
    "    return np.array([A, B, C, alpha, beta, gamma])\n",
    "\n",
    "\n",
    "def triclinic_vectors(dimensions):\n",
    "    \"\"\"Convert `[A,B,C,alpha,beta,gamma]` to a triclinic box representation.\n",
    "\n",
    "    Original `code by Tsjerk Wassenaar`_ posted on the Gromacs mailinglist.\n",
    "\n",
    "    If *dimensions* indicates a non-periodic system (i.e. all lengths\n",
    "    0) then null-vectors are returned.\n",
    "\n",
    "    .. _code by Tsjerk Wassenaar:\n",
    "       http://www.mail-archive.com/gmx-users@gromacs.org/msg28032.html\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    dimensions : [A, B, C, alpha, beta, gamma]\n",
    "        list of box lengths and angles (in degrees) such as\n",
    "        ``[A,B,C,alpha,beta,gamma]``\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    numpy.array\n",
    "        numpy 3x3 array B, with ``B[0]`` as first box vector,\n",
    "        ``B[1]``as second vector, ``B[2]`` as third box vector.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    The first vector is always pointing along the X-axis,\n",
    "    i.e., parallel to (1, 0, 0).\n",
    "\n",
    "\n",
    "    .. versionchanged:: 0.7.6\n",
    "       Null-vectors are returned for non-periodic (or missing) unit cell.\n",
    "\n",
    "    \"\"\"\n",
    "    B = np.zeros((3, 3))\n",
    "    x, y, z, a, b, c = dimensions[:6]\n",
    "\n",
    "    if np.all(dimensions[:3] == 0):\n",
    "        return B\n",
    "\n",
    "    B[0][0] = x\n",
    "    if a == 90. and b == 90. and c == 90.:\n",
    "        B[1][1] = y\n",
    "        B[2][2] = z\n",
    "    else:\n",
    "        a = np.deg2rad(a)\n",
    "        b = np.deg2rad(b)\n",
    "        c = np.deg2rad(c)\n",
    "        B[1][0] = y * np.cos(c)\n",
    "        B[1][1] = y * np.sin(c)\n",
    "        B[2][0] = z * np.cos(b)\n",
    "        B[2][1] = z * (np.cos(a) - np.cos(b) * np.cos(c)) / np.sin(c)\n",
    "        B[2][2] = np.sqrt(z * z - B[2][0] ** 2 - B[2][1] ** 2)\n",
    "    return B\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "for angle in np.linspace(1, 90, num=1000):\n",
    "    dimensions = [10, 15, 20, 90, 90, angle]\n",
    "    tri = triclinic_vectors(dimensions)\n",
    "    np.testing.assert_almost_equal(triclinic_box(*tri), dimensions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "for angles in [[90, 90, 90], [60, 60, 90], [60, 60, 60], [71.53, 109.47, 71.53]]:\n",
    "    dimensions = [10, 15, 20] + angles\n",
    "    tri = triclinic_vectors(dimensions)\n",
    "    np.testing.assert_almost_equal(triclinic_box(*tri), dimensions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}