{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating energy\n",
    "\n",
    "pyscal can also be used for quick energy calculations. The energy per atom can give an insight into the local environment. pyscal relies on the [python library interface](https://lammps.sandia.gov/doc/Python_library.html) of [LAMMPS](https://lammps.sandia.gov/) for calculation of energy. The python library interface, hence, is a requirement for energy calculation. The easiest way to set up the LAMMPS library is by installing LAMMPS through the conda package.\n",
    "```\n",
    "conda install -c conda-forge lammps\n",
    "```\n",
    "Alternatively, [this page](https://lammps.sandia.gov/doc/Python_head.html) provides information on how to compile manually.\n",
    "\n",
    "### Interatomic potentials\n",
    "\n",
    "An interatomic potential is also required for the calculation of energy. The potential can be of [any type that LAMMPS supports](https://lammps.sandia.gov/doc/pair_style.html). For this example, we will use an EAM potential for Mo which is provided in the file `Mo.set`.   \n",
    "\n",
    "We start by importing the necessary modules,  "
   ]
  },
  {
   "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": [
    "For this example, a [LAMMPS dump file](https://lammps.sandia.gov/doc/dump.html) will be read in and the energy will be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(\"conf.bcc.dump\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the energy can be calculated by,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.calculate_energy(species=['Mo'], pair_style='eam/alloy', \n",
    "                     pair_coeff='* * Mo.set Mo', mass=95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first keyword above is `species`, which specifies the atomic species. This is required for [ASE](https://wiki.fysik.dtu.dk/ase/) module which is used under the hood for convertion of files. `pair_style` species the type of potential used in LAMMPS. See documentation [here](https://lammps.sandia.gov/doc/pair_style.html). `pair_coeff` is another LAMMPS command which is documented well [here](https://lammps.sandia.gov/doc/pair_coeff.html). Also, the mass needs to be provided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the calculation is over, the energy can be accessed for each atom as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.743130736133679"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atoms[0].energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to find the energy averaged over the neighbors using the `averaged` keyword. However, a neighbor calculation should be done before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "sys.calculate_energy(species=['Mo'], pair_style='eam/alloy', \n",
    "                     pair_coeff='* * Mo.set Mo', mass=95, averaged=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.534395941639571"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atoms = sys.atoms\n",
    "atoms[0].avg_energy"
   ]
  },
  {
   "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 energy parameters for each of these configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.calculate_energy(species=['Al'], pair_style='eam/alloy', \n",
    "                     pair_coeff='* * Al.eam.fs Al', mass=26.98, averaged=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets gather the energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms\n",
    "solid_energy = [atom.energy for atom in atoms]\n",
    "solid_avg_energy = [atom.avg_energy for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can repeat the calculations for the liquid phase,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(\"../tests/conf.lqd.Al.dump\")\n",
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "sys.calculate_energy(species=['Al'], pair_style='eam/alloy', \n",
    "                     pair_coeff='* * Al.eam.fs Al', mass=26.98, averaged=True)\n",
    "atoms = sys.atoms\n",
    "liquid_energy = [atom.energy for atom in atoms]\n",
    "liquid_avg_energy = [atom.avg_energy for atom in atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Energy')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEGCAYAAAB8Ys7jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAS+klEQVR4nO3dbZCdZ33f8e+vlhk/gMd2WanCNhHpqE4ZGmy8BVPazAQFxk08lqatCWSAbepW01KICe0EAS8yGaapO8mQMNNAR2Mg24QIFAdGKskEq2syKR3isDKOsSNTpdQRxkJaCAQCLtTk3xfnlrVendXe+3DO6tr9fmZ27odznT3/S3v022uvcz+kqpAktedvrHcBkqSVMcAlqVEGuCQ1ygCXpEYZ4JLUqC3jfLHnPve5tWPHjnG+pCQ17+jRo1+tqomF+8ca4Dt27GB2dnacLylJzUvy58P2O4UiSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNGuuZmK06/p73PL2+821vW8dKJOksR+CS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDWqV4An+dkkjyR5OMmBJJckuTrJkSTHu+VVoy5WknTWkgGe5BrgZ4DJqnoRcBHwWmAfMFNVO4GZbluSNCZ9p1C2AJcm2QJcBjwB7Aamu8engT1rX54kaTFLBnhVfRn4ZeAEcBL4y6q6F9hWVSe7NieBrcOen2Rvktkks3Nzc2tXuSRtcn2mUK5iMNp+AfA84PIkr+/7AlW1v6omq2pyYmJi5ZVKkp6hzxTKjwH/p6rmqur/AR8D/gFwKsl2gG55enRlSpIW6hPgJ4Cbk1yWJMAu4BhwGJjq2kwBh0ZToiRpmCUvJ1tV9ye5B3gAeAr4HLAfeDZwMMkdDEL+9lEWKkl6pl7XA6+qnwd+fsHu7zIYjUuS1oFnYkpSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoPvfEvD7Jg/O+vpnkrUmuTnIkyfFuedU4CpYkDfS5K/0XquqGqroBuAn4DvBxYB8wU1U7gZluW5I0JsudQtkF/O+q+nMGd6qf7vZPA3vWsjBJ0vktN8BfCxzo1rdV1UmAbrl12BOS7E0ym2R2bm5u5ZVKkp6hd4AneRZwG/Dby3mBqtpfVZNVNTkxMbHc+iRJi1jOCPwfAw9U1alu+1SS7QDd8vRaFydJWtxyAvx1nJ0+ATgMTHXrU8ChtSpKkrS0XgGe5DLgVcDH5u2+C3hVkuPdY3etfXmSpMVs6dOoqr4D/M0F+77G4KgUSdI68ExMSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1Kj+t6R58ok9yR5NMmxJC9PcnWSI0mOd8urRl2sJOmsviPw9wK/X1U/BLwYOAbsA2aqaicw021LksZkyQBPcgXwI8AHAKrqe1X1DWA3MN01mwb2jKpISdK5+ozAfxCYAz6U5HNJ7k5yObCtqk4CdMutw56cZG+S2SSzc3Nza1a4JG12fQJ8C/AS4P1VdSPwbZYxXVJV+6tqsqomJyYmVlimJGmhPgH+OPB4Vd3fbd/DINBPJdkO0C1Pj6ZESdIwSwZ4VX0F+FKS67tdu4A/BQ4DU92+KeDQSCqUJA21pWe7twAfTvIs4IvATzMI/4NJ7gBOALePpkRJ0jC9AryqHgQmhzy0a23LkST15ZmYktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNarX9cCTPAZ8C/g+8FRVTSa5GvgosAN4DHhNVX19NGVKkhZazgj8R6vqhqo6c2OHfcBMVe0EZljGjY4lSau3mimU3cB0tz4N7Fl9OZKkvvoGeAH3JjmaZG+3b1tVnQTolluHPTHJ3iSzSWbn5uZWX7EkCeh/U+NXVNUTSbYCR5I82vcFqmo/sB9gcnKyVlCjJGmIXiPwqnqiW54GPg68FDiVZDtAtzw9qiIlSedaMsCTXJ7kOWfWgVcDDwOHgamu2RRwaFRFSpLO1WcKZRvw8SRn2v9WVf1+ks8CB5PcAZwAbh9dmZKkhZYM8Kr6IvDiIfu/BuwaRVGSpKV5JqYkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVG9AzzJRUk+l+QT3fbVSY4kOd4trxpdmZKkhZYzAr8TODZvex8wU1U7gZluW5I0Jr0CPMm1wE8Ad8/bvRuY7tangT1rW5ok6Xz6jsB/Ffg54K/n7dtWVScBuuXWYU9MsjfJbJLZubm5VRUrSTpryQBPcitwuqqOruQFqmp/VU1W1eTExMRKvoUkaYgl70oPvAK4LcmPA5cAVyT5TeBUku1VdTLJduD0KAuVJD3TkiPwqnpHVV1bVTuA1wL3VdXrgcPAVNdsCjg0siolSefoMwJfzF3AwSR3ACeA29empAvPU088sd4lSNI5lhXgVfUHwB90618Ddq19SZKkPjwTU5IatZoplA3tyaMrOuhGksbGEbgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1Kj+twT85Ikf5zkT5I8kuQXuv1XJzmS5Hi3vGr05UqSzugzAv8u8MqqejFwA3BLkpuBfcBMVe0EZrptSdKY9LknZlXVX3WbF3dfBewGprv908CekVQoSRqq1xx4kouSPMjgzvNHqup+YFtVnQTollsXee7eJLNJZufm5taqbkna9HoFeFV9v6puAK4FXprkRX1foKr2V9VkVU1OTEystE5J0gLLOgqlqr7B4KbGtwCnkmwH6Jan17w6SdKi+hyFMpHkym79UuDHgEeBw8BU12wKODSqIiVJ5+pzU+PtwHSSixgE/sGq+kSSzwAHk9wBnABuH2GdkqQFlgzwqnoIuHHI/q8Bu0ZRlCRpaZ6JKUmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqVJ9bql2X5FNJjiV5JMmd3f6rkxxJcrxbXjX6ciVJZ/QZgT8F/Luq+rvAzcC/TfJCYB8wU1U7gZluW5I0JksGeFWdrKoHuvVvAceAa4DdwHTXbBrYM6oiJUnnWtYceJIdDO6PeT+wrapOwiDkga2LPGdvktkks3Nzc6urVpL0tN4BnuTZwO8Ab62qb/Z9XlXtr6rJqpqcmJhYSY2SpCF6BXiSixmE94er6mPd7lNJtnePbwdOj6ZESdIwfY5CCfAB4FhVvWfeQ4eBqW59Cji09uVJkhazpUebVwBvAD6f5MFu3zuBu4CDSe4ATgC3j6ZESdIwSwZ4VX0ayCIP71rbciRJfXkmpiQ1ygCXpEYZ4JLUKANckhplgEtSo/ocRqh5njx69On1S2+6aR0rkbTZGeCLeOzAgfUuQZLOyykUSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqOWPJEnyQeBW4HTVfWibt/VwEeBHcBjwGuq6uujK1Pq7xd/48Gn19/5hhvWsRJptPqMwH8duGXBvn3ATFXtBGa6bUnSGC0Z4FX1h8BfLNi9G5ju1qeBPWtclyRpCSu9Fsq2qjoJUFUnk2xdrGGSvcBegOc///krfDnpmVMj8zlNos1q5B9iVtX+qpqsqsmJiYlRv5wkbRorHYGfSrK9G31vB06vZVHSciw2Mpc2upWOwA8DU936FHBobcqRJPW1ZIAnOQB8Brg+yeNJ7gDuAl6V5Djwqm5bkjRGS06hVNXrFnlo1xrXIo2Ux4dro/FMTElqlLdU04bmB5zayByBS1KjHIHrguYIWlqcAa5Nzw831SoDXJuSI3ttBM6BS1KjHIHP8+TRo+tdglbg8W/99TO2r32O4xJtDga4Ljh/9OWnhu6fH9SGtGSAax3ND+qbrzn7Vvyvhx56en2xoF446pY2IwN8FeZPuVx6003rWEmb5gf1fWMeUS82mveIFLXEANeaWWxEvVib+TbaiLrPv4W0Wr6zNNRiQQv9Aul8z1+OlQT7KH4ZvOl9Dzy9/r43vWTF32fhv4vhrtXw3aNNbzUfjq52pO1IXavhO0aaZ36Yzx91z7dWf11Iq2WAb3IrCaPVBNj8Dy5bNb8Pb9z9w2v2fR2Na7lW9S5JcgvwXuAi4O6q8s48F6j1HDVuhNBeT4sFu4GvFf/Uk1wE/BqDW6o9Dnw2yeGq+tO1Kk79+Cf9hWH+z+EZv7RWMEpfzc90sef2OTJotb8I/KUyXqv5F34p8GdV9UWAJB8BdgMbOsAfO3Dg6fUdrzt7t7nFjglfybHirQbyZhxpL7fPC9svdwpmsfdGn/fMStqsZsS/3MNKR/XLo0+/W/1lk6pa2ROTfwbcUlX/stt+A/CyqnrzgnZ7gb3d5vXAF1Ze7rp5LvDV9S5izDZjn2Fz9nsz9hna6vcPVNXEwp2r+bWTIfvO+W1QVfuB/at4nXWXZLaqJte7jnHajH2Gzdnvzdhn2Bj9Xs35y48D183bvhZ4YnXlSJL6Wk2AfxbYmeQFSZ4FvBY4vDZlSZKWsuIplKp6KsmbgU8yOIzwg1X1yJpVdmFpegpohTZjn2Fz9nsz9hk2QL9X/CGmJGl9eVV8SWqUAS5JjTLAh0jy7iQPJXkwyb1Jnneethcl+VyST4yzxrXWp89JrkvyqSTHkjyS5M71qHUt9f1ZJ7klyReS/FmSfeOucy0l+aUkj3b9/niSKxdp97Pdz/nhJAeSXDLuWtfSMvp9ZZJ7urbHkrx83LX25Rz4EEmuqKpvdus/A7ywqv71Im3fBkwCV1TVrWMsc0316XOS7cD2qnogyXOAo8Celi+f0LPfFwH/i3mXjQBe12q/k7wauK87EOE/AVTV2xe0uQb4NIN/jyeTHAR+r6p+fewFr5E+/e7aTQP/o6ru7o6wu6yqvjHmcntxBD7Emf/QncsZcoISQJJrgZ8A7h5HXaPUp89VdbKqHujWvwUcA64ZT4Wj0fNn/fRlI6rqe8CZy0Y0qaruraoz55f/EYNzOIbZAlyaZAtwGY2f59Gn30muAH4E+ED3nO9dqOENXk52UUn+A/BG4C+BH12k2a8CPwc8Z1x1jVLPPp9puwO4Ebh/5IWNWI9+XwN8ad7248DLxlDaOPwL4KMLd1bVl5P8MnACeBK4t6ruHXdxIzS038APAnPAh5K8mMFfmXdW1bfHWVxfm3YEnuS/d3N7C792A1TVu6rqOuDDwJuHPP9W4HRVHV342IVqtX2e932eDfwO8NYFI9gL0hr0u9dlIy4kS/W5a/Mu4CkG/V74/KsY/JXxAuB5wOVJXj+u+ldqtf1mMKh9CfD+qroR+DZw4X7mUVV+necL+AHg4SH7/yODkdhjwFeA7wC/ud71jrLP3WMXMzh5623rXecYf9YvBz45b/sdwDvWu95V9nUK+AyD+d1hj98OfGDe9huB96133WPo998CHpu3/Y+A313vuhf72rQj8PNJsnPe5m3AowvbVNU7quraqtrB4DIC91XVBT9CWUyfPicJg7nBY1X1nnHVNkp9+s0Gu2xEBjdieTtwW1V9Z5FmJ4Cbk1zW/dx3MfjMo1l9+l1VXwG+lOT6btcuLuBLZBvgw93V/dn1EPBq4E6AJM9L8nvrW9rI9OnzK4A3AK/sDrt7MMmPr1O9a2XJftfgg68zl404Bhysti8b8Z8ZfG5zpPsZ/hc4p8/3A/cADwCfZ5AVrZ96vmS/O28BPty9J24AfnH8pfbjYYSS1ChH4JLUKANckhplgEtSowxwSWqUAS5JjfJUem0ISb7P4HC3Mz5SVXetVz3SOHgYoTaEJH9VVc9e4++5pc5e/Ei64DiFog0tyWNJfiHJA0k+n+SHuv2XJ/lgks9mcD333d3+f57kt5P8N+De7kzEgxlcQ/qjSe5PMpnkjiS/Mu91/lWSDXF2qtphgGujuHTe2aEPJvnJeY99tapeArwf+PfdvncxuPzB32dwBcJfSnJ599jLgamqeiXwJuDrVfXDwLuBm7o2HwFuS3Jxt/3TwIdG1jtpCOfAtVE8WVU3LPLYx7rlUeCfdOuvZhDAZwL9EuD53fqRqvqLbv0fAu8FqKozp9xTVd9Och9wa5JjwMVVNX8OXho5A1ybwXe75fc5+54P8E+r6gvzGyZ5GYNLiDKv3WLuBt7J4AJYjr41dk6haLP6JPCW7kp7JLlxkXafBl7TtXkh8PfOPNBd8Ok64KeAAyOtVhrCANdGsXAOfKlDCN/N4NrmDyV5uNse5n3ARDd18nbgIQZ37jnjIPA/q+rrq6xfWjYPI5TOI4MbGl9cVf83yd8GZoC/U4N7Y5LkE8CvVNXMetapzck5cOn8LgM+1R1tEuDfVNX3klwJ/DHwJ4a31osjcElqlHPgktQoA1ySGmWAS1KjDHBJapQBLkmN+v/lo6rpqJXwqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin = -3.5\n",
    "xmax = -2.5\n",
    "bins = np.arange(xmin, xmax, 0.01)\n",
    "x = plt.hist(solid_energy, bins=bins, density=True, alpha=0.5, color=\"#EF9A9A\")\n",
    "x = plt.hist(solid_avg_energy, bins=bins, density=True, alpha=0.5, color=\"#B71C1C\")\n",
    "x = plt.hist(liquid_energy, bins=bins, density=True, alpha=0.5, color=\"#90CAF9\")\n",
    "x = plt.hist(liquid_avg_energy, bins=bins, density=True, alpha=0.5, color=\"#0D47A1\")\n",
    "plt.xlabel(r\"Energy\")"
   ]
  }
 ],
 "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
}