{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating coordination numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will read in a configuration from an MD simulation and then calculate the coordination number distribution.  \n",
    "This example assumes that you read the [basic example](https://pyscal.readthedocs.io/en/latest/examples.html#basic-examples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.core as pc\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in a file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is setting up a system. We can create atoms and simulation box using the ``pyscal.crystal_structures`` module. Let us start by importing the module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.crystal_structures as pcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms, box = pcs.make_crystal('bcc', lattice_constant= 4.00, repetitions=[6,6,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above function creates an bcc crystal of 6x6x6 unit cells with a lattice constant of 4.00 along with a simulation box that encloses the particles. We can then create a ``System`` and assign the atoms and box to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.atoms = atoms\n",
    "sys.box = box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by calculating the neighbors of each atom in the system. There are two ways to do this, using a ``cutoff`` method and using a ``voronoi`` polyhedra method. We will try with both of them. First we try with cutoff system - which has three sub options. We will check each of them in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cutoff method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cutoff method takes cutoff distance value and finds all atoms within the cutoff distance of the host atom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.find_neighbors(method='cutoff', cutoff=4.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets get all the atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "let us try accessing the coordination number of an atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atoms[0].coordination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we would expect for a bcc type lattice, we see that the atom has 14 neighbors (8 in the first shell and 6 in the second). Lets try a more interesting example by reading in a bcc system with thermal vibrations. Thermal vibrations lead to distortion in atomic positions, and hence there will be a distribution of coordination numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile('conf.dump')\n",
    "sys.find_neighbors(method='cutoff', cutoff=3.6)\n",
    "atoms = sys.atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can loop over all atoms and create a histogram of the results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "coord = [atom.coordination for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets plot and see the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Cutoff method')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXUUlEQVR4nO3de7hddX3n8fdHAiqKcslpRC4NKuigj442Wm+1KNZBBWEYqsFbsPjw6OBt1Cq0jmJHpoBWrfVWKjSIyEUGBUVH8wQiWhUMikJAJENFkgZy8IIofYDgd/5YK6ubk5zk5CR773Oy36/n2c9e67du370I+3PWZf9WqgpJkgAeNOwCJEkzh6EgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCtImJHl8kmuS3JXkLUkemuTLSe5M8oUh1LM4yQe20bpOSvK5bbEubT8MBc06SV6ZZHmS3yZZk+RrSZ47xWV/luSFW7C5dwGXV9UuVfUx4ChgHrBHVf35NMqfsiTHJPl2P7chTWQoaFZJ8nbgo8D/pvly3hf4JHB4nzb5h8CKCeM/rap1fdqeNFSGgmaNJI8E/gY4vqouqqrfVdV9VfXlqvrLdp4HnF5JclCSVe3w2TQh8uX2KONdbfvLkqxI8usky5L8p7b9MuD5wMfb+c8F3gu8oh0/diM1npTkC0k+155yujbJAUlOTLI2ya1JXtT7mZKc0R7xrE7ygSQ7tDV8GnhWu61f92xmtySXtuu/Mslje9b37CTfb09vfT/Js3um7Zfkm+1yS4C5W/vfRNsfQ0GzybOAhwBfnM7CVfUa4OfAYVX18Ko6LckBwLnA24Ax4Ks0obFTVb0A+Bbwpnb+o2mOUM5vx8+YZFOHAWcDuwE/BL5O8//aXjSh9o898y4G1gGPA54KvAh4fVXdALwB+G67rV17llkIvL9d/0rgZIAkuwOXAh8D9gA+DFyaZI92uc8DV9OEwf8CFm3J/tNoMBQ0m+wB3LGNT928Ari0qpZU1X3Ah4CHAs/e9GKb9K2q+npb5xdowuaUdv3nAfOT7JpkHvAS4G3tUc9a4CM0X/qb8sWquqpd/znAf27bXwrcVFVnV9W6qjoX+AlwWJJ9gacD/7Oq7qmqK4Avb8Vn1HZqzrALkLbAL4C5SeZsw2B4NHDL+pGq+n2SW2n+qp+u23uG/50myO7vGQd4eLvtHYE1SdbP/yDg1s2s/7ae4bvbdcGEz9K6heazPBr4VVX9bsK0fTazLY0YjxQ0m3wXuAc4YhPz/A7YuWf8UROmT+wW+N9oLh4DkObbeR9g9fTLnLJbaT7P3KratX09oqqeOEmtm/OAz9Lal+azrKG5FvGwCdOkBzAUNGtU1Z00F3o/keSIJDsn2THJi5Oc1s52DfCSJLsneRTNtYJetwOP6Rm/AHhpkoOT7Ai8g+aL+jv9/TRQVWuAbwB/l+QRSR6U5LFJ/rSn1r2T7DTFVX4VOKC9ZXdOklcABwJfqapbgOXA+5Ps1N7Ce9g2/kjaDhgKmlWq6u+AtwPvAcZp/tp+E/CldpazgR8BP6P5wj1/wir+FnhPe6fRO6vqRuDVwD8Ad9B8UR5WVff2+aOs91pgJ+B64FfAhcCe7bTLaG6HvS3JHZtbUVX9AjiUJth+QfMbi0Orav2yrwT+GPgl8D7gs9vuY2h7ER+yI0lazyMFSVLHUJAkdQwFSVLHUJAkdWb1j9fmzp1b8+fPH3YZkjSrXH311XdU1djGps3qUJg/fz7Lly8fdhmSNKskmfjL946njyRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJnVn9i2ZplFw0tnCg2zty/LyBbk8zg0cKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOvaRKmzHI3kntmVTD1rcjhSRnJlmb5Lqetg8m+UmSHyf5YpJde6admGRlkhuT/Jd+1SVJmlw/Tx8tBg6Z0LYEeFJVPRn4KXAiQJIDgYXAE9tlPplkhz7WJknaiL6FQlVdAfxyQts3qmpdO/o9YO92+HDgvKq6p6r+FVgJPKNftUmSNm6YF5r/AvhaO7wXcGvPtFVtmyRpgIYSCkn+GlgHnDONZY9LsjzJ8vHx8W1fnCSNsIGHQpJjgEOBV1VVtc2rgX16Ztu7bdtAVZ1eVQuqasHY2Fhfa5WkUTPQUEhyCPAu4GVVdXfPpEuAhUkenGQ/YH/gqkHWJknq4+8UkpwLHATMTbIKeB/N3UYPBpYkAfheVb2hqlYkuQC4nua00vFVdX+/apMkbVzfQqGqjt5I8xmbmP9k4OR+1SNJ2jy7uZAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdfoWCknOTLI2yXU9bbsnWZLkpvZ9t7Y9ST6WZGWSHyd5Wr/qkiRNrp9HCouBQya0nQAsrar9gaXtOMCLgf3b13HAp/pYlyRpEn0Lhaq6AvjlhObDgbPa4bOAI3raP1uN7wG7JtmzX7VJkjZu0NcU5lXVmnb4NmBeO7wXcGvPfKvatg0kOS7J8iTLx8fH+1epJI2goV1orqoCahrLnV5VC6pqwdjYWB8qk6TRNehQuH39aaH2fW3bvhrYp2e+vds2SdIADToULgEWtcOLgIt72l/b3oX0TODOntNMkqQBmdOvFSc5FzgImJtkFfA+4BTggiTHArcAL29n/yrwEmAlcDfwun7VJUmaXN9CoaqOnmTSwRuZt4Dj+1WLJGlq/EWzJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOn3rOlvaGheNLRzo9o4cP2+g25NmKo8UJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEmdoYRCkv+RZEWS65Kcm+QhSfZLcmWSlUnOT7LTMGqTpFE28FBIshfwFmBBVT0J2AFYCJwKfKSqHgf8Cjh20LVJ0qgb1umjOcBDk8wBdgbWAC8ALmynnwUcMaTaJGlkDTwUqmo18CHg5zRhcCdwNfDrqlrXzrYK2Gtjyyc5LsnyJMvHx8cHUbIkjYwphUKSw5JskwBJshtwOLAf8GjgYcAhU12+qk6vqgVVtWBsbGxblCRJak31i/4VwE1JTkvyhK3c5guBf62q8aq6D7gIeA6wa3s6CWBvYPVWbkeStIWmFApV9WrgqcD/AxYn+W57GmeXaWzz58Azk+ycJMDBwPXA5cBR7TyLgIunsW5J0laY8imhqvoNzYXg84A9gf8K/CDJm7dkg1V1ZbueHwDXtjWcDrwbeHuSlcAewBlbsl5J0tab0vMUkhwOHAM8Dvgs8IyqWptkZ5q/8v9hSzZaVe8D3jeh+WbgGVuyHknD4fMutl9TfcjOkTS/Ibiit7Gq7k7i7wkkaTsx1dNHt00MhCSnAlTV0m1elSRpKKYaCn+2kbYXb8tCJEnDt8nTR0neCPx34LFJftwzaRfgX/pZmCRp8DZ3TeHzwNeAvwVO6Gm/q6p+2beqJElDsblQqKr6WZLjJ05IsrvBIEnbl6kcKRxK0zdRAemZVsBj+lSXJGkINhkKVXVo+77fYMqRJA3TVDvEe06Sh7XDr07y4ST79rc0SdKgTfWW1E8Bdyd5CvAOmj6Qzu5bVZKkoZhqKKyrqqLp8vrjVfUJmttSJUnbkal2c3FXkhOBVwPPa5+tsGP/ypIkDcOWPE/hHuDYqrqN5nkHH+xbVZKkoZjSkUIbBB/uGf85TW+pkqTtyFTvPjoyyU1J7kzymyR3JflNv4uTJA3WVK8pnAYcVlU39LMYSdJwTfWawu0GgiRt/6Z6pLA8yfnAl2guOANQVRf1pSpJ0lBMNRQeAdwNvKinrQBDQZK2I1O9++h1/S5EkjR8U7376IAkS5Nc144/Ocl7+luaJGnQpnqh+Z+AE4H7AKrqx8DCfhUlSRqOqYbCzlV11YS2ddu6GEnScE01FO5I8liai8skOQpYM92NJtk1yYVJfpLkhiTPSrJ7kiXtj+SWJNltuuuXJE3PVEPheOAfgSckWQ28DXjDVmz374H/W1VPAJ4C3EDzDOilVbU/sJQHPhNakjQAm7z7KMnbe0a/ClxOEyS/A/4bPf0hTVWSRwLPA44BqKp7gXuTHA4c1M52FrAMePeWrl+SNH2bO1LYpX0tAN4I7AbsSnOU8LRpbnM/YBz45yQ/TPKZ9qlu86pq/Smp24B5G1s4yXFJlidZPj4+Ps0SJEkbs8lQqKr3V9X7abrKflpVvbOq3gH8ETDdx3HOoQmUT1XVU2mOOh5wqqh9oE9NUtPpVbWgqhaMjY1NswRJ0sZM9ZrCPODenvF7meQv+SlYBayqqivb8QtpQuL2JHsCtO9rp7l+SdI0TbWbi88CVyX5Yjt+BLB4OhusqtuS3Jrk8VV1I3AwcH37WgSc0r5fPJ31S5Kmb6rdXJyc5GvAn7RNr6uqH27Fdt8MnJNkJ+Bm4HU0Ry0XJDkWuAV4+VasX5I0DVM9UqCqfgD8YFtstKquobl4PdHB22L9kqTpmeo1BUnSCDAUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEmdoYVCkh2S/DDJV9rx/ZJcmWRlkvOT7DSs2iRpVA3zSOGtwA0946cCH6mqxwG/Ao4dSlWSNMKGEgpJ9gZeCnymHQ/wAuDCdpazgCOGUZskjbI5Q9ruR4F3Abu043sAv66qde34KmCvjS2Y5DjgOIB99923z2WOnovGFg5sW0eOnzewbUmamoEfKSQ5FFhbVVdPZ/mqOr2qFlTVgrGxsW1cnSSNtmEcKTwHeFmSlwAPAR4B/D2wa5I57dHC3sDqIdQmSSNt4EcKVXViVe1dVfOBhcBlVfUq4HLgqHa2RcDFg65NkkbdTPqdwruBtydZSXON4Ywh1yNJI2dYF5oBqKplwLJ2+GbgGcOsR5JG3Uw6UpAkDZmhIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqDDwUkuyT5PIk1ydZkeStbfvuSZYkual9323QtUnSqBvGkcI64B1VdSDwTOD4JAcCJwBLq2p/YGk7LkkaoDmD3mBVrQHWtMN3JbkB2As4HDione0sYBnw7kHXJ2l2uWhs4cC2deT4eQPb1rAM9ZpCkvnAU4ErgXltYADcBsybZJnjkixPsnx8fHwgdUrSqBhaKCR5OPB/gLdV1W96p1VVAbWx5arq9KpaUFULxsbGBlCpJI2OoYRCkh1pAuGcqrqobb49yZ7t9D2BtcOoTZJG2TDuPgpwBnBDVX24Z9IlwKJ2eBFw8aBrk6RRN/ALzcBzgNcA1ya5pm37K+AU4IIkxwK3AC8fQm2SNNKGcffRt4FMMvngQdYiSXogf9EsSeoYCpKkjqEgSeoYCpKkzjDuPtIEg/yZPozGT/UlTY9HCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeqM7EN2fLCNpG1pe/lO8UhBktSZcaGQ5JAkNyZZmeSEYdcjSaNkRoVCkh2ATwAvBg4Ejk5y4HCrkqTRMaNCAXgGsLKqbq6qe4HzgMOHXJMkjYxU1bBr6CQ5Cjikql7fjr8G+OOqelPPPMcBx7WjjwduHHihM8dc4I5hFzHDuE825D7Z0Kjvkz+sqrGNTZh1dx9V1enA6cOuYyZIsryqFgy7jpnEfbIh98mG3CeTm2mnj1YD+/SM7922SZIGYKaFwveB/ZPsl2QnYCFwyZBrkqSRMaNOH1XVuiRvAr4O7ACcWVUrhlzWTOZptA25TzbkPtmQ+2QSM+pCsyRpuGba6SNJ0hAZCpKkjqEwSyQ5M8naJNf1tO2eZEmSm9r33YZZ46BNsk8+mOQnSX6c5ItJdh1mjYO2sX3SM+0dSSrJ3GHUNiyT7ZMkb27/raxIctqw6ptpDIXZYzFwyIS2E4ClVbU/sLQdHyWL2XCfLAGeVFVPBn4KnDjoooZsMRvuE5LsA7wI+PmgC5oBFjNhnyR5Pk1vCU+pqicCHxpCXTOSoTBLVNUVwC8nNB8OnNUOnwUcMdCihmxj+6SqvlFV69rR79H81mVkTPLvBOAjwLuAkbuzZJJ98kbglKq6p51n7cALm6EMhdltXlWtaYdvA+YNs5gZ6C+Arw27iGFLcjiwuqp+NOxaZpADgD9JcmWSbyZ5+rALmilm1O8UNH1VVUlG7q/AyST5a2AdcM6waxmmJDsDf0Vz6kj/YQ6wO/BM4OnABUkeU96j75HCLHd7kj0B2ncPgYEkxwCHAq/yf3IeC+wH/CjJz2hOp/0gyaOGWtXwrQIuqsZVwO9pOskbeYbC7HYJsKgdXgRcPMRaZoQkh9CcO39ZVd097HqGraqurao/qKr5VTWf5svwaVV125BLG7YvAc8HSHIAsBOj3Wtqx1CYJZKcC3wXeHySVUmOBU4B/izJTcAL2/GRMck++TiwC7AkyTVJPj3UIgdskn0y0ibZJ2cCj2lvUz0PWORRZcNuLiRJHY8UJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0EjJcnP1vcSmuQ7W7GeY5I8umf8M0kO3BY19kOS+RvrOVWayFDQdivJJrtxqapnb8XqjwG6UKiq11fV9Vuxvhltc/tS2w9DQTNGkte2z0H4UZKz27b5SS5r25cm2Xcz7YuTfDrJlcBpSfZI8o22z/zPAOnZ3m/b94OSLEtyYdu//jlJ0k57b5LvJ7kuyelpHAUsAM5pfyD30Hb5Be0yRye5tl3m1N7tJTm5/XzfS7JBB4ZJTmr7/1+W5OYkb+n5vL3PjXhnkpPa4WVJPpJkeZIbkjw9yUVpnrPxgZ7Vz2k/2w3tZ925Xf6P2k7hrk7y9Z6uU5Yl+WiS5cBbt/I/r2aLqvLla+gv4Ik0zz+Y247v3r5/mebXptD0evqlzbQvBr4C7NCOfwx4bzv8Upquo9dv47ft+0HAnTT9Aj2I5tevz+2tox0+GzisHV4GLOiZtowmKB5N88yCMZpO1y4DjmjnqZ7lTwPes5H9cBLwHeDBNH3x/ALYEZgPXNcz3zuBk3q2fWo7/Fbg34A923WsAvZoly/gOe18Z7br2LHd3ljb/grgzJ71fnLY/zZ8DfblkYJmihcAX6iqOwCqan3/988CPt8Onw08dzPttOu5vx1+HvC5dp2XAr+aZPtXVdWqqvo9cA3NlyjA89vula9ta3ziZj7H04FlVTVezXMdzmlrALiXJrAAru7ZxkSXVtU97b5Yy9S6RL+kfb8WWFFVa6p5VsDNwD7ttFur6l/a4c/R7LPHA0+i7RYEeA8PfAbF+VPYtrYjnifU9uh301jmnp7h+2lOtTwE+CTNEcGt7emah2xFXfdV1fp+Ze5n8v//NqiFphvw3j/iJtaxfpnfT1j+9z3bmdinTdGcTltRVc+apJbp7EvNYh4paKa4DPjzJHtA8/zptv07wMJ2+FXAtzbTPtEVwCvbdb4Y2JLnWK//4r0jycOBo3qm3UXT8d5EVwF/mmRukh2Ao4FvbsE2J3M78AftNZIH03QNvqX2TbL+y/+VwLeBG4Gx9e1JdkyyuaMhbcc8UtCMUFUrkpwMfDPJ/cAPae7weTPwz0n+EhgHXtcuMln7RO8Hzk2ygiZIpvyM4qr6dZJ/Aq6jebLd93smLwY+neTfaU5lrV9mTZITgMtp/gq/tKq2ukvzqrovyd/QhM5q4CfTWM2NwPFJzgSuBz5VVfe2F84/luSRNN8JHwVWbG3Nmp3sJVWS1PH0kSSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSp8/8BkYf05WEX33oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nos, counts = np.unique(coord, return_counts=True)\n",
    "plt.bar(nos, counts, color=\"#AD1457\")\n",
    "plt.ylabel(\"density\")\n",
    "plt.xlabel(\"coordination number\")\n",
    "plt.title(\"Cutoff method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive cutoff methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``pyscal`` also has adaptive cutoff methods implemented. These methods remove the restriction on having the same cutoff. A distinct cutoff is selected for each atom during runtime. ``pyscal`` uses two distinct algorithms to do this - ``sann`` and ``adaptive``. Please check the [documentation](https://pyscal.readthedocs.io/en/latest/nearestneighbormethods.html) for a explanation of these algorithms.  For the purpose of this example, we will use the `adaptive` algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``adaptive algorithm``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.find_neighbors(method='cutoff', cutoff='adaptive', padding=1.5)\n",
    "atoms = sys.atoms\n",
    "coord = [atom.coordination for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Cutoff adaptive method')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAY+ElEQVR4nO3deZQlZZ3m8e/TFi4oikCKbFqooK1Ot9qlra22C+rBFYZBxRVplDOO+9pqOy591Ba03dqtUaFwA5QuFUWnRaRk2gUtcGETRUQoLCBRQVCHRX7zR7wZXpPMqiQr772Zld/POfdkrG/84mZVPBkR976RqkKSJIC/GHcBkqTFw1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBS0pSe6e5AdJrkry4iS3SvLFJFcm+ewCtL82yXMXotYZ2n5dko8Oo+1hSnJBkkctUFtDe3+1MAwFbbYkT0+yLsnVSTYk+UqSh8xx3Zt6wHk1cHJVbVNV7wP2B3YEtq+qJ8+j/KFI8vAk6wenVdXbqmpRHxCTrE7ylnHXofExFLRZkrwceA/wNrqD852ADwL7DGmTdwbOmjb+k6q6fkjbk5aXqvLla14v4HbA1cCTN7LMauAtA+MPB9a34U8ANwB/aO28uk1/Et2B/wpgLfCXbfrXgT8C/68tfzRwLXBdGz94hu0/APh2a2sD8H7g5gPzHw38GLiyzfsG8Nw2765tm78CLgc+BWw7sO4FwGuBs4HfAEcCtwRu3fbphlbX1cDOwJuAT7Z1vwK8cFqtPwT2a8P3AE4Efg2cCzxlI+/xWuAtwLfatr4IbN/q/S3wPWDlwPIztg0c0t7La6faGdjPVwI/au/TscAtB9p7HnBea+94YOe5vL++Fudr7AX4WrovYG/gemDFRpZZzSyh0MYvAB41ML4n8Lt2MNmK7nLReVMH8nYAfO7A8v2Bdpbt/w3wQGAFsBI4B3hpm7cDcBXdJaitgJe1/ZkKhbu1Om4BTACnAO+ZVvuZwG7AdsA3p/Z1+n5OrxV4NvDNgXn3pAuuW9CFykXAQa3u+9KF0j1n2ce17T26K11Qnw38BHhUW//jwJFt2Y22Pf33NbCf36ULtu3ae/g/27xHtvXv12r/N+CUuby/vhbny8tH2hzbA5fXwl66eSpwQlWdWFXXAe8EbgX83Xwaq6rTquo7VXV9VV0A/DvwsDb7ccBZVXVc29Z7gEsG1j2v1XFNVU0C7xpYd8r7q+qiqvo18FbgaXMs7XPAfZLcuY0/A1hTVdcATwAuqKojW93fB/4D2Ng9kyOr6mdVdSXdWcjPqupr7XfzWbqDP/NsG+B9VfXLtp9fBO4zUPcRVXV6q/21wIOSrGQT768WpxXjLkBL2q+AHZKsWMBg2Bn4xdRIVd2Q5CJgl/k0lmRPuoP5KmBrun/zpw1s66KBbVXb1tS6OwLvBR4KbEN3D+430zZx0cDwL1qbm1RVVyU5ATgAOJQuTJ7XZt8Z+NskVwyssoLucttsLh0Y/sMM47fZjLbhzw/mv+dP+7kzcPrUjKq6Osmv6H5fG31/tTh5pqDN8W3gGmDfjSzzO7qD8ZQ7Tps/vZveX9IduABIErrLMxfPs8YP0V3T3qOqbgu8Dkibt6G1PX1bU97W6vtvbd1nDqw7ZXD5O7X64cb7NZOjgacleRDdvYiT2/SLgG9U1bYDr9tU1fPn0OambKrtm9pt8vTf163pziAvZtPvrxYhQ0Hz1i5VvAH4QJJ9k2ydZKskj01yWFvsB8DjkmyX5I7AS6c1cylwl4HxzwCPT7JXkq2AV9AFz7fmWeY2dDdbr05yD2DwwHoCcK8k+yVZAbyYPw+tbehuuF6ZZBfgVTO0/4IkuybZDvgnupuwU/u1fZLbbaS2L9MdUP8ZOLaqbmjTvwTsmeRZ7f3cKsn9k/zlTdnxWWyq7em/j005GjgoyX2S3IIuSE9tl+o29f5qETIUtFmq6l+BlwOvBybp/hJ9IfD5tsgn6D5VcwHwVf500JzyL8Drk1yR5JVVdS7dX+T/RncD84nAE6vq2nmW+Erg6XQ3PD8yuP2qupzuWvrb6S6F7UF3s3jKm+luoF5Jd4BbM0P7n277dT7wM7pPAVFVP6Y7YJ7f9u1Gl5XaNfg1dDeEPz0w/SrgMXSXln5Jd+nmULobuZtlDm1/DLhnq/nzM7fyZ+19DfjfdPclNtDd7D6gzdvU+6tFKFU+ZEeajyQX0H2S5mvjrkVaKJ4pSJJ6hoIkqeflI0lSzzMFSVJvSX95bYcddqiVK1eOuwxJWlJOO+20y6tqYqZ5SzoUVq5cybp168ZdhiQtKUl+Mds8Lx9JknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknpL+hvN0pZszcQBQ2t7v8ljhta2ljbPFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPZ+8Jqk3rKe9+aS3pWNoZwpJjkhyWZIzB6Ztl+TEJD9tP2/fpifJ+5Kcl+RHSe43rLokSbMb5uWj1cDe06a9BjipqvYATmrjAI8F9mivQ4APDbEuSdIshhYKVXUK8Otpk/cBjmrDRwH7Dkz/eHW+A2ybZKdh1SZJmtmobzTvWFUb2vAlwI5teBfgooHl1rdpN5LkkCTrkqybnJwcXqWStAyN7dNHVVVAzWO9w6tqVVWtmpiYGEJlkrR8jToULp26LNR+XtamXwzsNrDcrm2aJGmERh0KxwMHtuEDgS8MTH92+xTSA4ErBy4zSZJGZGjfU0hyNPBwYIck64E3Am8HPpPkYOAXwFPa4l8GHgecB/weOGhYdUmSZje0UKiqp80ya68Zli3gBcOqRZI0N3ZzIUnqGQqSpJ6hIEnq2SGeNEfD6iwO7DBOi4dnCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeoZCpKknqEgSeqNJRSSvCzJWUnOTHJ0klsm2T3JqUnOS3JskpuPozZJWs5GHgpJdgFeDKyqqnsDNwMOAA4F3l1VdwN+Axw86tokabkb1+WjFcCtkqwAtgY2AI8EjmvzjwL2HVNtkrRsjTwUqupi4J3AhXRhcCVwGnBFVV3fFlsP7DLT+kkOSbIuybrJyclRlCxJy8Y4Lh/dHtgH2B3YGbg1sPdc16+qw6tqVVWtmpiYGFKVkrQ8jePy0aOAn1fVZFVdB6wBHgxs2y4nAewKXDyG2iRpWRtHKFwIPDDJ1kkC7AWcDZwM7N+WORD4whhqk6RlbRz3FE6lu6F8OnBGq+Fw4B+Blyc5D9ge+Nioa5Ok5W7FphdZeFX1RuCN0yafDzxgDOVIkhq/0SxJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqTeWLrOliSANRMHDKXd/SaPGUq7y4FnCpKknqEgSerNKRSSPDGJASJJW7i5HuifCvw0yWFJ7jHMgiRJ4zOnUKiqZwL3BX4GrE7y7SSHJNlmqNVJkkZqzpeEquq3wHHAMcBOwH8HTk/yoiHVJkkasbneU9gnyeeAtcBWwAOq6rHAXwOvGF55kqRRmuv3FPYD3l1VpwxOrKrfJzl44cuSJI3DXC8fXTI9EJIcClBVJy14VZKksZhrKDx6hmmPXchCJEnjt9HLR0meD/wv4K5JfjQwaxvgm8MsTJI0epu6p/Bp4CvAvwCvGZh+VVX9emhVSZLGYlOhUFV1QZIXTJ+RZDuDQZK2LHM5U3gCcBpQQAbmFXCX+Ww0ybbAR4F7t3b+ATgXOBZYCVwAPKWqfjOf9iVJ87PRG81V9YT2c/equkv7OfWaVyA07wX+T1Xdg+67DufQXZ46qar2AE7izy9XSZJGYK5fXntwklu34WcmeVeSO81ng0luB/w98DGAqrq2qq4A9gGOaosdBew7n/YlSfM314+kfgj4fZKpbzD/DPjEPLe5OzAJHJnk+0k+2gJnx6ra0Ja5BNhxppVbn0vrkqybnJycZwmSpJnMNRSur6qi+2v+/VX1AbqPpc7HCuB+wIeq6r7A75h2qahtq2ZauaoOr6pVVbVqYmJiniVIkmYy11C4KslrgWcCJ7RnK2w1z22uB9ZX1alt/Di6kLg0yU4A7edl82xfkjRPN+V5CtcAB1fVJcCuwDvms8G2/kVJ7t4m7QWcDRwPHNimHQh8YT7tS5Lmb04d4rUD+bsGxi8EPr4Z230R8KkkNwfOBw6iC6jPtA72fgE8ZTPalyTNw5xCIcl+wKHAHei+qxC6S/+3nc9Gq+oHwKoZZu01n/YkSQtjrl1nHwY8sarOGWYxkqTxmmsoXGogaLFZM3HA0Nreb/KYobUtLWZzDYV1SY4FPk93wxmAqlozlKokSWMx11C4LfB74DED0wowFCRpCzLXTx8dNOxCJEnjN9e+j/ZMclKSM9v4XyV5/XBLkySN2ly/vPYR4LXAdQBV9SNgeHf5JEljMddQ2Lqqvjtt2vULXYwkabzmGgqXJ7krrZO6JPsDGza+iiRpqZnrp49eABwO3CPJxcDPgWcMrSpJ0lhsNBSSvHxg9MvAyXRnF78D/gcD/SFJkpa+TZ0pTD0z4e7A/el6Lg3wLGD6PQZJ0hK30VCoqjcDJDkFuF9VXdXG3wScMPTqJEkjNdcbzTsC1w6MX8ssj8uUJC1dc73R/HHgu0k+18b3BVYPpSJJ0tjMtZuLtyb5CvDQNumgqvr+8MqSJI3DXM8UqKrTgdOHWIskaczmek9BkrQMGAqSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hIEnqjS0UktwsyfeTfKmN757k1CTnJTk2yc3HVZskLVfjPFN4CXDOwPihwLur6m7Ab4CDx1KVJC1jYwmFJLsCjwc+2sYDPBI4ri1yFN0zGyRJIzSuM4X3AK8Gbmjj2wNXVNX1bXw9sMtMKyY5JMm6JOsmJyeHX6kkLSMjD4UkTwAuq6rT5rN+VR1eVauqatXExMQCVydJy9ucH7KzgB4MPCnJ44BbArcF3gtsm2RFO1vYFbh4DLVJ0rI28jOFqnptVe1aVSuBA4CvV9UzgJOB/dtiBwJfGHVtkrTcLabvKfwj8PIk59HdY/jYmOuRpGVnHJePelW1Fljbhs8HHjDOeiRpuVtMZwqSpDEzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPUNBktQzFCRJvbE+jlNbljUTBwyl3f0mjxlKu5JuzDMFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVJv5KGQZLckJyc5O8lZSV7Spm+X5MQkP20/bz/q2iRpuRvHmcL1wCuq6p7AA4EXJLkn8BrgpKraAzipjUuSRmjkoVBVG6rq9DZ8FXAOsAuwD3BUW+woYN9R1yZJy91Y7ykkWQncFzgV2LGqNrRZlwA7zrLOIUnWJVk3OTk5kjolabkYWygkuQ3wH8BLq+q3g/OqqoCaab2qOryqVlXVqomJiRFUKknLx1hCIclWdIHwqapa0yZfmmSnNn8n4LJx1CZJy9k4Pn0U4GPAOVX1roFZxwMHtuEDgS+MujZJWu7G8ZCdBwPPAs5I8oM27XXA24HPJDkY+AXwlDHUJmkL5oOgNm3koVBV/wVkltl7jbIWSdKf8xvNkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6o2jl1RJWhaG1SsrDK9nVs8UJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9P320BfN5tJJuKs8UJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1PMjqSOyFDvGkrT8eKYgSeoZCpKk3qIKhSR7Jzk3yXlJXjPueiRpuVk09xSS3Az4APBoYD3wvSTHV9XZw9ie1/gl6cYW05nCA4Dzqur8qroWOAbYZ8w1SdKykqoadw0AJNkf2LuqntvGnwX8bVW9cNpyhwCHtNG7A+eOqMQdgMtHtK1xcP+Wvi19H92/hXPnqpqYacaiuXw0V1V1OHD4qLebZF1VrRr1dkfF/Vv6tvR9dP9GYzFdProY2G1gfNc2TZI0IospFL4H7JFk9yQ3Bw4Ajh9zTZK0rCyay0dVdX2SFwL/CdwMOKKqzhpzWYNGfslqxNy/pW9L30f3bwQWzY1mSdL4LabLR5KkMTMUJEk9Q2ETkrwsyVlJzkxydJJbjrumzZXkiCSXJTlzYNp2SU5M8tP28/bjrHFzzLJ/70jy4yQ/SvK5JNuOs8bNMdP+Dcx7RZJKssM4alsos+1jkhe13+NZSQ4bV32ba5Z/o/dJ8p0kP0iyLskDxlGbobARSXYBXgysqqp7090AH17/GKOzGth72rTXACdV1R7ASW18qVrNjffvRODeVfVXwE+A1466qAW0mhvvH0l2Ax4DXDjqgoZgNdP2Mckj6Ho5+OuquhfwzjHUtVBWc+Pf4WHAm6vqPsAb2vjIGQqbtgK4VZIVwNbAL8dcz2arqlOAX0+bvA9wVBs+Cth3pEUtoJn2r6q+WlXXt9Hv0H0PZkma5fcH8G7g1cCS//TILPv4fODtVXVNW+aykRe2QGbZvwJu24Zvx5iONYbCRlTVxXR/jVwIbACurKqvjreqodmxqja04UuAHcdZzJD9A/CVcRexkJLsA1xcVT8cdy1DtCfw0CSnJvlGkvuPu6AF9lLgHUkuojvujOVs1lDYiHZdfR9gd2Bn4NZJnjneqoavus8pL/m/NmeS5J+A64FPjbuWhZJka+B1dJcctmQrgO2ABwKvAj6TJOMtaUE9H3hZVe0GvAz42DiKMBQ27lHAz6tqsqquA9YAfzfmmobl0iQ7AbSfS/bUfDZJngM8AXhGbVlf0Lkr3R8uP0xyAd2lsdOT3HGsVS289cCa6nwXuIGuE7ktxYF0xxiAz9L1HD1yhsLGXQg8MMnW7S+SvYBzxlzTsBxP94+S9vMLY6xlwSXZm+56+5Oq6vfjrmchVdUZVXWHqlpZVSvpDp73q6pLxlzaQvs88AiAJHsCN2fL6jX1l8DD2vAjgZ+OowhDYSOq6lTgOOB04Ay692tRfBV9cyQ5Gvg2cPck65McDLwdeHSSn9KdIb19nDVujln27/3ANsCJ7SN/Hx5rkZthlv3bosyyj0cAd2kf4zwGOHCpnvHNsn/PA/41yQ+Bt/GnRwSMtrYl+p5KkobAMwVJUs9QkCT1DAVJUs9QkCT1DAVJUs9Q0LKS5IKpHkSTfGsz2nlOkp0Hxj+a5J4LUeMwJFk5U6+q0nSGgrZYrRPDWVXV5nw7/Tl0XZ9MtfXcqjp7M9pb1Db1XmrLYSho0Ujy7Pa8gx8m+USbtjLJ19v0k5LcaRPTVyf5cJJTgcOSbJ/kq63//Y8CGdje1e3nw5OsTXJc66v/U1N96iR5Q5LvtedpHJ7O/sAq4FPti3C3auuvaus8LckZbZ1DB7eX5K1t/76T5EadDiZ5U+trf22S85O8eGB/B/vef2WSN7XhtUne3frgPyfJ/ZOsSfdsjLcMNL+i7ds5bV+3buv/Tetg7rQk/znQ3cnaJO9Jsg54yWb+erVUVJUvX2N/Afeie87BDm18u/bzi3TfXIWud9PPb2L6auBLwM3a+PuAN7Thx9N19De1javbz4cDV9L1GfQXdN80fchgHW34E8AT2/BauudsMDhOd/ZwITBB14Hb14F92zI1sP5hwOtneB/eBHwLuAVdvz6/ArYCVgJnDiz3SuBNA9s+tA2/hK67hJ1aG+uB7dv6BTy4LXdEa2Ortr2JNv2pwBED7X5w3P82fI325ZmCFotHAp+tqssBqmqqr/kHAZ9uw58AHrKJ6bR2/tiG/x74ZGvzBOA3s2z/u1W1vqpuAH5AdxAFeES6rprPaDXeaxP7cX9gbXWdKE71xvr3bd61dIEFcNrANqY7oaquae/FZcytG/Pj288zgLOqakN1zx04H9itzbuoqr7Zhj9J957dHbg3rfsP4PX8+bMmjp3DtrUF8TqhtkS/m8c61wwM/5HuUsstgQ/SnRFc1C7XbM7jWK+rqql+Zf7I7P//blQLXXffg3/ETa9jap0bpq1/w8B2pvdpU3SX086qqgfNUst83kstYZ4paLH4OvDkJNtD98zoNv1b/OkRqM8A/u8mpk93CvD01uZjgZvy7OmpA+/lSW4D7D8w7yq6Dvam+y7wsCQ7JLkZ8DTgGzdhm7O5FLhDu0dyC7ouwG+qOyWZOvg/Hfgv4FxgYmp6kq2SbOpsSFswzxS0KFTVWUneCnwjyR+B79N9wudFwJFJXgVMAge1VWabPt2bgaOTnEUXJHN+fnFVXZHkI8CZdE+j+97A7NXAh5P8ge5S1tQ6G5K8BjiZ7q/wE6pqs7shr6rrkvwzXehcDPx4Hs2cC7wgyRHA2cCHquraduP8fUluR3dMeA9w1ubWrKXJXlIlST0vH0mSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSev8fCo+CFUaJCoAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nos, counts = np.unique(coord, return_counts=True)\n",
    "plt.bar(nos, counts, color=\"#AD1457\")\n",
    "plt.ylabel(\"density\")\n",
    "plt.xlabel(\"coordination number\")\n",
    "plt.title(\"Cutoff adaptive method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The adaptive method also gives similar results!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Voronoi method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voronoi method calculates the voronoi polyhedra of all atoms. Any atom that shares a voronoi face area with the host atom are considered neighbors. Voronoi polyhedra is calculated using the [Voro++](http://math.lbl.gov/voro++/) code. However, you dont need to install this specifically as it is linked to pyscal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.find_neighbors(method='voronoi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, let us get all atoms and find their coordination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms\n",
    "coord = [atom.coordination for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And visualise the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Voronoi method')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAY2klEQVR4nO3deZhldX3n8feHVUFkLRk2bUYRB31cmBZ1cEExBhUDwzAG10bJEJUQF1xAfQRNiIoZQZ+M5EEhDYggMa0QlyhPIxIXwAaRVaRFlm6BLqIISgI2fueP++uTa1FFF9V9763uer+e5z73nN/vLN9T1X0/dc6555xUFZIkAWww6gIkSbOHoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgK0hpK8oIkN4xo3fOSVJKN1tLyKsmT1saytG4yFDQrJfmXJB+ZpP2AJHesrQ/BtaGq/rWqdh/GupLcnOSlw1iX5iZDQbPV6cDrk2RC+xuAs6pq5SNZ2GwKEWk2MxQ0W30F2BZ4waqGJFsD+wNntPEtk5yRZDzJLUk+mGSD1ndoku8lOTHJvwHHJdmgTXNLkhVt3i3b9KsOwyxIcmuSu5J8oG/dmyY5Kckv2uukJJu2vn2SLJtqQ9py35bkxiT3JvmrJE9M8v0k9yQ5N8kmfdPvn+TKJHe3aZ7e2s8EHg/8c5LfJHlv32pe90jrbv3vSXJ763vzTH5RWs9UlS9fs/IFfBb4XN/4nwNX9o2fAZwHbAHMA34KHNb6DgVWAkcCGwGPBt4MLAX+K/AYYBFwZpt+HlBtnY8GngHcD/y31v8R4BLgccAY8H3gr1rfPsCyh9mOanU+FnhqW+7iVseWwHXAgjbts4AVwHOADYEFwM3Apq3/ZuClfctek7r3A+4EngZsDnyhLetJo/7d+xrda+QF+PI11Qt4PnA38Kg2/j3gnW14Q+ABYI++6f8cuKgNHwrcOmF5i4G39Y3vDvyuhcaqD9ed+/ovAw5pwz8DXtHX98fAzW14OqGwd9/45cD7+sb/L3BSGz551Yd2X/8NwIva8FShMJO6TwM+1tf3ZEPBl4ePNGtV1XeBu4ADkzwR2IveX7MA2wEbA7f0zXILsFPf+G0TFrnjJNNvBGzf13ZH3/B99PYoppp3x+luC72/yFf590nGV63nCcBR7dDR3UnuBnaZxrpmUveO/OHPqH86zVGGgma7M4A3Aq8HvllVqz5M76L3V/4T+qZ9PLC8b3ziLYB/Mcn0K/nDD+ipTDbvL6Yx3yN1G3B8VW3V99qsqs5u/Y/0tsYPV/ft9AKnv09znKGg2e4M4KXA/6H3jSQAqupB4Fzg+CRbJHkC8C7g8w+zrLOBdybZNcljgL8BvljT+ybT2cAHk4wl2Q740GrWNVOfBd6S5Dnp2TzJK5Ns0frvpHcuYroeru5zgUOT7JFkM+DYtbURWncZCprVqupmeidHNwfOn9B9JPBb4Cbgu/QOLZ32MIs7DTgTuBj4OfAfbRnT8dfAEuAq4Grgita2VlXVEnoB+HfAr+idGD+0b5KP0vuQvzvJu6exyCnrrqpvACcBF7b1XLh2tkLrslT5kB1JUo97CpKkjqEgSeoYCpKkjqEgSeqs0zcJ22677WrevHmjLkOS1imXX375XVU1NlnfOh0K8+bNY8mSJaMuQ5LWKUmmvHrdw0eSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpM46fUWz1j2Lxg4Z6PIPGj9noMuX1nfuKUiSOgMLhSSnJVmR5JpJ+o5KUu2ZsbRn0X46ydIkVyXZc1B1SZKmNsg9hYXAfhMbk+wCvAy4ta/55cBu7XU4cPIA65IkTWFgoVBVFwO/nKTrROC9QP/DoQ8AzqieS4CtkuwwqNokSZMb6jmFJAcAy6vqxxO6dgJu6xtf1tomW8bhSZYkWTI+Pj6gSiVpbhpaKCTZDHg/8KE1WU5VnVJV86tq/tjYpM+IkCTN0DC/kvpEYFfgx0kAdgauSLIXsBzYpW/anVubJGmIhranUFVXV9XjqmpeVc2jd4hoz6q6AzgfeGP7FtJzgV9X1e3Dqk2S1DPIr6SeDfwA2D3JsiSHPczkXwduApYCnwXeNqi6JElTG9jho6p6zWr65/UNF3DEoGqRJE2PVzRLkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpM7BQSHJakhVJrulr+0SSnyS5KsmXk2zV13dMkqVJbkjyx4OqS5I0tUHuKSwE9pvQdgHwtKp6OvBT4BiAJHsAhwBPbfN8JsmGA6xNkjSJgYVCVV0M/HJC27eqamUbvQTYuQ0fAJxTVfdX1c+BpcBeg6pNkjS5UZ5TeDPwjTa8E3BbX9+y1vYQSQ5PsiTJkvHx8QGXKElzy0hCIckHgJXAWY903qo6parmV9X8sbGxtV+cJM1hGw17hUkOBfYH9q2qas3LgV36Jtu5tUmShmiooZBkP+C9wIuq6r6+rvOBLyT5JLAjsBtw2TBr0/pv0dghA13+QePnDHT50jAMLBSSnA3sA2yXZBlwLL1vG20KXJAE4JKqektVXZvkXOA6eoeVjqiqBwdVmyRpcgMLhap6zSTNpz7M9McDxw+qHknS6nlFsySpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjoDC4UkpyVZkeSavrZtklyQ5Mb2vnVrT5JPJ1ma5Kokew6qLknS1Aa5p7AQ2G9C29HA4qraDVjcxgFeDuzWXocDJw+wLknSFAYWClV1MfDLCc0HAKe34dOBA/vaz6ieS4CtkuwwqNokSZMb9jmF7avq9jZ8B7B9G94JuK1vumWt7SGSHJ5kSZIl4+Pjg6tUkuagkZ1orqoCagbznVJV86tq/tjY2AAqk6S5a9ihcOeqw0LtfUVrXw7s0jfdzq1NkjREww6F84EFbXgBcF5f+xvbt5CeC/y67zCTJGlINhrUgpOcDewDbJdkGXAs8DHg3CSHAbcAr26Tfx14BbAUuA9406DqkiRNbWChUFWvmaJr30mmLeCIQdUiSZqegYWCZq9FY4cMdPkHjZ8z0OVLGhxvcyFJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqTOtEIhyauSGCCStJ6b7gf9nwI3JjkhyVMGWZAkaXSmFQpV9XrgWcDPgIVJfpDk8CRbDLQ6SdJQTfuQUFXdA3wJOAfYAfifwBVJjhxQbZKkIZvuOYUDknwZuAjYGNirql4OPAM4anDlSZKGabpPXjsIOLGqLu5vrKr72vOWJUnrgekePrpjYiAk+ThAVS1e61VJkkZiuqHwR5O0vXymK03yziTXJrkmydlJHpVk1ySXJlma5ItJNpnp8iVJM/OwoZDkrUmuBp6S5Kq+18+Bq2aywiQ7AX8JzK+qpwEbAocAH6d3iOpJwK8AD0tJ0pCt7pzCF4BvAB8Fju5rv7eqfrmG6310kt8BmwG3Ay8BXtv6TweOA05eg3VIkh6h1R0+qqq6GTgCuLfvRZJtZrLCqloO/C1wK70w+DVwOXB3Va1sky0Ddpps/nZ9xJIkS8bHx2dSgiRpCqsLhS+098uBJe398r7xRyzJ1sABwK7AjsDmwH7Tnb+qTqmq+VU1f2xsbCYlSJKm8LCHj6pq//a+61pc50uBn1fVOECSRcDewFZJNmp7CzsDy9fiOiVJ0zDdi9f2TrJ5G359kk8mefwM13kr8NwkmyUJsC9wHfBt4OA2zQLgvBkuX5I0Q9P9SurJwH1JVl3B/DPgzJmssKoupXe7jCuAq1sNpwDvA96VZCmwLXDqTJYvSZq56V7RvLKqKskBwN9V1alrciVzVR0LHDuh+SZgr5kuU5K05qYbCvcmOQZ4PfDC9myFjQdXliRpFB7J8xTuBw6rqjvonQj+xMCqkiSNxLT2FFoQfLJv/FbgjEEVJUkajel+++igJDcm+XWSe5Lcm+SeQRcnSRqu6Z5TOAF4VVVdP8hiJEmjNd1zCncaCJK0/pvunsKSJF8EvkLvhDMAVbVoIFVJkkZiuqHwWOA+4GV9bQUYCpK0Hpnut4/eNOhCJEmjN91vHz05yeIk17Txpyf54GBLkyQN23RPNH8WOAb4HUBVXUXvaWmSpPXIdENhs6q6bELbykmnlCSts6YbCncleSK9k8skOZjeU9MkSeuR6X776Ah6t7d+SpLlwM+B1w2sKknSSDxsKCR5V9/o1+k9CGcD4LfA/6LvfkiSpHXf6vYUtmjvuwPPpvc0tABvACaeY5AkreNW94zmDwMkuRjYs6rubePHAV8beHWSpKGa7jmF7YEH+sYfaG2SpmHR2OC/wX3Q+DkDX4fWf9MNhTOAy5J8uY0fCCwcSEWSpJGZ1ldSq+p44E3Ar9rrTVX10ZmuNMlWSb6U5CdJrk/yvCTbJLmgPbfhgiRbz3T5kqSZme51ClTVFVX1qfb60Rqu91PAv1TVU4BnANcDRwOLq2o3YHEblyQN0bRDYW1JsiXwQuBUgKp6oKruBg4ATm+TnU7vEJUkaYiGHgrArsA48A9JfpTkc0k2B7avqlVXSd+BJ7IlaehGEQobAXsCJ1fVs+hdCPcHh4qqqmi31JgoyeFJliRZMj4+PvBiJWkuGUUoLAOWVdWlbfxL9ELiziQ7ALT3FZPNXFWnVNX8qpo/NjY2lIIlaa4YeihU1R3AbUl2b037AtcB5wMLWtsCeldPS5KGaLrXKaxtRwJnJdkEuIne1103AM5NchhwC/DqEdUmSXPWSEKhqq4E5k/Ste+wa5Ek/adRnFOQJM1ShoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqTOyUEiyYZIfJflqG981yaVJlib5YpJNRlWbJM1Vo9xTeDtwfd/4x4ETq+pJwK+Aw0ZSlSTNYSMJhSQ7A68EPtfGA7wE+FKb5HTgwFHUJklz2aj2FE4C3gv8vo1vC9xdVSvb+DJgp8lmTHJ4kiVJloyPjw++UkmaQ4YeCkn2B1ZU1eUzmb+qTqmq+VU1f2xsbC1XJ0lz20YjWOfewJ8keQXwKOCxwKeArZJs1PYWdgaWj6A2SZrThr6nUFXHVNXOVTUPOAS4sKpeB3wbOLhNtgA4b9i1SdJcN5uuU3gf8K4kS+mdYzh1xPVI0pwzisNHnaq6CLioDd8E7DXKeiRprptNewqSpBEzFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJnZHe5kLS4C0aO2Tg6zho/JyBr0PD4Z6CJKnjnsII+JebpNnKPQVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUmfooZBklyTfTnJdkmuTvL21b5PkgiQ3tveth12bJM11o9hTWAkcVVV7AM8FjkiyB3A0sLiqdgMWt3FJ0hANPRSq6vaquqIN3wtcD+wEHACc3iY7HThw2LVJ0lw30nMKSeYBzwIuBbavqttb1x3A9lPMc3iSJUmWjI+PD6VOSZorRhYKSR4D/BPwjqq6p7+vqgqoyearqlOqan5VzR8bGxtCpZI0d4wkFJJsTC8QzqqqRa35ziQ7tP4dgBWjqE2S5rJRfPsowKnA9VX1yb6u84EFbXgBcN6wa5OkuW4Ut87eG3gDcHWSK1vb+4GPAecmOQy4BXj1CGqTpDlt6KFQVd8FMkX3vsOsRZL0h7yiWZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUGcXFa5LmiEVjhwx8HQeNnzPwdcwl7ilIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSp48VrktZLXjg3M3M2FPwHI0kP5eEjSVLHUJAkdWZdKCTZL8kNSZYmOXrU9UjSXDKrzikk2RD4f8AfAcuAHyY5v6quG21lkjR96/I5y9m2p7AXsLSqbqqqB4BzgANGXJMkzRmpqlHX0ElyMLBfVf1ZG38D8Jyq+ou+aQ4HDm+juwM3DL3Q0dgOuGvURYyA2z23uN3D8YSqGpusY1YdPpqOqjoFOGXUdQxbkiVVNX/UdQyb2z23uN2jN9sOHy0Hdukb37m1SZKGYLaFwg+B3ZLsmmQT4BDg/BHXJElzxqw6fFRVK5P8BfBNYEPgtKq6dsRlzRZz7pBZ43bPLW73iM2qE82SpNGabYePJEkjZChIkjqGwiyU5LQkK5Jc09f2iSQ/SXJVki8n2WqUNQ7CZNvd13dUkkqy3ShqG6SptjvJke13fm2SE0ZV36BM8e/8mUkuSXJlkiVJ9hpljYOQZJck305yXfvdvr21b5PkgiQ3tvetR1GfoTA7LQT2m9B2AfC0qno68FPgmGEXNQQLeeh2k2QX4GXArcMuaEgWMmG7k7yY3tX8z6iqpwJ/O4K6Bm0hD/19nwB8uKqeCXyoja9vVgJHVdUewHOBI5LsARwNLK6q3YDFbXzoDIVZqKouBn45oe1bVbWyjV5C7xqO9cpk292cCLwXWC+/FTHFdr8V+FhV3d+mWTH0wgZsiu0u4LFteEvgF0Mtagiq6vaquqIN3wtcD+xE74+A09tkpwMHjqI+Q2Hd9GbgG6MuYhiSHAAsr6ofj7qWIXsy8IIklyb5TpJnj7qgIXkH8Ikkt9HbO1of94g7SeYBzwIuBbavqttb1x3A9qOoyVBYxyT5AL3dz7NGXcugJdkMeD+9wwhzzUbANvQOL7wHODdJRlvSULwVeGdV7QK8Ezh1xPUMTJLHAP8EvKOq7unvq961AiPZMzYU1iFJDgX2B15Xc+MCkycCuwI/TnIzvUNmVyT5LyOtajiWAYuq5zLg9/Rumra+WwAsasP/SO/OyeudJBvTC4SzqmrV9t6ZZIfWvwMwkkOGhsI6Isl+9I6r/0lV3Tfqeoahqq6uqsdV1byqmkfvg3LPqrpjxKUNw1eAFwMkeTKwCXPj7qG/AF7Uhl8C3DjCWgai7fGdClxfVZ/s6zqfXijS3s8bdm3gFc2zUpKzgX3o/WV4J3AsvWOrmwL/1ia7pKreMpICB2Sy7a6qU/v6bwbmV9V69eE4xe/7TOA04JnAA8C7q+rCUdU4CFNs9w3Ap+gdPvsP4G1VdfmoahyEJM8H/hW4mt4eIPQOk14KnAs8HrgFeHVVTfbFi8HWZyhIklbx8JEkqWMoSJI6hoIkqWMoSJI6hoIkqWMoaE5JcvOqO60m+f4aLOfQJDv2jX+u3dRsVkoyb7K7z0oTGQpabyV52MfNVtX/WIPFHwp0oVBVf1ZV163B8ma11f0stf4wFDRrJHlje17Ej5Oc2drmJbmwtS9O8vjVtC9M8vdJLgVOSLJtkm+1+9Z/Dkjf+n7T3vdJclGSL7XnF5y16j5DST6U5IdJrklySnoOBuYDZ7X7/j+6zT+/zfOaJFe3eT7ev74kx7ftuyTJQ254luS49pyBi5LclOQv+7a3/7kD705yXBu+KMmJ7fkD1yd5dpJF7b78f923+I3atl3ftnWzNv9/bzfduzzJN/tutXBRkpOSLAHevoa/Xq0rqsqXr5G/gKfSe07Edm18m/b+z8CCNvxm4CuraV8IfBXYsI1/GvhQG34lvZuMrVrHb9r7PsCv6d1baQPgB8Dz++tow2cCr2rDF9G7upr+cXp7D7cCY/Suyr0QOLBNU33znwB8cJKfw3HA9+ldvb4dvSvYNwbmAdf0Tfdu4Li+dX+8Db+d3q0idmjLWAZs2+YvYO823WltGRu39Y219j8FTutb7mdG/W/D13Bf7ilotngJ8I/VbmFR/3l5//OAL7ThM4Hnr6adtpwH2/ALgc+3ZX4N+NUU67+sqpZV1e+BK+l9iAK8uN2++upW41NXsx3PBi6qqvHqPf/irFYD9G5X8dU2fHnfOib6WlXd334WK5jeLZTPb+9XA9dW75799wM3Abu0vtuq6ntt+PP0fma7A08DLkhyJfBB/vBZHV+cxrq1HvE4odZHv53BPPf3DT9I71DLo4DP0NsjuK0drnnUGtT1u6padV+ZB5n6/99DaqF3u/T+P+Im1rFqnt9PmP/3feuZeE+bonc47dqqet4UtczkZ6l1mHsKmi0uBP53km2h97za1v594JA2/Dp6NxJ7uPaJLgZe25b5cuCRPPd21QfvXe3e9wf39d0LbDHJPJcBL0qyXZINgdcA33kE65zKncDj2jmSTendQv2RenySVR/+rwW+S+8GdGOr2pNsnGR1e0Naj7mnoFmhqq5NcjzwnSQPAj+i9w2fI4F/SPIeYBx4U5tlqvaJPgycneRaekEy7ec8V9XdST4LXEPvSVg/7OteCPx9kn+ndyhr1Ty3Jzka+Da9v8K/VlVrfAvkqvpdko/QC53lwE9msJgb6D0P+DTgOuDkqnqgnTj/dJIt6X0mnARcu6Y1a93kXVIlSR0PH0mSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOv8f/DP8Sc21TtcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nos, counts = np.unique(coord, return_counts=True)\n",
    "plt.bar(nos, counts, color=\"#AD1457\")\n",
    "plt.ylabel(\"density\")\n",
    "plt.xlabel(\"coordination number\")\n",
    "plt.title(\"Voronoi method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All methods find the coordination number, and the results are comparable. Cutoff method is very sensitive to the choice of cutoff radius, but voronoi method can slightly overestimate the neighbors due to thermal vibrations. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (myenv)",
   "language": "python",
   "name": "myenv"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}