{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding the MSD uncertainties from a VASP file\n",
    "\n",
    "Using ``kinisi`` to obtain the mean-squared displacement and uncertainty in a VASP Xdatcar type file is straightforward and involves using the `DiffusionAnalyzer` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipp as sc\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pymatgen.io.vasp import Xdatcar\n",
    "from kinisi.analyze import DiffusionAnalyzer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, the `params` dictionary describes some details about our simulation, this is unpacked when passed to the relevant class method. \n",
    "You can see other keyword arguments for a given class method in the [relevant documentation](./analyze.html#kinisi.diffusion_analyzer.DiffusionAnalyzer). \n",
    "We outline the `params` included in this example below:\n",
    "- `specie`: the atomic/ionic species to calculate the MSD for\n",
    "- `time_step`: the simulation timestep in your molecular dynamics simulation, for this example, this XDATCAR is from a VASP *ab-initio* molecular dynamics simulation and the simulation timestep is 2 femtoseconds. Note, that we use the fs unit from the `scipp` library. \n",
    "- `step_skip`: the frequency with which the data was written in your simulation trajectory. This XDATCAR contains every 50th molecular dynamics step (`NBLOCK=50` in VASP parlance).\n",
    "- `progress`: is you want a progress meter to be printed to the screen, this is not necessary for this documentation example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'specie': 'Li',\n",
    "          'time_step': 2.0 * sc.Unit('fs'),\n",
    "          'step_skip': 50 * sc.Unit('dimensionless'),\n",
    "          'progress': False\n",
    "          }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, for this example, we have a simulation that had a timestep of 2 femtoseconds but we only stored the trajectory information every 50 steps and we want to investigate only the lithium motion. \n",
    "\n",
    "The next step is to read the data file in to the appropriate `analyzer`. \n",
    "Below, the XDATCAR file is read into a `pymatgen` `Xdatcar` object, called `xd`. \n",
    "A `DiffusionAnalyzer` object is then created using the `from_xdatcar` class method, passing also the `params` dictionary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd = Xdatcar('./example_XDATCAR.gz')\n",
    "diff = DiffusionAnalyzer.from_xdatcar(xd, **params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `DiffusionAnalyzer` will determine the mean-squared displacements and variances (using the variance rescaling approach detailed in the [methodology](./methodology.html)).\n",
    "The `diff.msd` object is a `scipp.DataArray`, which can be plotted as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.errorbar(diff.dt.values, diff.msd.values, np.sqrt(diff.msd.variances))\n",
    "ax.set_xlabel(f'Time / {diff.dt.unit}')\n",
    "ax.set_ylabel(f'MSD / {diff.msd.unit}')\n",
    "ax.set_xlim(0, None)\n",
    "ax.set_ylim(0, None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to delve deeper into the data, this can be achieved by investigating the `scipp.DataArray` directly. \n",
    "This can be interpreted as a HTML object, as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff.da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by giving the appropriate units to the `time_step` input, the units of the time interval (the duration in which the particle travelled) can be correctly calculated. "
   ]
  }
 ],
 "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": 4
}