{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aca27212-c841-4b41-8a63-9767690291af",
   "metadata": {},
   "source": [
    "# Investigation of an Arrhenius relationship\n",
    "\n",
    "In addition to being able to determine the mean-squared displacement and diffusion coefficient from a given simulation, `kinisi` also includes [tools](./arrhenius.html) to investigate Arrhenius relationships. \n",
    "In this tutorial, we will look at how we can take advantage of these tools to study short, approximately 50 ps, simulations of lithium lanthanum zirconium oxide (LLZO).\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "Warning\n",
    "\n",
    "The warnings that are being ignored are related to the parsing of the files by `MDAnalysis` and lead to unnecessary print out to the screen that we want to avoid in the web documentation.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "980851d1-ed3a-4bbd-96fd-459b0c35a13d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import MDAnalysis as mda\n",
    "import matplotlib.pyplot as plt\n",
    "from kinisi.analyze import DiffusionAnalyzer\n",
    "from kinisi.arrhenius import StandardArrhenius\n",
    "import warnings\n",
    "np.random.seed(42)\n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f71efc1-d673-4c23-b300-779c276e9e37",
   "metadata": {},
   "source": [
    "To investigate this we will loop through a series of four temperatures and append each diffusion coefficient to a list. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d75aafe-222f-4042-b866-23c1ae623cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "temperatures = np.array([500, 600, 700, 800])\n",
    "D = []\n",
    "analyzers = []"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974c97c3-e678-48a1-a96a-3d1f971c11ee",
   "metadata": {},
   "source": [
    "To read these simulations we will use [MDAnalysis](https://userguide.mdanalysis.org/stable/index.html) (however, it is also possible to use data from a [VASP simulation](./vasp_d.html)).\n",
    "The parser, bootstrap, and diffusion parameters are all defined for all simulations, here we only consider the diffusive regime to begin after 5 ps.\n",
    "Additionally, we include in the `p_params` a `sub_sample_atoms` key, this defines the sampling frequency of atoms to be used in the analysis and the `sub_sample_traj` key, which defined the sampling frequency for the trajectory.\n",
    "This facility can be particularly useful for large simulations where `kinisi` might encounter issues related to out-of-memory problems. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6eb8f8c-1d22-4437-aa2e-322f0607fe9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_params = {'specie': 'LI',\n",
    "            'time_step': 5.079648e-4,\n",
    "            'step_skip': 100,\n",
    "            'min_dt': 0.001,\n",
    "            'sub_sample_atoms': 2,\n",
    "            'sub_sample_traj': 2,\n",
    "            'progress': False}\n",
    "d_params = {'progress': False}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f059cb7-8710-4f3f-a14c-9e6b2038eea1",
   "metadata": {},
   "source": [
    "File parsing and diffusion determination is then performed in a loop here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea84799-d980-4ed3-8b65-3eff77a5cab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in temperatures:\n",
    "    u = mda.Universe(f'_static/traj_{t}.gro', f'_static/traj_{t}.xtc')\n",
    "    d = DiffusionAnalyzer.from_universe(u, p_params)\n",
    "    d.diffusion(10, d_params)\n",
    "    D.append(d.D)\n",
    "    analyzers.append(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c204d962-a8f8-4989-9d84-e9a6ae52730f",
   "metadata": {},
   "source": [
    "The list of diffusion coefficient objects (which are `uravu.distribution.Distribution` type objects) and array of temperatures can then be passed to the `kinisi.arrhenius.StandardArrhenius` class, where we use the `bounds` keyword argument to give a minimum and maximum value for the activation energy and preexponential factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b9cfb47-30f1-4da6-b3c6-a2a70ef5e80c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = StandardArrhenius(temperatures, D, bounds=((0.01, 0.2), (1e-5, 1e-1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "507098fb-3afe-40b7-977d-1db0d8eb20a4",
   "metadata": {},
   "source": [
    "Having created the object, we can determine the maximum likelihood values for the parameters of activation energy and the preexponential factor. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d75002d-8fe9-453a-b593-62666fb51255",
   "metadata": {},
   "outputs": [],
   "source": [
    "s.max_likelihood('mini')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb2a858-0ed6-436b-a650-bfba1bf4a8b2",
   "metadata": {},
   "source": [
    "After determining the maximum likelihood values, we can use Markov chain Monte Carlo (MCMC) to sample the probability distributions for these. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9651df37-a7ea-4197-a69a-d842a66071ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "s.mcmc()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7aa14fcf-acbc-4515-b654-7b16b9710591",
   "metadata": {},
   "source": [
    "We can then visualise the probability distributions for the parameters as histograms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8ddacad-2cf4-4fe0-9715-385d176d3e58",
   "metadata": {},
   "outputs": [],
   "source": [
    "from corner import corner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ff6962-2e0c-4099-ad4c-500654e33406",
   "metadata": {},
   "outputs": [],
   "source": [
    "corner(s.flatchain, labels=['$E_a$/eV', '$A$/eV'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7240f6b8-cbbc-4826-b7ee-299465739943",
   "metadata": {},
   "source": [
    "It is also possible to plot these probability distributions as Arrhenius relations on the data measured values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afbbfa87-d83b-4686-bb2f-ae9c692dc7f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "credible_intervals = [[16, 84], [2.5, 97.5], [0.15, 99.85]]\n",
    "alpha = [0.6, 0.4, 0.2]\n",
    "\n",
    "plt.errorbar(1000/s.x, s.y.n, s.y.ci(), marker='o', ls='', color='k', zorder=10)\n",
    "for i, ci in enumerate(credible_intervals):\n",
    "    plt.fill_between(1000/s.x,\n",
    "                     *np.percentile(s.distribution, ci, axis=1),\n",
    "                     alpha=alpha[i],\n",
    "                     color='#0173B2',\n",
    "                     lw=0)\n",
    "plt.yscale('log')\n",
    "plt.xlabel('$1000T^{-1}$/K$^{-1}$')\n",
    "plt.ylabel('$D$/cm$^2$s$^{-1}$')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:kinisi]",
   "language": "python",
   "name": "conda-env-kinisi-py"
  },
  "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}