diff --git a/docs/source/vasp_d.ipynb b/docs/source/vasp_d.ipynb index 8b14c07..8e816a1 100644 --- a/docs/source/vasp_d.ipynb +++ b/docs/source/vasp_d.ipynb @@ -9,7 +9,7 @@ "# Diffusion coefficient from a VASP file\n", "\n", "[Previously](./vasp_msd.html), we looked at obtaining accurate estimates for the mean-squared displacement with `kinisi`. \n", - "Here, we show using the same class to evaluate the diffusion coefficient, using the `kinisi` [methodology](./methodology.html)." + "Here, we show that the same `DiffusionAnalyzer` can be used to evaluate the diffusion coefficient, using the `kinisi` [methodology](./methodology.html)." ] }, { @@ -19,9 +19,10 @@ "outputs": [], "source": [ "import numpy as np\n", - "from kinisi.analyze import DiffusionAnalyzer\n", + "import matplotlib.pyplot as plt\n", + "import scipp as sc\n", "from pymatgen.io.vasp import Xdatcar\n", - "np.random.seed(42)\n", + "from kinisi.analyze import DiffusionAnalyzer\n", "rng = np.random.RandomState(42)" ] }, @@ -29,7 +30,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "There the `p_params` dictionary describes details about the simulation, and are documented in the [parser module](./parser.html) (also see the [MSD Notebook](./vasp_msd.html))." + "As wil the [previous example](./vasp_msd.html), `params` dictionary will describe the details about the simulation." ] }, { @@ -38,17 +39,19 @@ "metadata": {}, "outputs": [], "source": [ - "p_params = {'specie': 'Li',\n", - " 'time_step': 2.0,\n", - " 'step_skip': 50,\n", - " 'progress': False}" + "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": [ - "While the `b_params` dictionary describes the details of the bootstrapping process, these are documented in the [diffusion module](./diffusion.html#kinisi.diffusion.MSDBootstrap). Here, we indicate that we only want to investigate diffusion in the *xy*-plane of the system." + "In this example, we will add an additional key-value pair to the dictionary. \n", + "This argument means that the diffusion coefficient will only be calculated in the *xy* plane of the simulation box. " ] }, { @@ -57,35 +60,14 @@ "metadata": {}, "outputs": [], "source": [ - "u_params = {'dimension': 'xy', \n", - " 'progress': False}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "xd = Xdatcar('./example_XDATCAR.gz')\n", - "diff = DiffusionAnalyzer.from_Xdatcar(xd, parser_params=p_params, uncertainty_params=u_params)" + "params['dimension']= 'xy'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the above cells, we parse and determine the uncertainty on the mean-squared displacement as a function of the timestep. \n", - "We should visualise this, to check that we are observing diffusion in our material and to determine the timescale at which this diffusion begins. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt" + "As with the previous example, we now use the `from_xdatcar` class method to construct the `DiffusionAnalyzer` object. " ] }, { @@ -94,17 +76,16 @@ "metadata": {}, "outputs": [], "source": [ - "plt.errorbar(diff.dt, diff.msd, diff.msd_std)\n", - "plt.ylabel('MSD/Å$^2$')\n", - "plt.xlabel(r'$\\Delta t$/ps')\n", - "plt.show()" + "xd = Xdatcar('./example_XDATCAR.gz')\n", + "diff = DiffusionAnalyzer.from_xdatcar(xd, **params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can visualise this on a log-log scale, which helps to reveal the diffusive regime (the region where the gradient stops changing)." + "In the above cells, we parse and determine the uncertainty on the mean-squared displacement as a function of the timestep. \n", + "We should visualise this, to check that we are observing diffusion in our material and to determine the timescale at which this diffusion begins. " ] }, { @@ -113,12 +94,13 @@ "metadata": {}, "outputs": [], "source": [ - "plt.errorbar(diff.dt, diff.msd, diff.msd_std)\n", - "plt.axvline(2, color='g')\n", - "plt.ylabel('MSD/Å$^2$')\n", - "plt.xlabel(r'$\\Delta t$/ps')\n", - "plt.yscale('log')\n", - "plt.xscale('log')\n", + "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()" ] }, @@ -126,8 +108,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The green line at 2 ps appears to be a reasonable estimate of the start of the diffusive regime. \n", - "Therefore, we want to pass `2` as the argument to the diffusion analysis below. " + "We can visualise this on a log-log scale, which helps to reveal the diffusive regime (the region where the gradient stops changing)." ] }, { @@ -136,26 +117,24 @@ "metadata": {}, "outputs": [], "source": [ - "diff.diffusion(1, diffusion_params={'random_state': rng, 'progress': False})" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This method estimates the correlation matrix between the timesteps and uses likelihood sampling to find the diffusion coefficient, $D$ and intercept (the units are cm\\ :sup:`2`\\ /s and Å\\ :sup:`2` respectively).\n", - "Now we can probe these objects.\n", + "fig, ax = plt.subplots()\n", "\n", - "We can get the median and 95 % confidence interval using, " + "ax.errorbar(diff.dt.values, diff.msd.values, np.sqrt(diff.msd.variances))\n", + "ax.axvline(3000, color='g')\n", + "ax.set_xlabel(f'Time ({diff.dt.unit})')\n", + "ax.set_ylabel(f'MSD ({diff.msd.unit})')\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "plt.show()" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "diff.D.n, diff.D.ci()" + "The green line at 3000 fs appears to be a reasonable estimate of the start of the diffusive regime. \n", + "Therefore, we want to pass `3000 * sc.Units('fs')` as the argument to the diffusion analysis below. \n", + "At this stage, we pass the `random_state` argument to ensure reproducibility. " ] }, { @@ -164,14 +143,16 @@ "metadata": {}, "outputs": [], "source": [ - "diff.intercept.n, diff.intercept.ci()" + "start_of_diffusion = 3000 * sc.Unit('fs')\n", + "diff.diffusion(start_of_diffusion, progress=False, random_state=rng)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Or we can get all of the sampled values from one of these objects. " + "This method estimates the correlation matrix between the timesteps and uses posterior sampling to find the self-diffusion coefficient, $D*$ and intercept.\n", + "We can find the mean of the marginal posterior samples of $D*$: " ] }, { @@ -180,38 +161,14 @@ "metadata": {}, "outputs": [], "source": [ - "diff.D.samples" + "diff.D.mean()" ] }, { "cell_type": "markdown", - "metadata": { - "tags": [] - }, - "source": [ - "Having completed the analysis, we can save the object for use later (such as downstream by a plotting script). " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "nbsphinx": "hidden", - "tags": [] - }, - "source": [ - "These hidden cells exist to remove any existing `example.hdf` file on builiding." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "nbsphinx": "hidden", - "tags": [] - }, - "outputs": [], + "metadata": {}, "source": [ - "!rm example.hdf" + "And the standard deviation: " ] }, { @@ -220,15 +177,15 @@ "metadata": {}, "outputs": [], "source": [ - "diff.save('example.hdf')" + "diff.D.std(ddof=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The data is saved in a HDF5 file format which helps us efficiently store our data. \n", - "We can then load the data with the following class method. " + "The same for the intercept can be found from `diff.intercept`. \n", + "A histogram of the marginal posterior probability distribution for $D*$ can be plotted as shown below. " ] }, { @@ -237,14 +194,21 @@ "metadata": {}, "outputs": [], "source": [ - "loaded_diff = DiffusionAnalyzer.load('example.hdf')" + "fig, ax = plt.subplots()\n", + "\n", + "ax.hist(diff.D.values, bins=40, density=True)\n", + "ax.set_xlabel(f'D* ({diff.D.unit})')\n", + "ax.set_ylabel(f'p(D*) ({(1 / diff.D.unit).unit})')\n", + "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can plot the data with the credible intervals from the $D$ and $D_{\\text{offset}}$ distribution." + "It is also possible to plot the posterior distribution of the models on the data. \n", + "We represent this distribution with 1σ-, 2σ-, and 3σ-credible intervals. \n", + "Here we remove the measured error bars, in the interest of clarity." ] }, { @@ -256,35 +220,17 @@ "credible_intervals = [[16, 84], [2.5, 97.5], [0.15, 99.85]]\n", "alpha = [0.6, 0.4, 0.2]\n", "\n", - "plt.plot(loaded_diff.dt, loaded_diff.msd, 'k-')\n", + "fig, ax = plt.subplots()\n", + "ax.plot(diff.dt.values, diff.msd.values, 'k-')\n", "for i, ci in enumerate(credible_intervals):\n", - " plt.fill_between(loaded_diff.dt,\n", - " *np.percentile(loaded_diff.distribution, ci, axis=1),\n", - " alpha=alpha[i],\n", - " color='#0173B2',\n", - " lw=0)\n", - "plt.ylabel('MSD/Å$^2$')\n", - "plt.xlabel(r'$\\Delta t$/ps')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Additionally, we can visualise the distribution of the diffusion coefficient that has been determined." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.hist(loaded_diff.D.samples, density=True)\n", - "plt.axvline(loaded_diff.D.n, c='k')\n", - "plt.xlabel('$D$/cm$^2$s$^{-1}$')\n", - "plt.ylabel('$p(D$/cm$^2$s$^{-1})$')\n", + " ax.fill_between(diff.dt.values,\n", + " *np.percentile(diff.distributions, ci, axis=1),\n", + " alpha=alpha[i],\n", + " color='#0173B2',\n", + " lw=0)\n", + "ax.set_ylabel('MSD ' + ax.get_ylabel())\n", + "ax.set_xlim(0, None)\n", + "ax.set_ylim(0, None)\n", "plt.show()" ] }, @@ -292,7 +238,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Or the joint probability distribution for the diffusion coefficient and intercept." + "Finally, the joint posterior probability distribution for the diffusion coefficient and intercept can be visualised with the `corner` library. " ] }, { @@ -310,52 +256,15 @@ "metadata": {}, "outputs": [], "source": [ - "corner(loaded_diff._diff.flatchain, labels=['$D$/cm$^2$s$^{-1}$', 'intercept/Å$^2$'])\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Above, the posterior distribution is plotted as a function of $\\Delta t$.\n", - "`kinisi` can also compute the [posterior predictive distribution](https://en.wikipedia.org/wiki/Posterior_predictive_distribution) of mean-squared displacement data that would be expected from simulations, given the set of linear models for the _ensemble_ MSD described by the posterior distribution, above.\n", - "Calculating this posterior predictive distribution involves an additional sampling procedure that samples probable simulated MSDs for each of the probable linear models in the (previously sampled) posterior distribution.\n", - "The posterior predictive distribution captures uncertainty in the true ensemble MSD and uncertainty in the simulated MSD data expected for each ensemble MSD model, and so has a larger variance than the ensemble MSD posterior distribution." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ppd = loaded_diff.posterior_predictive({'n_posterior_samples': 128,\n", - " 'progress': False})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.plot(loaded_diff.dt, loaded_diff.msd, 'k-')\n", - "for i, ci in enumerate(credible_intervals):\n", - " plt.fill_between(loaded_diff.dt[7:],\n", - " *np.percentile(ppd.T, ci, axis=1),\n", - " alpha=alpha[i],\n", - " color='#0173B2',\n", - " lw=0)\n", - "plt.ylabel('MSD/Å$^2$')\n", - "plt.xlabel('$\\Delta t$/ps')\n", + "corner(np.array([i.values for i in diff.flatchain.values()]).T, \n", + " labels=['/'.join([k, str(v.unit)]) for k, v in diff.flatchain.items()])\n", "plt.show()" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "kinisi-scipp", "language": "python", "name": "python3" }, @@ -369,7 +278,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/docs/source/vasp_msd.ipynb b/docs/source/vasp_msd.ipynb index 274a6ff..31d2030 100644 --- a/docs/source/vasp_msd.ipynb +++ b/docs/source/vasp_msd.ipynb @@ -15,17 +15,24 @@ "metadata": {}, "outputs": [], "source": [ + "import scipp as sc\n", "import numpy as np\n", - "from kinisi.analyze import DiffusionAnalyzer\n", - "np.random.seed(1)" + "import matplotlib.pyplot as plt\n", + "from pymatgen.io.vasp import Xdatcar\n", + "from kinisi.analyze import DiffusionAnalyzer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "There the `params` dictionary describes details about the simulation, and are fully documented in the [parser module](./parser.html).\n", - "The `'specie'` is the atomic/ionic species that you want to calculate the MSD for, the `'time_step'` is the simulation timestep in your molecular dynamics, and the `'step_skip'` is the frequency which with the data was written in your simulation trajectory." + "Below, the `params` dictionary describes some details about our simulation, this are 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. " ] }, { @@ -35,18 +42,21 @@ "outputs": [], "source": [ "params = {'specie': 'Li',\n", - " 'time_step': 2.0,\n", - " 'step_skip': 50,\n", + " 'time_step': 2.0 * sc.Unit('fs'),\n", + " 'step_skip': 50 * sc.Unit('dimensionless'),\n", " 'progress': False\n", - " }\n", - "u_params = {'progress': False}" + " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Therefore, for this example, we have a simulation that had a timestep of 2 femtoseconds (note the [FAQs about units](https://kinisi.readthedocs.io/en/latest/faq.html)) but we only stored the trajectory information every 50 steps and we want to investigate only the lithium motion. " + "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. " ] }, { @@ -55,17 +65,16 @@ "metadata": {}, "outputs": [], "source": [ - "msd = DiffusionAnalyzer.from_file('example_XDATCAR.gz', parser_params=params, uncertainty_params=u_params)" + "xd = Xdatcar('./example_XDATCAR.gz')\n", + "diff = DiffusionAnalyzer.from_xdatcar(xd, **params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The `DiffusionAnalyzer` will perform the bootstrapping process to obtain the displacements and uncertainties (this is detailed in the [methodology](./methodology.html). \n", - "To find out how to determine the diffusion coefficient with `DiffusionAnalyzer` continue [here](https://kinisi.readthedocs.io/en/latest/vasp_d.html). \n", - "\n", - "Then the MSD (`msd`) as a function of timestep (`dt`) can be plotted." + "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. " ] }, { @@ -74,7 +83,22 @@ "metadata": {}, "outputs": [], "source": [ - "import matplotlib.pyplot as plt" + "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. " ] }, { @@ -83,16 +107,20 @@ "metadata": {}, "outputs": [], "source": [ - "plt.errorbar(msd.dt, msd.msd, msd.msd_std)\n", - "plt.ylabel('MSD/Å$^2$')\n", - "plt.xlabel(r'$\\Delta t$/ps')\n", - "plt.show()" + "diff.msd" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that by given 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": "Python 3 (ipykernel)", + "display_name": "kinisi-scipp", "language": "python", "name": "python3" }, @@ -106,7 +134,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/kinisi/analyzer.py b/kinisi/analyzer.py index b95104a..885c86f 100644 --- a/kinisi/analyzer.py +++ b/kinisi/analyzer.py @@ -108,7 +108,7 @@ def _from_universe(cls, distance_unit=distance_unit, progress=progress) return cls(p) - + @property def n_atoms(self) -> int: """ @@ -116,6 +116,20 @@ def n_atoms(self) -> int: """ return self.trajectory.displacements.sizes['atom'] + @property + def intercept(self) -> sc.Variable: + """ + :return: The intercept of the linear relationship. + """ + return self.diff.intercept + + @property + def dt(self) -> sc.Variable: + """ + :return: The time intervals used for the mean-squared displacement. + """ + return self.msd_da.coords['time interval'] + def _flatten_list(this_list: list) -> list: """ diff --git a/kinisi/conductivity_analyzer.py b/kinisi/conductivity_analyzer.py index 521b70c..29b351c 100644 --- a/kinisi/conductivity_analyzer.py +++ b/kinisi/conductivity_analyzer.py @@ -10,11 +10,13 @@ from typing import Union, List import numpy as np import scipp as sc +from scipp.typing import VariableLikeType from kinisi.diffusion import Diffusion from kinisi.displacement import calculate_mscd from kinisi.parser import Parser from kinisi.analyzer import Analyzer + class ConductivityAnalyzer(Analyzer): """ The :py:class:`kinisi.analyze.ConductivityAnalyzer` class performs analysis of conductivity in materials. @@ -28,17 +30,17 @@ class ConductivityAnalyzer(Analyzer): def __init__(self, trajectory: Parser) -> None: super().__init__(trajectory) - self.mscd = None + self.msd_da = None @classmethod - def from_xdatcar(cls, + def from_xdatcar(cls, trajectory: Union['pymatgen.io.vasp.outputs.Xdatcar', List['pymatgen.io.vasp.outputs.Xdatcar']], specie: Union['pymatgen.core.periodic_table.Element', 'pymatgen.core.periodic_table.Specie'], - time_step: sc.Variable, - step_skip: sc.Variable, - ionic_charge: sc.Variable, + time_step: VariableLikeType, + step_skip: VariableLikeType, + ionic_charge: VariableLikeType, dtype: Union[str, None] = None, - dt: sc.Variable = None, + dt: VariableLikeType = None, dimension: str = 'xyz', distance_unit: sc.Unit = sc.units.angstrom, progress: bool = True) -> 'ConductivityAnalyzer': @@ -72,11 +74,16 @@ def from_xdatcar(cls, :returns: The :py:class:`ConductivityAnalyzer` object with the mean-squared charge displacement calculated. """ - p = super()._from_xdatcar(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, progress) - p.mscd = calculate_mscd(p.trajectory, ionic_charge, progress) + p = super()._from_xdatcar(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, + progress) + p.msd_da = calculate_mscd(p.trajectory, ionic_charge, progress) return p - def conductivity(self, start_dt: sc.Variable, temperature: sc.Variable, volume: sc.Variable, diffusion_params: Union[dict, None] = None) -> None: + def conductivity(self, + start_dt: VariableLikeType, + temperature: VariableLikeType, + volume: VariableLikeType, + diffusion_params: Union[dict, None] = None) -> None: """ Calculation of the conductivity. Keyword arguments will be passed of the :py:func:`bayesian_regression` method. @@ -88,7 +95,7 @@ def conductivity(self, start_dt: sc.Variable, temperature: sc.Variable, volume: """ if diffusion_params is None: diffusion_params = {} - self.diff = Diffusion(msd=self.mscd, n_atoms=self.n_atoms) + self.diff = Diffusion(msd=self.msd_da, n_atoms=self.n_atoms) self.diff._conductivity(start_dt=start_dt, temperature=temperature, volume=volume, **diffusion_params) @property @@ -98,6 +105,31 @@ def distributions(self) -> np.array: plotting of credible intervals. """ if self.diff.intercept is not None: - return self.diff.gradient.values * self.mscd.coords['time interval'].values[:, np.newaxis] + self.diff.intercept.values + return self.diff.gradient.values * self.msd_da.coords[ + 'time interval'].values[:, np.newaxis] + self.diff.intercept.values else: - return self.diff.gradient.values * self.mscd.coords['time interval'].values[:, np.newaxis] \ No newline at end of file + return self.diff.gradient.values * self.msd_da.coords['time interval'].values[:, np.newaxis] + + @property + def sigma(self) -> VariableLikeType: + """ + :return: The conductivity. + """ + return self.diff.sigma + + @property + def mscd(self) -> VariableLikeType: + """ + :return: The mean-squared charge displacement. + """ + return self.msd_da.data + + @property + def flatchain(self) -> sc.DataGroup: + """ + :returns: The flatchain of the MCMC samples. + """ + flatchain = {'sigma': self.sigma} + if self.intercept is not None: + flatchain['intercept'] = self.intercept + return sc.DataGroup(**flatchain) diff --git a/kinisi/diffusion_analyzer.py b/kinisi/diffusion_analyzer.py index 76d3f9c..aeb0e4f 100644 --- a/kinisi/diffusion_analyzer.py +++ b/kinisi/diffusion_analyzer.py @@ -10,6 +10,7 @@ from typing import Union, List import numpy as np import scipp as sc +from scipp.typing import VariableLikeType from kinisi.displacement import calculate_msd from kinisi.diffusion import Diffusion from kinisi.parser import Parser, PymatgenParser, MDAnalysisParser @@ -26,16 +27,16 @@ class DiffusionAnalyzer(Analyzer): def __init__(self, trajectory: Parser) -> None: super().__init__(trajectory) - self.msd = None + self.msd_da = None @classmethod def from_xdatcar(cls, trajectory: Union['pymatgen.io.vasp.outputs.Xdatcar', List['pymatgen.io.vasp.outputs.Xdatcar']], specie: Union['pymatgen.core.periodic_table.Element', 'pymatgen.core.periodic_table.Specie'], - time_step: sc.Variable, - step_skip: sc.Variable, + time_step: VariableLikeType, + step_skip: VariableLikeType, dtype: Union[str, None] = None, - dt: sc.Variable = None, + dt: VariableLikeType = None, dimension: str = 'xyz', distance_unit: sc.Unit = sc.units.angstrom, progress: bool = True) -> 'DiffusionAnalyzer': @@ -69,17 +70,17 @@ def from_xdatcar(cls, """ p = super()._from_xdatcar(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, progress) - p.msd = calculate_msd(p.trajectory, progress) + p.msd_da = calculate_msd(p.trajectory, progress) return p @classmethod def from_universe(cls, trajectory: 'MDAnalysis.core.universe.Universe', specie: str = None, - time_step: sc.Variable = None, - step_skip: sc.Variable = None, + time_step: VariableLikeType = None, + step_skip: VariableLikeType = None, dtype: Union[str, None] = None, - dt: sc.Variable = None, + dt: VariableLikeType = None, dimension: str = 'xyz', distance_unit: sc.Unit = sc.units.angstrom, progress: bool = True) -> 'DiffusionAnalyzer': @@ -89,22 +90,44 @@ def from_universe(cls, :param trajectory: The :py:class:`MDAnalysis """ - p = super()._from_universe(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, progress) - p.msd = calculate_msd(p.trajectory, progress) + p = super()._from_universe(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, + progress) + p.msd_da = calculate_msd(p.trajectory, progress) return p - def diffusion(self, start_dt: sc.Variable, diffusion_params: Union[dict, None] = None) -> None: + def diffusion(self, + start_dt: VariableLikeType, + cond_max: float = 1e16, + fit_intercept: bool = True, + n_samples: int = 1000, + n_walkers: int = 32, + n_burn: int = 500, + n_thin: int = 10, + progress: bool = True, + random_state: np.random.mtrand.RandomState = None) -> None: """ Calculate the diffusion coefficient using the mean-squared displacement data. :param start_dt: The time at which the diffusion regime begins. - :param diffusion_params: The keyword arguements for the diffusion calculation - (see :py:func:`diffusion.bayesian_regression`). Optional, defaults to :py:attr:`None`. + :param cond_max: The maximum condition number of the covariance matrix. Optional, default is :py:attr:`1e16`. + :param fit_intercept: Whether to fit an intercept. Optional, default is :py:attr:`True`. + :param n_samples: The number of MCMC samples to take. Optional, default is :py:attr:`1000`. + :param n_walkers: The number of walkers to use in the MCMC. Optional, default is :py:attr:`32`. + :param n_burn: The number of burn-in samples to discard. Optional, default is :py:attr:`500`. + :param n_thin: The thinning factor for the MCMC samples. Optional, default is :py:attr:`10`. + :param progress: Whether to show the progress bar. Optional, default is :py:attr:`True`. + :param random_state: The random state to use for the MCMC. Optional, default is :py:attr:`None`. """ - if diffusion_params is None: - diffusion_params = {} - self.diff = Diffusion(self.msd) - self.diff._diffusion(start_dt, **diffusion_params) + self.diff = Diffusion(self.msd_da) + self.diff._diffusion(start_dt, + cond_max=cond_max, + fit_intercept=fit_intercept, + n_samples=n_samples, + n_walkers=n_walkers, + n_burn=n_burn, + n_thin=n_thin, + progress=progress, + random_state=random_state) @property def distributions(self) -> np.array: @@ -113,6 +136,31 @@ def distributions(self) -> np.array: plotting of credible intervals. """ if self.diff.intercept is not None: - return self.diff.gradient.values * self.msd.coords['time interval'].values[:, np.newaxis] + self.diff.intercept.values + return self.diff.gradient.values * self.msd_da.coords[ + 'time interval'].values[:, np.newaxis] + self.diff.intercept.values else: - return self.diff.gradient.values * self.msd.coords['time interval'].values[:, np.newaxis] \ No newline at end of file + return self.diff.gradient.values * self.msd_da.coords['time interval'].values[:, np.newaxis] + + @property + def D(self) -> VariableLikeType: + """ + :return: The self-diffusion coefficient. + """ + return self.diff.D + + @property + def msd(self) -> VariableLikeType: + """ + :return: The mean-squared displacement data. + """ + return self.msd_da.data + + @property + def flatchain(self) -> sc.DataGroup: + """ + :returns: The flatchain of the MCMC samples. + """ + flatchain = {'D*': self.D} + if self.intercept is not None: + flatchain['intercept'] = self.intercept + return sc.DataGroup(**flatchain) diff --git a/kinisi/jump_diffusion_analyzer.py b/kinisi/jump_diffusion_analyzer.py index 91e2dd7..a49dadb 100644 --- a/kinisi/jump_diffusion_analyzer.py +++ b/kinisi/jump_diffusion_analyzer.py @@ -10,6 +10,7 @@ from typing import Union, List import numpy as np import scipp as sc +from scipp.typing import VariableLikeType from kinisi.displacement import calculate_mstd from kinisi.diffusion import Diffusion from kinisi.parser import Parser, PymatgenParser @@ -31,16 +32,16 @@ class JumpDiffusionAnalyzer(Analyzer): def __init__(self, trajectory: Parser) -> None: super().__init__(trajectory) - self.mstd = None + self.msd_da = None @classmethod def from_xdatcar(cls, trajectory: Union['pymatgen.io.vasp.outputs.Xdatcar', List['pymatgen.io.vasp.outputs.Xdatcar']], specie: Union['pymatgen.core.periodic_table.Element', 'pymatgen.core.periodic_table.Specie'], - time_step: sc.Variable, - step_skip: sc.Variable, + time_step: VariableLikeType, + step_skip: VariableLikeType, dtype: Union[str, None] = None, - dt: sc.Variable = None, + dt: VariableLikeType = None, dimension: str = 'xyz', distance_unit: sc.Unit = sc.units.angstrom, progress: bool = True) -> 'JumpDiffusionAnalyzer': @@ -74,10 +75,10 @@ def from_xdatcar(cls, """ p = super()._from_xdatcar(trajectory, specie, time_step, step_skip, dtype, dt, dimension, distance_unit, progress) - p.mstd = calculate_mstd(p.trajectory, progress) + p.msd_da = calculate_mstd(p.trajectory, progress) return p - def jump_diffusion(self, start_dt: sc.Variable, diffusion_params: Union[dict, None] = None) -> None: + def jump_diffusion(self, start_dt: VariableLikeType, diffusion_params: Union[dict, None] = None) -> None: """ Calculate the jump diffusion coefficient using the mean-squared total displacement data. @@ -87,10 +88,9 @@ def jump_diffusion(self, start_dt: sc.Variable, diffusion_params: Union[dict, No """ if diffusion_params is None: diffusion_params = {} - self.diff = Diffusion(msd=self.mstd, n_atoms=self.n_atoms) + self.diff = Diffusion(msd=self.msd_da, n_atoms=self.n_atoms) self.diff._jump_diffusion(start_dt, **diffusion_params) - @property def distributions(self) -> np.array: """ @@ -98,6 +98,31 @@ def distributions(self) -> np.array: plotting of credible intervals. """ if self.diff.intercept is not None: - return self.diff.gradient.values * self.mstd.coords['time interval'].values[:, np.newaxis] + self.diff.intercept.values + return self.diff.gradient.values * self.msd_da.coords[ + 'time interval'].values[:, np.newaxis] + self.diff.intercept.values else: - return self.diff.gradient.values * self.mstd.coords['time interval'].values[:, np.newaxis] \ No newline at end of file + return self.diff.gradient.values * self.msd_da.coords['time interval'].values[:, np.newaxis] + + @property + def D_J(self) -> VariableLikeType: + """ + :return: The jump diffusion coefficient. + """ + return self.diff.D_J + + @property + def mstd(self) -> VariableLikeType: + """ + :return: The mean-squared total displacement. + """ + return self.msd_da.data + + @property + def flatchain(self) -> sc.DataGroup: + """ + :return: The flatchain of the MCMC samples. + """ + flatchain = {'D_J': self.D_J} + if self.intercept is not None: + flatchain['intercept'] = self.intercept + return sc.DataGroup(**flatchain)