{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Zmatrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import chemcoord as cc\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "water = cc.Cartesian.read_xyz('water_dimer.xyz', start_index=1).get_zmat()\n",
    "small = cc.Cartesian.read_xyz('MIL53_small.xyz', start_index=1).get_zmat()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <caption>Zmat</caption>\n",
       "<thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atom</th>\n",
       "      <th>b</th>\n",
       "      <th>bond</th>\n",
       "      <th>a</th>\n",
       "      <th>angle</th>\n",
       "      <th>d</th>\n",
       "      <th>dihedral</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O</td>\n",
       "      <td>$\\vec{0}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>56.385853</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O</td>\n",
       "      <td>2</td>\n",
       "      <td>2.351206</td>\n",
       "      <td>1</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>3</td>\n",
       "      <td>-16.755013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>1</td>\n",
       "      <td>-180.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>5</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>2</td>\n",
       "      <td>163.244987</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  atom       b      bond    a       angle    d    dihedral\n",
       "1    O  origin  0.000000  e_z    0.000000  e_x    0.000000\n",
       "2    H       1  0.910922  e_z   56.385853  e_x   -0.000000\n",
       "3    H       1  0.910922    2  107.000024  e_x   -0.000000\n",
       "4    O       2  2.351206    1  132.466298    3  -16.755013\n",
       "5    H       4  0.910922    2  132.466298    1 -180.000000\n",
       "6    H       4  0.910922    5  107.000024    2  163.244987"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the normal zmatrix with the only difference being that the upper right triangle is filled with references to the origin and canonical unit vectors. Chemically speaking this fixes translational and rotational degrees of freedom and allows to preserve the information about the absolute position in space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we concentrate on operations you can do with zmatrices.\n",
    "\n",
    "Keep in mind, that there is an own tutorial dedicated to transformation from cartesian space to internal coordinates and back."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The slicing operations are the same as for ``pandas.DataFrames``. (http://pandas.pydata.org/pandas-docs/stable/indexing.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the ``'x'`` axis is of particular interest you can slice it out with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.000000\n",
       "2    0.910922\n",
       "3    0.910922\n",
       "4    2.351206\n",
       "5    0.910922\n",
       "6    0.910922\n",
       "Name: bond, dtype: float64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water['bond']\n",
    "# or explicit label based indexing\n",
    "water.loc[:, 'bond']\n",
    "# or explicit integer based indexing\n",
    "water.iloc[:, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is very easy to e.g. count the atoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H    4\n",
       "O    2\n",
       "Name: atom, dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water['atom'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Returned type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The indexing behaves like Indexing and Selecting data in\n",
    "[Pandas](<http://pandas.pydata.org/pandas-docs/stable/indexing.html>).\n",
    "You can slice with `Zmat.loc[key]`, `Zmat.iloc[keys]`, and `Zmat[key]`.\n",
    "The only question is about the return type.\n",
    "If the information in the columns is enough to draw a molecule,\n",
    "an instance of the own class (e.g. `Zmat`)\n",
    "is returned.\n",
    "If the information in the columns is not enough to draw a molecule, there \n",
    "are two cases to consider:\n",
    "\n",
    "* A `pandas.Series` instance is returned for one dimensional slices\n",
    "* A `pandas.DataFrame` instance is returned in all other cases:\n",
    "\n",
    "        molecule.loc[:, ['atom', 'b', 'bond', 'a', 'angle', 'd', 'dihedral']] returns a `Zmat`.\n",
    "\n",
    "        molecule.loc[:, ['atom', 'b']]`` returns a `pandas.DataFrame`.\n",
    "\n",
    "        molecule.loc[:, 'atom']`` returns a `pandas.Series`.\n",
    "\n",
    "If rows are omitted, there is never a `Zmat` instance returned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There exist four different methods to perform assignments:\n",
    "`safe_loc`, `unsafe_loc`, `safe_iloc`, and `unsafe_iloc`.\n",
    "\n",
    "As in pandas `safe_loc` and `unsafe_loc` are used for label based indexing and\n",
    "`safe_iloc` and `unsafe_iloc` for row based indexing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The safe methods assert that assignments do not lead to zmatrices, which can't be transformed back to cartesian coordinates. They also insert dummy atoms where necessary.\n",
    "\n",
    "The unsafe methods are wrappers around their `pandas.DataFrame` counterparts and are a lot faster than the safe assignments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Unless there is a good (performance) reason, it is recommended to use the safe assignment methods!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Symbolic evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to use symbolic expressions from sympy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()\n",
    "d = sympy.Symbol('d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "symb_water = water.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "symb_water.safe_loc[4, 'bond'] = d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <caption>Zmat</caption>\n",
       "<thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atom</th>\n",
       "      <th>b</th>\n",
       "      <th>bond</th>\n",
       "      <th>a</th>\n",
       "      <th>angle</th>\n",
       "      <th>d</th>\n",
       "      <th>dihedral</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O</td>\n",
       "      <td>$\\vec{0}$</td>\n",
       "      <td>0</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>56.385853</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O</td>\n",
       "      <td>2</td>\n",
       "      <td>$d$</td>\n",
       "      <td>1</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>3</td>\n",
       "      <td>-16.755013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>1</td>\n",
       "      <td>-180.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>5</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>2</td>\n",
       "      <td>163.244987</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  atom       b      bond    a       angle    d    dihedral\n",
       "1    O  origin         0  e_z    0.000000  e_x    0.000000\n",
       "2    H       1  0.910922  e_z   56.385853  e_x   -0.000000\n",
       "3    H       1  0.910922    2  107.000024  e_x   -0.000000\n",
       "4    O       2         d    1  132.466298    3  -16.755013\n",
       "5    H       4  0.910922    2  132.466298    1 -180.000000\n",
       "6    H       4  0.910922    5  107.000024    2  163.244987"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "symb_water"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <caption>Zmat</caption>\n",
       "<thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atom</th>\n",
       "      <th>b</th>\n",
       "      <th>bond</th>\n",
       "      <th>a</th>\n",
       "      <th>angle</th>\n",
       "      <th>d</th>\n",
       "      <th>dihedral</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O</td>\n",
       "      <td>$\\vec{0}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>56.385853</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O</td>\n",
       "      <td>2</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>1</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>3</td>\n",
       "      <td>-16.755013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>1</td>\n",
       "      <td>-180.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>5</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>2</td>\n",
       "      <td>163.244987</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  atom       b      bond    a       angle    d    dihedral\n",
       "1    O  origin  0.000000  e_z    0.000000  e_x    0.000000\n",
       "2    H       1  0.910922  e_z   56.385853  e_x   -0.000000\n",
       "3    H       1  0.910922    2  107.000024  e_x   -0.000000\n",
       "4    O       2  2.000000    1  132.466298    3  -16.755013\n",
       "5    H       4  0.910922    2  132.466298    1 -180.000000\n",
       "6    H       4  0.910922    5  107.000024    2  163.244987"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "symb_water.subs(d, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cc.xyz_functions.view([symb_water.subs(d, i).get_cartesian() for i in range(2, 5)])\n",
    "\n",
    "# Uncomment if viewer cannot open molden files\n",
    "# for i in range(2, 5):\n",
    "#     symb_water.subs(d, i).get_cartesian().view()\n",
    "#     time.sleep(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Binary operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mathematical Operations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general rule is that mathematical operations using the binary operators\n",
    "``(+ - * /)`` and the unary operators ``(+ - abs)``\n",
    "are only applied to the ``['bond', 'angle', 'dihedral']`` columns.\n",
    "\n",
    "**Addition/Subtraction/Multiplication/Division**:\n",
    "The most common case is to add another Zmat instance.\n",
    "In this case it is tested, if the used references are the same.\n",
    "Afterwards the addition in the ``['bond', 'angle', 'dihedral']`` columns\n",
    "is performed.\n",
    "If you add a scalar to a Zmat it is added elementwise onto the\n",
    "``['bond', 'angle', 'dihedral']`` columns.\n",
    "If you add a 3-dimensional vector, list, tuple... the first element of this\n",
    "vector is added elementwise to the ``'bond'`` column of the\n",
    "Zmat instance and so on.\n",
    "The third possibility is to add a matrix with\n",
    "``shape=(len(Zmat), 3)`` which is again added elementwise.\n",
    "The same rules are true for subtraction, division and multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "distortion = water.copy()\n",
    "distortion.unsafe_loc[:, ['bond', 'angle', 'dihedral']] = 0\n",
    "distortion.safe_loc[4, 'bond'] = d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <caption>Zmat</caption>\n",
       "<thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atom</th>\n",
       "      <th>b</th>\n",
       "      <th>bond</th>\n",
       "      <th>a</th>\n",
       "      <th>angle</th>\n",
       "      <th>d</th>\n",
       "      <th>dihedral</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O</td>\n",
       "      <td>$\\vec{0}$</td>\n",
       "      <td>0</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>56.385853</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O</td>\n",
       "      <td>2</td>\n",
       "      <td>$d + 2.3512055093207$</td>\n",
       "      <td>1</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>3</td>\n",
       "      <td>-16.755013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>1</td>\n",
       "      <td>-180.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>5</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>2</td>\n",
       "      <td>163.244987</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  atom       b                 bond    a       angle    d    dihedral\n",
       "1    O  origin                    0  e_z    0.000000  e_x    0.000000\n",
       "2    H       1             0.910922  e_z   56.385853  e_x    0.000000\n",
       "3    H       1             0.910922    2  107.000024  e_x    0.000000\n",
       "4    O       2  d + 2.3512055093207    1  132.466298    3  -16.755013\n",
       "5    H       4             0.910922    2  132.466298    1 -180.000000\n",
       "6    H       4             0.910922    5  107.000024    2  163.244987"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water + distortion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <caption>Zmat</caption>\n",
       "<thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atom</th>\n",
       "      <th>b</th>\n",
       "      <th>bond</th>\n",
       "      <th>a</th>\n",
       "      <th>angle</th>\n",
       "      <th>d</th>\n",
       "      <th>dihedral</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O</td>\n",
       "      <td>$\\vec{0}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>$\\vec{e_z}$</td>\n",
       "      <td>56.385853</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H</td>\n",
       "      <td>1</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>$\\vec{e_x}$</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>O</td>\n",
       "      <td>2</td>\n",
       "      <td>5.351206</td>\n",
       "      <td>1</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>3</td>\n",
       "      <td>-16.755013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>2</td>\n",
       "      <td>132.466298</td>\n",
       "      <td>1</td>\n",
       "      <td>-180.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H</td>\n",
       "      <td>4</td>\n",
       "      <td>0.910922</td>\n",
       "      <td>5</td>\n",
       "      <td>107.000024</td>\n",
       "      <td>2</td>\n",
       "      <td>163.244987</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  atom       b      bond    a       angle    d    dihedral\n",
       "1    O  origin  0.000000  e_z    0.000000  e_x    0.000000\n",
       "2    H       1  0.910922  e_z   56.385853  e_x    0.000000\n",
       "3    H       1  0.910922    2  107.000024  e_x    0.000000\n",
       "4    O       2  5.351206    1  132.466298    3  -16.755013\n",
       "5    H       4  0.910922    2  132.466298    1 -180.000000\n",
       "6    H       4  0.910922    5  107.000024    2  163.244987"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(water + distortion).subs(d, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(water + distortion).subs(d, 3).get_cartesian().view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Movements without small angular approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "cartesians = cc.xyz_functions.read_molden('MIL53_ht_lt_movement.molden')\n",
    "STEPS = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A typical approximation is the small angular approximation when calculating activation barriers for the transition from one allotrope to another.\n",
    "\n",
    "Let's import the following allotropes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m1, m2 = cartesians[0], cartesians[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest approach for obtaining structures for the movement from ``m1`` to ``m2``, is to linearize the movement in cartesian space.\n",
    "\n",
    "For geometric reasons this leads to shortened bond lengths, when bond angles change. The larger the change in angles, the larger this effect which leads to an overestimation of the activation energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cc.xyz_functions.view([m1 + (m2 - m1) * i / STEPS for i in range(STEPS)])\n",
    "\n",
    "# Uncomment if viewer cannot open molden files\n",
    "# for molecule in [m1 + (m2 - m1) * i / STEPS for i in range(STEPS)]:\n",
    "#     molecule.view()\n",
    "#     time.sleep(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another approach is to build two Zmatrices with the exact same references/construction table and linearise the movement in the space of internal coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zm1 = m1.get_zmat()\n",
    "c_table = zm1.loc[:, ['b', 'a', 'd']]\n",
    "zm2 = m2.get_zmat(c_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The term ``zm2 - zm1`` is not convertible to cartesian coordinates.\n",
    "For this reason we have to switch of the testing for the validity of zmatrices when using operators.\n",
    "\n",
    "The ``minimize_dihedrals`` method chooses the minimal absolute value representation of an angle equivalence class.\n",
    "So it will move e.g. by -30° instead of 270°."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with cc.TestOperators(False):\n",
    "    zmats = [zm1 + (zm2 - zm1).minimize_dihedrals() * i / STEPS for i in range(STEPS)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting structures are a lot more reasonable for the interconversion of allotropes.\n",
    "(Look for example at the ``C-H`` distance in methyl groups)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cc.xyz_functions.view([x.get_cartesian() for x in zmats])\n",
    "\n",
    "# Uncomment if viewer cannot open molden files\n",
    "# for molecule in [x.get_cartesian() for x in zmats]:\n",
    "#     molecule.view()\n",
    "#     time.sleep(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}