{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entropy parameters\n",
    "\n",
    "In this example, the entropy parameters are calculated and used for distinction of solid and liquid. For a description of entropy parameters, see [here](http://pyscal.com/en/latest/methods/entropy_parameters/entropy_parameters.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.core as pc\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have two test configurations for Al at 900 K, one is fcc structured and the other one is in liquid state. We calculate the entropy parameters for each of these configurations. First we start by reading in the fcc configuration. For entropy parameters, the values of the integration limit $r_m$ is chosen as 1.4, based on the [original publication](https://aip.scitation.org/doi/10.1063/1.4998408). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(\"../tests/conf.fcc.Al.dump\")\n",
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values of $r_m$ is in units of lattice constant, so we need to calculate the lattice constant first. Since is a cubic box, we can do this by,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat = (sys.box[0][1]-sys.box[0][0])/5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate the entropy parameter and its averaged version. Averaging can be done in two methods, as a simple average over neighbors or using a switching function. We will use a simple averaging over the neighbors. The `local` keyword allows to use a local density instead of the global one. However, this only works if the neighbors were calculated using a cutoff method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.calculate_entropy(1.4*lat, averaged=True, local=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculated values are stored for each atom. This can be accessed as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms\n",
    "solid_entropy = [atom.entropy for atom in atoms]\n",
    "solid_avg_entropy = [atom.avg_entropy for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can quickly repeat the calculation for the liquid structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(\"../tests/conf.lqd.Al.dump\")\n",
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "lat = (sys.box[0][1]-sys.box[0][0])/5\n",
    "sys.calculate_entropy(1.4*lat, local=True, averaged=True)\n",
    "atoms = sys.atoms\n",
    "liquid_entropy = [atom.entropy for atom in atoms]\n",
    "liquid_avg_entropy = [atom.avg_entropy for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$s_s^i$')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAENCAYAAAASUO4dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASZklEQVR4nO3df4zkd13H8eeL6xnKD0Ox23r9xQleENLoFTcVU2OAUlMLodSooYnljOgRsbGkGC2QCGhUUChoouABDSc/Ck2QtDYIlAoCBop7cD2uXvXAnNj27C2iUoSIbd/+Md+T7d7u7XdnZ3bms/t8JJOZ73e+M/PK3N5rP/uZ73y/qSokSe161KQDSJLWxiKXpMZZ5JLUOItckhpnkUtS4yxySWrcKev5Yqeffnpt3759PV9Skpq3b9++r1XVzHL3r2uRb9++nbm5ufV8SUlqXpJ/Odn9Tq1IUuMscklqnEUuSY2zyCWpcRa5JDXOIpekxlnkktQ4i1ySGreuXwhq3eHrrz9h3Y5rr51AEkn6LkfkktQ4i1ySGmeRS1LjLHJJapxFLkmNs8glqXEWuSQ1rneRJ9mS5ItJbu2Wn5jktiSHu+vTxhdTkrSc1YzIrwEOLVi+Dri9qnYAt3fLkqR11qvIk5wDPA94x4LVlwN7u9t7gReONpokqY++I/K3AL8JPLxg3ZlVdRSguz5jxNkkST2sWORJng8cq6p9w7xAkt1J5pLMzc/PD/MUkqST6DMivwh4QZIjwPuB5yR5D3B/km0A3fWxpR5cVXuqaraqZmdmZkYUW5J03IpFXlWvrKpzqmo78CLgb6rqF4BbgF3dZruAm8eWUpK0rLXsR/564JIkh4FLumVJ0jpb1fHIq+qTwCe72/8OXDz6SJKk1fCbnZLUOItckhpnkUtS4yxySWqcRS5JjbPIJalxFrkkNc4il6TGWeSS1DiLXJIaZ5FLUuMscklqnEUuSY2zyCWpcas6jK1OdPj6609Yt+PaayeQRNJm5YhckhrX5+TLj07y+SR3Jrkryeu69a9Ncm+S/d3lsvHHlSQt1mdq5X+A51TVN5NsBT6T5K+7+95cVW8cXzxJ0kpWLPKqKuCb3eLW7lLjDCVJ6q/XHHmSLUn2A8eA26rqju6uq5McSHJDktOWeezuJHNJ5ubn50cUW5J0XK8ir6qHqmoncA5wYZLzgbcCTwF2AkeBNy3z2D1VNVtVszMzMyOKLUk6blV7rVTVfwKfBC6tqvu7gn8YeDtw4RjySZJW0GevlZkkT+hunwo8F7g7ybYFm10BHBxPREnSyfTZa2UbsDfJFgbFf1NV3Zrk3Ul2Mvjg8wjw0vHFlCQtp89eKweAC5ZYf9VYEkmSVsVvdkpS4yxySWqcRS5JjbPIJalxFrkkNc4il6TGWeSS1DiLXJIaZ5FLUuMscklqnEUuSY2zyCWpcRa5JDXOIpekxlnkktS4PmcIenSSzye5M8ldSV7XrX9iktuSHO6ulzz5siRpvPqMyP8HeE5V/QiDEy1fmuSZwHXA7VW1A7i9W5YkrbMVi7wGvtktbu0uBVwO7O3W7wVeOJaEkqST6jVHnmRLkv3AMeC2qroDOLOqjgJ012cs89jdSeaSzM3Pz48qtySp06vIq+qhqtoJnANcmOT8vi9QVXuqaraqZmdmZobNKUlaxqr2Wqmq/wQ+CVwK3J9kG0B3fWzk6SRJK+qz18pMkid0t08FngvcDdwC7Oo22wXcPK6QkqTlndJjm23A3iRbGBT/TVV1a5LPAjcleQnwVeDnxphTkrSMFYu8qg4AFyyx/t+Bi8cRSpLUn9/slKTGWeSS1DiLXJIaZ5FLUuMscklqnEUuSY2zyCWpcRa5JDXOIpekxlnkktQ4i1ySGmeRS1LjLHJJapxFLkmNs8glqXEWuSQ1rs+p3s5N8okkh5LcleSabv1rk9ybZH93uWz8cSVJi/U51duDwCuq6gtJHg/sS3Jbd9+bq+qN44snSVpJn1O9HQWOdrcfSHIIOHvcwSRJ/axqjjzJdgbn77yjW3V1kgNJbkhy2jKP2Z1kLsnc/Pz8msJKkk7Uu8iTPA74IPDyqvoG8FbgKcBOBiP2Ny31uKraU1WzVTU7MzMzgsiSpIV6FXmSrQxK/L1V9ZcAVXV/VT1UVQ8DbwcuHF9MSdJy+uy1EuCdwKGqun7B+m0LNrsCODj6eJKklfTZa+Ui4CrgS0n2d+teBVyZZCdQwBHgpWNJKEk6qT57rXwGyBJ3fXj0cSRJq+U3OyWpcRa5JDXOIpekxlnkktQ4i1ySGmeRS1LjLHJJapxFLkmNs8glqXEWuSQ1ziKXpMZZ5JLUOItckhpnkUtS4yxySWpcnzMEnZvkE0kOJbkryTXd+icmuS3J4e56yZMvS5LGq8+I/EHgFVX1NOCZwK8leTpwHXB7Ve0Abu+WJUnrbMUir6qjVfWF7vYDwCHgbOByYG+32V7gheMKKUla3qrmyJNsBy4A7gDOrKqjMCh74IxlHrM7yVySufn5+bWllSSdoHeRJ3kc8EHg5VX1jb6Pq6o9VTVbVbMzMzPDZJQknUSvIk+ylUGJv7eq/rJbfX+Sbd3924Bj44koSTqZPnutBHgncKiqrl9w1y3Aru72LuDm0ceTJK3klB7bXARcBXwpyf5u3auA1wM3JXkJ8FXg58YTcXo8eN99J6w75ayzJpBEkr5rxSKvqs8AWebui0cbR5K0Wn6zU5IaZ5FLUuMscklqnEUuSY2zyCWpcRa5JDXOIpekxlnkktQ4i1ySGmeRS1LjLHJJapxFLkmNs8glqXEWuSQ1rs/xyDetb+/bN+kIkrSiPmcIuiHJsSQHF6x7bZJ7k+zvLpeNN6YkaTl9plbeBVy6xPo3V9XO7vLh0caSJPW1YpFX1aeAr69DFknSENbyYefVSQ50Uy+njSyRJGlVhi3ytwJPAXYCR4E3Lbdhkt1J5pLMzc/PD/lykqTlDFXkVXV/VT1UVQ8DbwcuPMm2e6pqtqpmZ2Zmhs0pSVrGUEWeZNuCxSuAg8ttK0karxX3I09yI/As4PQk9wCvAZ6VZCdQwBHgpWPMKEk6iRWLvKquXGL1O8eQRZI0BL/ZuUYP3nffpCNI2uQ81ookNc4il6TGWeSS1DiLXJIa54ed0oj9/rv3n7DuVVftnEASbRaOyCWpcY7IO55EQlKrHJFLUuMscklqnFMrJ3HkxhsnHUGSVuSIXJIa54hcmhB3U9SoOCKXpMZZ5JLUOKdWxmCpfdJP/dEfnUASSZtBnzME3QA8HzhWVed3654IfADYzuAMQT9fVf8xvpjS9FpqrltaT32mVt4FXLpo3XXA7VW1A7i9W5YkTcCKRV5VnwK+vmj15cDe7vZe4IUjziVJ6mnYOfIzq+ooQFUdTXLGchsm2Q3sBjjvvPOGfDmpbU6/aJzGvtdKVe2pqtmqmp2ZmRn3y0nSpjNskd+fZBtAd31sdJEkSasx7NTKLcAu4PXd9c0jSyRNMadINI1WHJEnuRH4LPDUJPckeQmDAr8kyWHgkm5ZkjQBK47Iq+rKZe66eMRZJElD8JudY7DU4W+f5jc7JY3JpixyT+smaSPxoFmS1DiLXJIaZ5FLUuMscklqnEUuSY2zyCWpcRa5JDVuU+5HLi3mGe3VMkfkktQ4i1ySGufUirSMaThkrVM+6sMil6bINPzyUHucWpGkxq1pRJ7kCPAA8BDwYFXNjiKUJKm/UUytPLuqvjaC55E2rHseePiEdec83j+INRr+JElS49Za5AV8LMm+JLtHEUiStDprnVq5qKruS3IGcFuSu6vqUws36Ap+N8B55523xpfbOJY6S9Gpng5uU3G6RaOypp+aqrqvuz4GfAi4cIlt9lTVbFXNzszMrOXlpA3vngcefsRF6mPoEXmSxwKPqqoHuts/BfzOyJJJEzbsiNkC1npby9TKmcCHkhx/nvdV1UdGkmoD8oTPG4PTIZpGQxd5Vf0z8CMjzCKtm8/d++CkI0gj41f0p4gfgKoPj7+ixSxybUp/cfOBRyw7PaKWWeRq1lLTI888e2P9SI97Tn4zvIebgcMQSWqcv3ol2t9lsM+8+Vo+4F38WEft08V/DU3cuP+8f9mffWFkzyVNI4tcWqP1Hs27L7sWs8ilDWrYqRT3sW+PRa6xWu+51c06jdL6HL/WxiLXVHJUON3cbXG6+M5LGps+v5D9BbB2voNTblq/tu+IWaPi6H7tNvy75VEH187SVh/+nEzOhi/yjWhaR+nrbfHxUvRIi9+fF1/+wxNKonGzyDeRvn/CjnNk5ahNGj2LfINYPErfjCN0adzz7dM6n7+mBEkuBf4Y2AK8o6peP5JUE3DkxhsnHWEiHCFvHktNRU3rdMs0/PXYkrWcs3ML8KfAJcA9wN8nuaWq/mFU4YYxrR9uLv5Fsf3KK0f6/HMPn/GI5a0b8AfcOfHRa6nch7WWXwrTMNruYy0pLwS+3J3yjSTvBy4HJlrkrVhcvMuZfdSxoR+72P/ed7TXdlvP2jbU84+SpT1dpuGD01GOvjfaSD5VNdwDk58FLq2qX+6WrwJ+rKquXrTdbmB3t/hU4B+HjztWpwNfm3SINWg5f8vZoe38LWeHtvOvJvuTqmpmuTvXMiLPEutO+K1QVXuAPWt4nXWRZK6qZiedY1gt5285O7Sdv+Xs0Hb+UWZfy7Ev7wHOXbB8DnDf2uJIklZrLUX+98COJD+Q5HuAFwG3jCaWJKmvoadWqurBJFcDH2Ww++ENVXXXyJKtv6mf/llBy/lbzg5t5285O7Sdf2TZh/6wU5I0HTw/lCQ1ziKXpMZt2iJP8rtJDiTZn+RjSc5aZrsjSb7UbTe33jmX0jd7t+2WJF9Mcut6ZjyZPvmTPDrJ55PcmeSuJK+bRNal9Mx/bpJPJDnU5b9mElkXW8XP/Q1JjiU5uN4ZT2YV+S9N8o9JvpzkuvXOuZQkf5Tk7i7/h5I8YZntrklysPu5eXmvJ6+qTXkBvnfB7V8H3rbMdkeA0yedd5js3f3XAu8Dbp107tXkZ/A9hcd1t7cCdwDPnHT2VeTfBjyju/144J+Ap7eQvbvvJ4FnAAcnnXmI934L8BXgycD3AHdOyXv/U8Ap3e03AG9YYpvzgYPAYxjsjPJxYMdKz71pR+RV9Y0Fi49liS8zTau+2ZOcAzwPeMd65OqrT/4a+Ga3uLW7TMW/Uc/8R6vqC93tB4BDwNnrk3B5fX92qupTwNfXJdQq9Mz//4cPqarvAMcPHzJRVfWxqjp+bIDPMfjuzWJPAz5XVd/qtv1b4IqVnruNI8KMSZLfA14M/Bfw7GU2K+BjSQr48xp8U3XiemZ/C/CbDEaEU6VP/u7AbPuAHwT+tKruWL+EJ9fz/T++7XbgAgZ/VUzcarJPox75zwb+dcHyPcCPrUO01fgl4ANLrD8I/F6S7wO+DVwGrDilu6FH5Ek+3s01Lb5cDlBVr66qc4H3Alcv8zQXVdUzgJ8Gfi3JT7aQPcnzgWNVNZHDQY7iva+qh6pqJ4ORy4VJzm8pf/c8jwM+CLx80Why6rNPygjy9zp8yDislL3b5tXAgwzyPzJk1SEG0y63AR9hMC208hG+Jj1vNA0X4En0mAsEXgv8xqTz9skO/AGDkcgR4N+AbwHvmXTeNbz3r5m2936l/Aymgz4KXDvpnMO898D2Pv8205Yf+HHgowuWXwm8ctJ5uyy7gM8Cj+m5/e8DL1tpuw09Ij+ZJDsWLL4AuHuJbR6b5PHHbzP4sGLin+L3yV5Vr6yqc6pqO4PDJ/xNVf3COkU8qZ7v/czxT/WTnAo8d6ntJqFn/gDvBA5V1fXrlW0lfbJPs575p/LwIRmciOe3gBdU1bdOst0Z3fV5wM8AK5/1ZtK/oSb4m/GDDEr5APBXwNnd+rOAD3e3n8zgT5s7gbuAV086d9/si7Z/FtO110qf9/6HgS922xwEfnvSuVeZ/ycY/Dl/ANjfXS5rIXu3fCNwFPhfBn/ZvWTS2VeZ/zIGewp9ZYr+336Zwdz98Z+Hty2T/dMMzutwJ3Bxn+f2K/qS1LhNO7UiSRuFRS5JjbPIJalxFrkkNc4il6TGWeSS1DiLXFogye9OOoO0Wha51Eny/WzyA8mpTRa59F0XMPjGndQUi1ybVpJdSfZ1Z2z5NLATi1wN8s9IbUrdwdB+C9hZVd/pDtD1JuDwZJNJq+exVrQpJXkMgwOh3QrsraqpOB+rNAynVrQp1eAwoucDfwfsSfKyCUeShubUijalJDuq6jDw/iRPBx496UzSsJxa0aaU5F0MziTz3wymWH4FmAF+m8G5ID9SVR+fWEBpFRyRa1Oqql9cvC7JDwHfAf6kqr667qGkITkilxboyvwPgV+tqnsnnUfqwxG51EnyBmAL8FXg2ITjSL05Ipekxrn7oSQ1ziKXpMZZ5JLUOItckhpnkUtS4yxySWqcRS5JjbPIJalxFrkkNe7/AHAkkrFV865kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin = -3.55\n",
    "xmax = -2.9\n",
    "bins = np.arange(xmin, xmax, 0.01)\n",
    "x = plt.hist(solid_entropy, bins=bins, density=True, alpha=0.5, color=\"#EF9A9A\")\n",
    "x = plt.hist(solid_avg_entropy, bins=bins, density=True, alpha=0.5, color=\"#B71C1C\")\n",
    "x = plt.hist(liquid_entropy, bins=bins, density=True, alpha=0.5, color=\"#90CAF9\")\n",
    "x = plt.hist(liquid_avg_entropy, bins=bins, density=True, alpha=0.5, color=\"#0D47A1\")\n",
    "plt.xlabel(r\"$s_s^i$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distributions of $s_s^i$ given in light red and light blue are fairly distinct but show some overlap. The averaged entropy parameter, $\\bar{s}_s^i$ show distinct peaks which can distinguish solid and liquid very well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.47074287018929, -3.47074287018929)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(solid_entropy), np.mean(solid_avg_entropy)"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}