{ "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", "
\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", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "980851d1-ed3a-4bbd-96fd-459b0c35a13d", "metadata": {}, "outputs": [], "source": [ "# This initialises the duecredit framework\n", "%env DUECREDIT_ENABLE=yes\n", "%env DUECREDIT_FILE=arrhenius.duecredit.p\n", "import duecredit\n", "duecredit.due.activate()\n", "\n", "import numpy as np\n", "import scipp as sc\n", "import MDAnalysis as mda\n", "import matplotlib.pyplot as plt\n", "from kinisi.analyze import DiffusionAnalyzer\n", "from kinisi.arrhenius import Arrhenius\n", "from scipy.stats import uniform\n", "import warnings\n", "np.random.seed(42)\n", "warnings.filterwarnings(\"ignore\", category=UserWarning)\n", "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)" ] }, { "cell_type": "markdown", "id": "3f71efc1-d673-4c23-b300-779c276e9e37", "metadata": {}, "source": [ "To investigate this we will loop through a series of four temperatures and store each mean and variance of the diffusion coefficient in a dictionary. " ] }, { "cell_type": "code", "execution_count": null, "id": "7d75aafe-222f-4042-b866-23c1ae623cf3", "metadata": {}, "outputs": [], "source": [ "temperatures = np.array([500, 600, 700, 800])\n", "D = {'mean': [], 'var': []}" ] }, { "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 any `kinisi.Parser`).\n", "The parameters are defined for all simulations, here we only consider the diffusive regime to begin after 10 ps." ] }, { "cell_type": "code", "execution_count": null, "id": "a6eb8f8c-1d22-4437-aa2e-322f0607fe9e", "metadata": {}, "outputs": [], "source": [ "params = {'specie': 'LI',\n", " 'time_step': 5.079648e-4 * sc.Unit('ps'),\n", " 'step_skip': 100 * sc.units.dimensionless,\n", " '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, **params)\n", " d.diffusion(10 * sc.Unit('ps'), progress=False)\n", " D['mean'].append(sc.mean(d.D).value)\n", " D['var'].append(sc.var(d.D, ddof=1).value)\n", " unit = sc.mean(d.D).unit" ] }, { "cell_type": "markdown", "id": "4b50aa36", "metadata": {}, "source": [ "We can then use the diffusion coefficients and their variances can be used along with the array of temperatures to create a `scipp.DataArray`. " ] }, { "cell_type": "code", "execution_count": null, "id": "af3fcfde", "metadata": {}, "outputs": [], "source": [ "td = sc.DataArray(\n", " data=sc.array(dims=['temperature'], \n", " values=D['mean'], \n", " variances=D['var'],\n", " unit=unit,),\n", " coords={\n", " 'temperature': sc.Variable(dims=['temperature'], \n", " values=temperatures, \n", " unit='K')\n", " })\n", "td" ] }, { "cell_type": "markdown", "id": "c204d962-a8f8-4989-9d84-e9a6ae52730f", "metadata": {}, "source": [ "With the `DataArray` created, we can now create the `Arrhenius` object. " ] }, { "cell_type": "code", "execution_count": null, "id": "5b9cfb47-30f1-4da6-b3c6-a2a70ef5e80c", "metadata": {}, "outputs": [], "source": [ "priors = [uniform(0.1, 0.1), uniform(1e-5, 9e-5)]\n", "\n", "s = Arrhenius(td, priors=priors)" ] }, { "cell_type": "markdown", "id": "507098fb-3afe-40b7-977d-1db0d8eb20a4", "metadata": {}, "source": [ "When you create a `TemperatureDependent` object, of which `Arrhenius` is a subclass, you can optionally provide priors as a list of `scipy.stats` objects. \n", "For example, if we wanted the activation energy to be restricted between 0.1 and 0.2 eV and the preexponential factor between 10-5 and 10-4 cm2/s, then we would use the following: \n", "\n", "```python\n", "priors = [uniform(0.1, 0.1), uniform(1e-5, 9e-5)]\n", "\n", "s = Arrhenius(td, priors=priors)\n", "```\n", "\n", "So there are two uniform distributions, one for the activation energy and another for the preexponential factor.\n", "The first value in `uniform` is `loc` (the minimum), and the second is `scale` (the width, i.e., max - min). \n", "We recommend looking at the `scipy.stats` documentation to clarify how to use these objects. \n", "If you do not set `priors`, they are automatically set to be uniform with a min and max of 0.5 and 1.5 times the maximum-likelihood values.\n", "\n", "The `Arrhenius` object automatically estimates the maximum-likelihood values for the activation energy and preexponential factor and stores them in an appropriate `sc.DataGroup`. " ] }, { "cell_type": "code", "execution_count": null, "id": "3d75002d-8fe9-453a-b593-62666fb51255", "metadata": {}, "outputs": [], "source": [ "s" ] }, { "cell_type": "markdown", "id": "8fb2a858-0ed6-436b-a650-bfba1bf4a8b2", "metadata": {}, "source": [ "However, `kinisi` also allows us to estimate the full distribution of these parameters with Markov chain Monte Carlo (MCMC) sampling. " ] }, { "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(np.array([i.values for i in s.flatchain.values()]).T, \n", " labels=[' / '.join([k, str(v.unit)]) for k, v in s.flatchain.items()])\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/td.coords['temperature'].values, \n", " td.data.values, \n", " np.sqrt(td.data.variances), \n", " marker='o', \n", " ls='', \n", " color='k', \n", " zorder=10)\n", "\n", "for i, ci in enumerate(credible_intervals):\n", " plt.fill_between(1000/td.coords['temperature'].values,\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()" ] }, { "cell_type": "markdown", "id": "91e79f94", "metadata": {}, "source": [ "Finally, the activation energy and preexponential factor are `kinisi.samples.Samples` objects and are stored in the `Arrhenius` object. " ] }, { "cell_type": "code", "execution_count": null, "id": "8e93e891", "metadata": {}, "outputs": [], "source": [ "s.activation_energy" ] }, { "cell_type": "code", "execution_count": null, "id": "d5fc4945", "metadata": {}, "outputs": [], "source": [ "s.preexponential_factor" ] }, { "cell_type": "code", "execution_count": null, "id": "36ff9170", "metadata": {}, "outputs": [], "source": [ "from duecredit import due\n", "\n", "due.dump()" ] } ], "metadata": { "kernelspec": { "display_name": "kinisi-docs", "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.12.15" } }, "nbformat": 4, "nbformat_minor": 5 }