{
"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": [
"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",
"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": [
"s = Arrhenius(td, bounds=[[0.1 * sc.Unit('eV'), 0.2 * sc.Unit('eV')], \n",
" [1e-5 * td.data.unit, 1e-4 * td.data.unit]])"
]
},
{
"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 bounds as a list of lists. \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",
"Arrhenius(td, bounds=[[0.1 * sc.Unit('eV'), 0.2 * sc.Unit('eV')], \n",
" [1e-5 * td.data.unit, 1e-4 * td.data.unit]])\n",
"```\n",
"\n",
"If you do not set bounds, they are automatically set to 0.5 and 1.5 of the best fit values.\n",
"\n",
"The `Arrhenius` object automatically estimates the maximum likelihood values for the activation energy and preexponential factor and stores them is 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"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "kinisi-dev",
"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.11.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}