{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating coordination numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will read in a snapshot from an MD simulation and then calculate the coordination number distribution.  \n",
    "This example assumes that you read the basic example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pybop.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": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile('conf.dump')"
   ]
  },
  {
   "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.  \n",
    "We will try with both of them. First we try with cutoff system."
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.get_neighbors(method='cutoff', cutoff=3.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets get all the atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.get_allatoms()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "lets try accessing the coordination number of an atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atoms[0].get_coordination()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can loop over all atoms and create a histogram of the results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "coord = [atom.get_coordination() for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets plot and see the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Cutoff method')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF4pJREFUeJzt3Xm0ZWV95vHvIwUqimGoCyKDhTQY0WVaUxqnGJRogzLFxraIQ4m4WNo4tRqR1nZIJAG0cUjiQISAiAwSFBRtZZWUxKhgASogIjQilJTUxQFR0kDhr/84u+Lh1lt1Tw3nnnPrfj9rnXX2fvfw/s6uqvPUHs7eqSokSZrqQaMuQJI0ngwISVKTASFJajIgJElNBoQkqcmAkCQ1GRDSOiR5bJKrktyV5A1JHprkC0nuTPLZEdRzWpL3baJ1vSfJpzfFurR5MiA06yT5yyTLkvwmyYokX07yrAGXvTnJn69Hd28DllbVNlX1EeAwYCdgh6p68QaUP7Akr0zyjWH2Ia2LAaFZJcmbgQ8Bf0vvi3p34KPAIUPq8tHAtVPGf1RVq4bUnzQ2DAjNGkn+APhr4OiqOr+qfltV91XVF6rqr7p5HnAIJsm+SZZ3w2fQC5QvdHsfb+vaD05ybZJfJVma5HFd+9eA5wD/0M1/FvAu4CXd+JGNGt+T5LNJPt0dlro6yd5Jjk2yMsmtSZ7f/5mSnNLtCf00yfuSbNHV8HHg6V1fv+rrZrskF3XrvyzJnn3re0aS73SHwL6T5Bl90/ZI8vVuuYuB+Rv7Z6LNmwGh2eTpwEOAz23IwlX1cuAW4KCqenhVnZhkb+As4E3ABPAlegGyVVU9F/hX4HXd/IfT23M5pxs/ZS1dHQScAWwHXAV8hd6/tV3oBdwn+uY9HVgF/CfgScDzgVdX1XXAa4BvdX1t27fM4cB7u/XfCBwHkGR74CLgI8AOwEnARUl26Jb7DHAFvWD4G2Dxem1AzTkGhGaTHYA7NvHhnZcAF1XVxVV1H/AB4KHAM9a92Dr9a1V9pavzs/SC5/hu/WcDC5Jsm2Qn4ADgTd3e0Ergg8CiadZ/flVd3q3/TOA/d+0vBG6oqjOqalVVnQX8EDgoye7AU4D/VVX3VNWlwBc24jNqDpg36gKk9fBzYH6SeZswJB4F/GT1SFX9Lsmt9P63v6Fu7xv+d3qhdn/fOMDDu763BFYkWT3/g4Bbp1n/z/qG7+7WBVM+S+cn9D7Lo4BfVtVvp0zbbZq+NIe5B6HZ5FvA/wMOXcc8vwW27ht/5JTpU29ffBu9E88ApPdNvRvw0w0vc2C3AvcA86tq2+71iKp6/Fpqnc4DPktnd3qfZQW9cxcPmzJNWisDQrNGVd1J7yTxPyY5NMnWSbZMckCSE7vZvgu8IMn2SR5J79xCv9uBx/SNnwu8MMl+SbYE3kLvS/ubw/00UFUrgK8C/zvJI5I8KMmeSf6sr9Zdk2w14Cq/BOzdXQY8L8lLgH2AL1bVT4BlwHuTbNVdFnzQJv5I2swYEJpVquok4M3AO4FJev8Lfx3w+W6WM4DvATfT+/I9Z8oq/g54Z3fF0lur6nrgZcDfA3fQ+9I8qKruHfJHWe0VwFbAD4BfAucBO3fTvkbvEtufJbljuhVV1c+BA+mF3M/p/YbjwKpavexfAn8C/AJ4N/CpTfcxtDmKDwySJLW4ByFJajIgJElNBoQkqcmAkCQ1zeofys2fP78WLFgw6jIkaVa54oor7qiqienmm9UBsWDBApYtWzbqMiRpVkky9Rf3TR5ikiQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNc3qX1JLc8n5E4tmtL8XTZ49o/1p/LgHIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCbv5ipNYybvouodVDVOhrYHkeTUJCuTXNPX9v4kP0zy/SSfS7Jt37Rjk9yY5Pok/2VYdUmSBjPMQ0ynAftPabsYeEJVPRH4EXAsQJJ9gEXA47tlPppkiyHWJkmaxtACoqouBX4xpe2rVbWqG/02sGs3fAhwdlXdU1U/Bm4Enjqs2iRJ0xvlSepXAV/uhncBbu2btrxrkySNyEgCIsk7gFXAmaubGrPVWpY9KsmyJMsmJyeHVaIkzXkzHhBJFgMHAi+tqtUhsBzYrW+2XYHbWstX1clVtbCqFk5MTAy3WEmaw2Y0IJLsDxwDHFxVd/dNuhBYlOTBSfYA9gIun8naJEkPNLTfQSQ5C9gXmJ9kOfBuelctPRi4OAnAt6vqNVV1bZJzgR/QO/R0dFXdP6zaJEnTG1pAVNXhjeZT1jH/ccBxw6pHkrR+vNWGJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoaWkAkOTXJyiTX9LVtn+TiJDd079t17UnykSQ3Jvl+kicPqy5J0mCGuQdxGrD/lLa3A0uqai9gSTcOcACwV/c6CvjYEOuSJA1gaAFRVZcCv5jSfAhwejd8OnBoX/unqufbwLZJdh5WbZKk6c30OYidqmoFQPe+Y9e+C3Br33zLu7Y1JDkqybIkyyYnJ4darCTNZeNykjqNtmrNWFUnV9XCqlo4MTEx5LIkae6a6YC4ffWho+59Zde+HNitb75dgdtmuDZJUp+ZDogLgcXd8GLggr72V3RXMz0NuHP1oShJ0mjMG9aKk5wF7AvMT7IceDdwPHBukiOBW4AXd7N/CXgBcCNwN3DEsOqSJA1maAFRVYevZdJ+jXkLOHpYtUiS1t+4nKSWJI0ZA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNQ7vdt7Qxzp9YNKP9vWjy7BntT5oN3IOQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNIwmIJP8jybVJrklyVpKHJNkjyWVJbkhyTpKtRlGbJKlnxgMiyS7AG4CFVfUEYAtgEXAC8MGq2gv4JXDkTNcmSfq9UR1imgc8NMk8YGtgBfBc4Lxu+unAoSOqTZLECAKiqn4KfAC4hV4w3AlcAfyqqlZ1sy0Hdmktn+SoJMuSLJucnJyJkiVpThooIJIcmGSThEmS7YBDgD2ARwEPAw5ozFqt5avq5KpaWFULJyYmNkVJkqSGQb/0FwE3JDkxyeM2ss8/B35cVZNVdR9wPvAMYNvukBPArsBtG9mPJGkjDBQQVfUy4EnA/wX+Ocm3ukM922xAn7cAT0uydZIA+wE/AC4BDuvmWQxcsAHrliRtIgMfNqqqXwP/ApwN7Az8BXBlktevT4dVdRm9k9FXAld3NZwMHAO8OcmNwA7AKeuzXknSpjXQ8yCSHAwcAewJnAE8tapWJtkauA74+/XptKreDbx7SvNNwFPXZz2SRsPndcwNgz4w6DB6v1G4tL+xqu5O8qpNX5YkadQGPcS0Ymo4JDkBoKqWbPKqJEkjN2hAPK/R1ro0VZK0mVjnIaYkrwX+O7Bnku/3TdoG+LdhFiZJGq3pzkF8Bvgy8HfA2/va76qqXwytKknSyE0XEFVVNyc5euqEJNsbEpK0+RpkD+JAevdKKiB90wp4zJDqkiSN2DoDoqoO7N73mJlyJEnjYtCb9T0zycO64ZclOSnJ7sMtTZI0SoNe5vox4O4kfwS8DfgJvV9US5I2U4MGxKqqKnq36f5wVX2Y3qWukqTN1KC32rgrybHAy4BnJ9kC2HJ4ZUmSRm3QPYiXAPcAR1bVz+g97e39Q6tKkjRyA+1BdKFwUt/4LcCnhlWUJGn0Br2K6UVJbkhyZ5JfJ7krya+HXZwkaXQGPQdxInBQVV03zGIkSeNj0HMQtxsOkjS3DLoHsSzJOcDn6Z2sBqCqzh9KVZKkkRs0IB4B3A08v6+tAANCkjZTg17FdMSwC5EkjZdBr2LaO8mSJNd0409M8s7hliZJGqVBT1L/E3AscB9AVX0fWDSsoiRJozdoQGxdVZdPaVu1qYuRJI2PQQPijiR70jsxTZLDgBUb2mmSbZOcl+SHSa5L8vQk2ye5uPtB3sVJttvQ9UuSNt6gAXE08AngD5P8FHgT8JqN6PfDwP+pqj8E/gi4jt4zr5dU1V7AEh74DGxJ0gxb51VMSd7cN/ol4BJ6ofJb4L/Sd3+mQSV5BPBs4JUAVXUvcG+SQ4B9u9lOB5YCx6zv+iVJm8Z0exDbdK+FwGuB7YBt6e097LOBfT4GmAT+OclVST7ZPa1up6paAdC979haOMlRSZYlWTY5ObmBJUiSprPOgKiq91bVe4H5wJOr6q1V9Rbgj4FdN7DPecCTgY9V1ZPo7Y0MfDipqk6uqoVVtXBiYmIDS5AkTWfQcxC7A/f2jd8LLNjAPpcDy6vqsm78PHqBcXuSnQG695UbuH5J0iYw6K02zgAuT/I5elcy/QW98wTrrap+luTWJI+tquuB/YAfdK/FwPHd+wUbsn5J0qYx6K02jkvyZeBPu6Yjquqqjej39cCZSbYCbgKOoLc3c26SI4FbgBdvxPolSRtp0D0IqupK4MpN0WlVfZfeie+p9tsU65ckbbxBz0FIkuYYA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTSMLiCRbJLkqyRe78T2SXJbkhiTnJNlqVLVJkka7B/FG4Lq+8ROAD1bVXsAvgSNHUpUkCRhRQCTZFXgh8MluPMBzgfO6WU4HDh1FbZKknnkj6vdDwNuAbbrxHYBfVdWqbnw5sEtrwSRHAUcB7L777kMuc+45f2LRjPX1osmzZ6wvSetvxvcgkhwIrKyqK/qbG7NWa/mqOrmqFlbVwomJiaHUKEkazR7EM4GDk7wAeAjwCHp7FNsmmdftRewK3DaC2iRJnRnfg6iqY6tq16paACwCvlZVLwUuAQ7rZlsMXDDTtUmSfm+cfgdxDPDmJDfSOydxyojrkaQ5bVQnqQGoqqXA0m74JuCpo6xHkvR747QHIUkaIwaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkppmPCCS7JbkkiTXJbk2yRu79u2TXJzkhu59u5muTZL0e6PYg1gFvKWqHgc8DTg6yT7A24ElVbUXsKQblySNyLyZ7rCqVgAruuG7klwH7AIcAuzbzXY6sBQ4ZqbrkzS7nD+xaMb6etHk2TPW1zgY6TmIJAuAJwGXATt14bE6RHZcyzJHJVmWZNnk5ORMlSpJc87IAiLJw4F/Ad5UVb8edLmqOrmqFlbVwomJieEVKElz3EgCIsmW9MLhzKo6v2u+PcnO3fSdgZWjqE2S1DOKq5gCnAJcV1Un9U26EFjcDS8GLpjp2iRJvzfjJ6mBZwIvB65O8t2u7X8CxwPnJjkSuAV48QhqkyR1RnEV0zeArGXyfjNZiyRp7fwltSSpyYCQJDUZEJKkJgNCktQ0iquYNMVM3ioA5t7tAiRtGPcgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqmrMPDPIhPZI2pc3xO8U9CElS09gFRJL9k1yf5MYkbx91PZI0V41VQCTZAvhH4ABgH+DwJPuMtipJmpvGKiCApwI3VtVNVXUvcDZwyIhrkqQ5KVU16hr+Q5LDgP2r6tXd+MuBP6mq1/XNcxRwVDf6WOD6GS90fMwH7hh1EWPGbbImt8ma5vo2eXRVTUw307hdxZRG2wMSrKpOBk6emXLGW5JlVbVw1HWME7fJmtwma3KbDGbcDjEtB3brG98VuG1EtUjSnDZuAfEdYK8keyTZClgEXDjimiRpThqrQ0xVtSrJ64CvAFsAp1bVtSMua5x5qG1NbpM1uU3W5DYZwFidpJYkjY9xO8QkSRoTBoQkqcmAmCWSnJpkZZJr+tq2T3Jxkhu69+1GWeNMW8s2eX+SHyb5fpLPJdl2lDXOtNY26Zv21iSVZP4oahuVtW2TJK/vbutzbZITR1XfODMgZo/TgP2ntL0dWFJVewFLuvG55DTW3CYXA0+oqicCPwKOnemiRuw01twmJNkNeB5wy0wXNAZOY8o2SfIcendpeGJVPR74wAjqGnsGxCxRVZcCv5jSfAhwejd8OnDojBY1Yq1tUlVfrapV3ei36f2WZs5Yy98TgA8Cb2PKD0/ngrVsk9cCx1fVPd08K2e8sFnAgJjddqqqFQDd+44jrmfcvAr48qiLGLUkBwM/rarvjbqWMbI38KdJLkvy9SRPGXVB42isfgchbSpJ3gGsAs4cdS2jlGRr4B3A80ddy5iZB2wHPA14CnBukseU1/0/gHsQs9vtSXYG6N7dTQaSLAYOBF7qP3j2BPYAvpfkZnqH3K5M8siRVjV6y4Hzq+dy4Hf0buCnPgbE7HYhsLgbXgxcMMJaxkKS/YFjgIOr6u5R1zNqVXV1Ve1YVQuqagG9L8YnV9XPRlzaqH0eeC5Akr2BrZjbd3dtMiBmiSRnAd8CHptkeZIjgeOB5yW5gd4VKsePssaZtpZt8g/ANsDFSb6b5OMjLXKGrWWbzGlr2SanAo/pLn09G1js3uaavNWGJKnJPQhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJpTkty8+m6mSb65Eet5ZZJH9Y1/Msk+m6LGYUiyoHWHV2ldDAhttpKs81YyVfWMjVj9K4H/CIiqenVV/WAj1jfWptuW2jwZEBobSV7RPcfhe0nO6NoenWRJ174kye7TtJ+W5KQklwAnJNkhyVeTXJXkE0D6+vtN975vkqVJzuueJXFmknTT3pXkO0muSXJyeg4DFgJndj/Ge2i3/MJumcOTXN0tc0J/f0mO6z7ft5Ps1NgG7+meX7A0yU1J3tC1P2APoHu2w3u64aVJPpjk0iTXJXlKkvPTe07I+/pWPy/J6d02O6+7TxNJ/ri7Yd0VSb7Sd/uWpUn+NsnXgTdu5B+vZqOq8uVr5C/g8cD1wPxufPvu/Qv0fuUKvbuzfn6a9tOALwJbdOMfAd7VDb+Q3u2uV/fxm+59X+BOevcpehC9X90+q7+ObvgM4KBueCmwsG/aUnqh8Sh6z1yYoHdDuK8Bh3bzVN/yJwLvbGyH9wDfBB5M795APwe2BBYA1/TN91bgPX19n9ANvxG4Ddi5W8dyYIdu+QKe2c13areOLbv+Jrr2lwCn9q33o6P+u+FrdC/3IDQungucV1V3AFTV6vv3Px34TDd8BvCsadoBPltV93fDzwY+3a3zIuCXa+n/8qpaXlW/A75L7wsV4DndLaGv7mp8/DSf4ynA0qqarN5zKc7sagC4l154AVzR18dUF1XVPd22WAmssafRcGH3fjVwbVWtqN6zDm4Cduum3VpV/9YNf5reNnss8AS6W5MA7+SBz9A4Z4C+tZnyuKLGRRjsYTZrm6e//bcDLtPvnr7h++kdjnkI8FF6ewq3dod0HjLNerKOafdV1epa7mft//7WqIXercv7/0M3tY7Vy/xuyvK/6+tn6naort5rq+rpa6ll6rbUHOIehMbFEuC/JdkBes/b7tq/CSzqhl8KfGOa9qku7aaT5AB6zwAY1Oov4TuSPBw4rG/aXfRuCjjVZcCfJZmfZAvgcODr69Hn2twO7NidU3kwvduZr6/dk6wOgsPpbbPrgYnV7Um2TDLdXpLmCPcgNBaq6tokxwFfT3I/cBW9K4XeAJya5K+ASeCIbpG1tU/1XuCsJFfS+6Ie+JnMVfWrJP9E77DNzcB3+iafBnw8yb/TO9y1epkVSY4FLqH3v/MvVdVG34a9qu5L8tf0AujHwA83YDXXAYu7k/U3AB+rqnu7k+4fSfIH9L4TPgRcu7E1a/bzbq6SpCYPMUmSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpKb/D+zX/2keXaGpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": [
    "Now finally we can compare with the voronoi method"
   ]
  },
  {
   "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.  \n",
    "Voronoi polyhedra is calculated using the Voro++ code. However, you dont need to install this specifically as it is linked to pybop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.get_neighbors(method='voronoi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, lets get all atoms and find their coordination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.get_allatoms()\n",
    "coord = [atom.get_coordination() for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And visualise the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Voronoi method')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGRtJREFUeJzt3Xm8X3V95/HXm1VBFDAXhtWgRdweLkxErRsVnaJSoQxOYVyC0KEqdUeF6kPQliraAvqYkRmUlIAIIkahVkd5BJC6AAZEtohQQIgEchlkUVow+Jk/zon8vJwkP5L8ltz7ej4e5/E753u2z7k3+b3v2VNVSJI01QajLkCSNJ4MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQlpLSV6e5PoRrXt2kkqy0TpaXiX5o3WxLK3/DAiNpSTfSfKJjvZ9k9yxrr4Q14Wq+teq2m0Y60pyS5JXD2NdkgGhcXUq8JYkmdL+FuCMqlr+WBY2ToEirS8MCI2rbwBbAy9f0ZBkK2Af4LR2+ElJTksymeQXST6aZIN23MFJfpDkhCR3A8ck2aCd5hdJlrXzPqmdfsWhmrlJbk1yV5KP9Kx70yQnJrm97U5Msmk7bs8kS1a2Ie1y35nkhiT3J/nbJE9L8qMk9yU5O8kmPdPvk+TKJPck+WGS57btpwM7A/+c5NdJPtSzmjc91rrb8R9MsrQdd8ia/KI0jVWVnd1YdsAXgC/2DP8VcGXP8GnAucAWwGzg58Ch7biDgeXAu4CNgMcDhwA3Ak8FngAsAE5vp58NVLvOxwPPAx4EntmO/wRwCbANMAH8EPjbdtyewJJVbEcB5wFPBJ7dLndhW8eTgOuAue20uwPLgBcBGwJzgVuATdvxtwCv7ln22tS9N3An8Bxgc+DL7bL+aNS/e7vx6EZegJ3dyjrgZcC9wOPb4R8A72v7N2y/CJ/VM/1fARe1/QcDt05Z3kLgnT3DuwG/bQNkxRftjj3jLwMObPv/DXhdz7g/BW5p+/sJiJf2DF8OfLhn+B+BE9v+k1Z8gfeMvx54Zdu/soBYk7rnAZ/qGfd0A8Kut/MQk8ZWVX0fmAT2TfJU4IU0f+UCzAI2AX7RM8svgB16hm+bssjtO6bfCNi2p+2Onv4HaPY0Vjbv9v1uC81f6iv8e8fwivU8BfhAe3jpniT3ADv1sa41qXt7/vBn1DudZEBo7J0GvJXm5PR3q2rFF+tdNH/9P6Vn2p2BX/YMT31U8e0d0y/nD7+sV6Zr3tv7mO+xug04tqq27Ok2q6oz2/GP9fHLq6p7KU349I6Tfs+A0Lg7DXg18D+A+Ssaq+ph4Gzg2CRbJHkK8H7gS6tY1pnA+5LskuQJwN8DX6n+rog6E/hokokks4CPrWZda+oLwNuTvCiNzZO8PskW7fg7ac5d9GtVdZ8NHJzkWUk2A45eVxuh6cGA0FirqltoTqxuTnOit9e7gN8ANwHfpzn8NG8Vi5sHnA5cDNwM/Ee7jH78HbAIuAq4GriibVunqmoRTRj+T+BXNCfVD+6Z5JM0X/j3JDmij0WutO6q+jZwInBBu54L1s1WaLpIlS8MkiQ9mnsQkqROBoQkqZMBIUnqZEBIkjqt1w8wmzVrVs2ePXvUZUjSeuXyyy+/q6omVjfdeh0Qs2fPZtGiRaMuQ5LWK0n6umveQ0ySpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjoZEJKkTuv1ndRa/yyYOHCgy99/8qyBLl+aSdyDkCR1GlhAJJmXZFmSazrGHZGk2nfk0r5793NJbkxyVZLdB1WXJKk/g9yDOBXYe2pjkp2A1wC39jS/Fti17Q4DThpgXZKkPgwsIKrqYuDujlEnAB8Cel+GvS9wWjUuAbZMst2gapMkrd5Qz0EkeQPwy6r66ZRROwC39Qwvadu6lnFYkkVJFk1OTg6oUknS0AIiyWbAR4CPdY3uaKuONqrq5KqaU1VzJiZW+74LSdIaGuZlrk8DdgF+mgRgR+CKJHvQ7DHs1DPtjsDtQ6xNkjTF0PYgqurqqtqmqmZX1WyaUNi9qu4AzgPe2l7N9GLg3qpaOqzaJEmPNsjLXM8EfgTslmRJkkNXMfm3gJuAG4EvAO8cVF2SpP4M7BBTVR20mvGze/oLOHxQtUiSHjvvpJYkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1GlgAZFkXpJlSa7paftMkp8luSrJ15Ns2TPuqCQ3Jrk+yZ8Oqi5JUn8GuQdxKrD3lLbzgedU1XOBnwNHASR5FnAg8Ox2ns8n2XCAtUmSVmNgAVFVFwN3T2n7blUtbwcvAXZs+/cFzqqqB6vqZuBGYI9B1SZJWr1RnoM4BPh2278DcFvPuCVt26MkOSzJoiSLJicnB1yiJM1cIwmIJB8BlgNnrGjqmKy65q2qk6tqTlXNmZiYGFSJkjTjbTTsFSaZC+wD7FVVK0JgCbBTz2Q7ArcPuzZJ0iOGGhBJ9gY+DLyyqh7oGXUe8OUkxwPbA7sClw2zNk1/CyYOHOjy9588a6DLl4ZtYAGR5ExgT2BWkiXA0TRXLW0KnJ8E4JKqentVXZvkbOA6mkNPh1fVw4OqTZK0egMLiKo6qKP5lFVMfyxw7KDqkSQ9Nt5JLUnqZEBIkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSp08ACIsm8JMuSXNPTtnWS85Pc0H5u1bYnyeeS3JjkqiS7D6ouSVJ/BrkHcSqw95S2I4GFVbUrsLAdBngtsGvbHQacNMC6JEl9GFhAVNXFwN1TmvcF5rf984H9etpPq8YlwJZJthtUbZKk1Rv2OYhtq2opQPu5Tdu+A3Bbz3RL2rZHSXJYkkVJFk1OTg60WEmaycblJHU62qprwqo6uarmVNWciYmJAZclSTPXsAPizhWHjtrPZW37EmCnnul2BG4fcm2SpB7DDojzgLlt/1zg3J72t7ZXM70YuHfFoShJ0mhsNKgFJzkT2BOYlWQJcDTwKeDsJIcCtwJvbCf/FvA64EbgAeBtg6pLktSfgQVEVR20klF7dUxbwOGDqkWS9NgNLCA0vhZMHDjQ5e8/edZAly9pOMblKiZJ0pgxICRJnQwISVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUqa+ASLJPEsNEkmaQfr/0DwRuSPLpJM8cZEGSpPHQV0BU1ZuBFwD/BvxTkh8lOSzJFgOtTpI0Mn0fNqqq+4CvAWcB2wF/DlyR5F0Dqk2SNEL9noN4Q5KvAxcAGwN7VNVrgecBRwywPknSiPT7RrkDgBOq6uLexqp6IMkh674sSdKo9XuIaenUcEhyHEBVLVznVUmSRq7fgHhNR9tr13SlSd6X5Nok1yQ5M8njkuyS5NIkNyT5SpJN1nT5kqS1t8qASPKOJFcDz0hyVU93M3DVmqwwyQ7Au4E5VfUcYEOay2iPozmMtSvwK+DQNVm+JGndWN05iC8D3wY+CRzZ035/Vd29lut9fJLfApsBS4FXAf+9HT8fOAY4aS3WIUlaC6s7xFRVdQtwOHB/T0eSrddkhVX1S+AfgFtpguFe4HLgnqpa3k62BNiha/72/otFSRZNTk6uSQmSpD6sLiC+3H5eDixqPy/vGX7MkmwF7AvsAmwPbE73+Yzqmr+qTq6qOVU1Z2JiYk1KkCT1YZWHmKpqn/Zzl3W4zlcDN1fVJECSBcAfA1sm2ajdi9gRuH0drlOS9Bj1e6PcS5Ns3va/OcnxSXZew3XeCrw4yWZJAuwFXAdcSHO/BcBc4Nw1XL4kaR3o9zLXk4AHkjwP+BDwC+D0NVlhVV0KnANcAVzd1nAy8GHg/UluBJ4MnLImy5ckrRv93km9vKoqyb7AZ6vqlCRz13SlVXU0cPSU5puAPdZ0mZKkdavfgLg/yVHAm4FXJNmQ5plMkqRpqt9DTH8BPAgcWlV30FyC+pmBVSVJGrm+9iDaUDi+Z/hW4LRBFSVJGr1+r2Lav31G0r1J7ktyf5L7Bl2cJGl0+j0H8Wngz6pq8SCLkSSNj37PQdxpOEjSzNLvHsSiJF8BvkFzshqAqlowkKokSSPXb0A8EXgA+C89bQUYEJI0TfV7FdPbBl2IJGm89HsV09OTLExyTTv83CQfHWxpkqRR6vck9ReAo4DfAlTVVTRvgZMkTVP9BsRmVXXZlLblnVNKkqaFfgPiriRPo32JT5IDaN4GJ0mapvq9iulwmkdyPyPJL4GbgTcNrCpJ0sitMiCSvL9n8Fs0L/XZAPgN8F/peT6TJGl6Wd0exBbt527AC2ne8hbgLcDFA6xLkjRiq3sn9ccBknwX2L2q7m+HjwG+OvDqJEkj0+85iJ2Bh3qGHwJmr/NqpGlqwcTgrwrff/Ksga9DM0u/AXE6cFmSr9NcyfTnwPyBVSVJGrm+LnOtqmOBtwG/Au4B3lZVn1zTlSbZMsk5SX6WZHGSlyTZOsn57Xsnzk+y1ZouX5K09vq9D4KquqKqPtt2P1nL9X4W+L9V9QzgecBi4EhgYVXtCixshyVJI9J3QKwrSZ4IvAI4BaCqHqqqe4B9eeSw1Xxgv2HXJkl6xNADAngqMAn8U5KfJPliks2BbatqKUD7uc0IapMktUYREBsBuwMnVdULaG666/twUpLDkixKsmhycnJQNUrSjDeKgFgCLKmqS9vhc2gC484k2wG0n8u6Zq6qk6tqTlXNmZiYGErBkjQTDT0gquoO4LYku7VNewHXAecBc9u2uTR3bUuSRqTf+yDWtXcBZyTZBLiJ5hLaDYCzkxwK3Aq8cUS1SZIYUUBU1ZXAnI5Rew27FklSt1Gcg5AkrQcMCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUaWUAk2TDJT5J8sx3eJcmlSW5I8pUkm4yqNknSaPcg3gMs7hk+DjihqnYFfgUcOpKqJEnAiAIiyY7A64EvtsMBXgWc004yH9hvFLVJkhqj2oM4EfgQ8Lt2+MnAPVW1vB1eAuzQNWOSw5IsSrJocnJy8JVK0gw19IBIsg+wrKou723umLS65q+qk6tqTlXNmZiYGEiNkiTYaATrfCnwhiSvAx4HPJFmj2LLJBu1exE7ArePoDZJUmvoexBVdVRV7VhVs4EDgQuq6k3AhcAB7WRzgXOHXZsk6RHjdB/Eh4H3J7mR5pzEKSOuR5JmtFEcYvq9qroIuKjtvwnYY5T1SJIeMU57EJKkMWJASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjqN9FEbkgZvwcSBA1/H/pNnDXwdGj73ICRJndyDGAH/opO0PnAPQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ2GHhBJdkpyYZLFSa5N8p62fesk5ye5of3cati1SZIeMYo9iOXAB6rqmcCLgcOTPAs4ElhYVbsCC9thSdKIDD0gqmppVV3R9t8PLAZ2APYF5reTzQf2G3ZtkqRHjPQcRJLZwAuAS4Ftq2opNCECbLOSeQ5LsijJosnJyWGVKkkzzsgCIskTgK8B762q+/qdr6pOrqo5VTVnYmJicAVK0gw3koBIsjFNOJxRVQva5juTbNeO3w5YNoraJEmNUVzFFOAUYHFVHd8z6jxgbts/Fzh32LVJkh4xisd9vxR4C3B1kivbtr8BPgWcneRQ4FbgjSOoTZLUGnpAVNX3gaxk9F7DrEWStHLeSS1J6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROo7hRTtIMsWDiwIGvY//Jswa+jpnKPQhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ28UU7StORNemtvxgaE/3gkadU8xCRJ6mRASJI6jV1AJNk7yfVJbkxy5KjrkaSZaqzOQSTZEPhfwGuAJcCPk5xXVdeNtjJJ6t90Occ5bnsQewA3VtVNVfUQcBaw74hrkqQZKVU16hp+L8kBwN5V9Zft8FuAF1XVX/dMcxhwWDu4G3D90AsdjVnAXaMuYgTc7pnF7R6Op1TVxOomGqtDTEA62v4gwarqZODk4ZQzPpIsqqo5o65j2NzumcXtHi/jdohpCbBTz/COwO0jqkWSZrRxC4gfA7sm2SXJJsCBwHkjrkmSZqSxOsRUVcuT/DXwHWBDYF5VXTvissbFjDus1nK7Zxa3e4yM1UlqSdL4GLdDTJKkMWFASJI6GRBjKMm8JMuSXNPT9pkkP0tyVZKvJ9lylDUOQtd294w7IkklmTWK2gZpZdud5F3tY2euTfLpUdU3KCv5d/78JJckuTLJoiR7jLLGQUiyU5ILkyxuf7fvadu3TnJ+khvaz61GXasBMZ5OBfae0nY+8Jyqei7wc+CoYRc1BKfy6O0myU40j1+5ddgFDcmpTNnuJH9C8xSB51bVs4F/GEFdg3Yqj/59fxr4eFU9H/hYOzzdLAc+UFXPBF4MHJ7kWcCRwMKq2hVY2A6PlAExhqrqYuDuKW3frarl7eAlNPeITCtd2906AfgQU26anC5Wst3vAD5VVQ+20ywbemEDtpLtLuCJbf+TmIb3QVXV0qq6ou2/H1gM7EDzB8H8drL5wH6jqfARBsT66RDg26MuYhiSvAH4ZVX9dNS1DNnTgZcnuTTJ95K8cNQFDcl7gc8kuY1mr2k67in/XpLZwAuAS4Ftq2opNCECbDO6yhoGxHomyUdodlHPGHUtg5ZkM+AjNIcaZpqNgK1oDkF8EDg7SdejaKabdwDvq6qdgPcBp4y4noFJ8gTga8B7q+q+UdfTxYBYjySZC+wDvKlmxg0sTwN2AX6a5Baaw2pXJPlPI61qOJYAC6pxGfA7mge6TXdzgQVt/1dpnvA87STZmCYczqiqFdt7Z5Lt2vHbASM/rGhArCeS7A18GHhDVT0w6nqGoaqurqptqmp2Vc2m+dLcvaruGHFpw/AN4FUASZ4ObMLMeMrp7cAr2/5XATeMsJaBaPcETwEWV9XxPaPOowlI2s9zh13bVN5JPYaSnAnsSfMX453A0TTHYjcF/l872SVV9faRFDggXdtdVaf0jL8FmFNV0+qLciW/79OBecDzgYeAI6rqglHVOAgr2e7rgc/SHGL7D+CdVXX5qGochCQvA/4VuJpmzxDgb2jOQ5wN7Exzxd4bq6rroo2hMSAkSZ08xCRJ6mRASJI6GRCSpE4GhCSpkwEhSepkQGhGSXLLiifCJvnhWizn4CTb9wx/sX3g2lhKMrvrKbnSqhgQmraSrPKVulX1x2ux+IOB3wdEVf1lVV23Fssba6v7WWp6MiA0NpK8tX3fxU+TnN62PSXJwrZ9YZKdV9N+apLjk1wIHJfkyUm+m+QnSf4PkJ71/br93DPJRUnOad+5ccaK5x4l+ViSHye5JsnJaRwAzAHOaN9b8Ph2/jntPAclubqd57je9SU5tt2+S5Js2/EzOKZ9T8JFSW5K8u62/Q/2ANr3YxzT9l+U5IQkF7fvGHhhkgXtewX+rmfxGyWZ3/7MzmmfdUWS/9w+EPDyJN/pedzDRUn+Psn3gPes5a9X66OqsrMbeQc8m+Yu2lnt8Nbt5z8Dc9v+Q4BvrKb9VOCbwIbt8OeAj7X9r6d5nPSKdfy6/dwTuJfmWU8bAD8CXtZbR9t/OvBnbf9FNHd10ztMs1dxKzBBczfwBcB+7TTVM/+ngY92/ByOAX5Ic9f8LJo75zcGZgPX9Ex3BHBMz7qPa/vfQ/O4iu3aZSwBntzOX8BL2+nmtcvYuF3fRNv+F8C8nuV+ftT/NuxG17kHoXHxKuCcah+jUY88YuAlwJfb/tOBl62mHeCrVfVw2/8K4EvtMv8F+NVK1n9ZVS2pqt8BV9J8oQL8SfvI7avbGp+9mu14IXBRVU1W8/6OM9oaoHlkxjfb/st71jHVv1TVg+3PYhnwqD2NDue1n1cD11bzzoEHgZuAndpxt1XVD9r+L9H8zHYDngOcn+RK4KP84btGvtLHujVNeVxR4yL090KglU3T2/6bPufp9WBP/8M0h2MeB3yeZk/htvaQzuNWs5xVPZL7t1W1opaHWfn/v0fVQvOI994/6KbWsWKe302Z/3c965n6c6i23mur6iUrqWXqz1IziHsQGhcLgf+W5MnQvJ+3bf8hcGDb/ybg+6tpn+ridjxJXkvzjoV+rfgSvqt9dv8BPePuB7bomOdS4JVJZiXZEDgI+N5jWOfK3Als055T2ZTmse+P1c5JVgTBQTQ/s+uBiRXtSTZOsrq9JM0Q7kFoLFTVtUmOBb6X5GHgJzRXCr0bmJfkg8Ak8LZ2lpW1T/Vx4MwkV9B8Uff9XuuquifJF2gO29wC/Lhn9KnA/07y7zSHu1bMszTJUcCFNH+df6uq1vqxzVX12ySfoAmgm4GfrcFiFgNz25P1NwAnVdVD7Un3zyV5Es13wonAtWtbs9Z/Ps1VktTJQ0ySpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnq9P8BHDsZxvseyS0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": [
    "Both 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 [conda env:pathsampling]",
   "language": "python",
   "name": "conda-env-pathsampling-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}