# Jump diffusion coefficient from a VASP file#

Previously, we looked at obtaining accurate estimates for the mean-squared displacement with kinisi. Here, we show using the same class to evaluate the diffusion coefficient, using the kinisi methodology.

:

import numpy as np
from kinisi.analyze import JumpDiffusionAnalyzer
from pymatgen.io.vasp import Xdatcar
np.random.seed(42)
rng = np.random.RandomState(42)


There the p_params dictionary describes details about the simulation, and are documented in the parser module (also see the MSD Notebook).

:

p_params = {'specie': 'Li',
'time_step': 2.0,
'step_skip': 50}


While the u_params dictionary describes the details of the uncertainty calculation process, these are documented in the diffusion module. Here, we indicate that we only want to investigate diffusion in the xy-plane of the system.

:

u_params = {'dimension': 'xy'}

:

xd = Xdatcar('./example_XDATCAR.gz')
diff = JumpDiffusionAnalyzer.from_Xdatcar(xd, parser_params=p_params, uncertainty_params=u_params)

Reading Trajectory: 100%|██████████| 140/140 [00:00<00:00, 4501.98it/s]
Getting Displacements: 100%|██████████| 100/100 [00:00<00:00, 17149.01it/s]
Finding Means and Variances: 100%|██████████| 50/50 [00:00<00:00, 177.79it/s]


In the above cells, we parse and determine the uncertainty on the mean-squared displacement as a function of the timestep. We should visualise this, to check that we are observing diffusion in our material and to determine the timescale at which this diffusion begins.

:

import matplotlib.pyplot as plt

:

plt.errorbar(diff.dt, diff.mstd, diff.mstd_std)
plt.ylabel('MSD/Å$^2$')
plt.xlabel('$\Delta t$/ps')
plt.show() The vertical green line indicates the start of the diffusive regime, based on the maximum of the non-Gaussian parameter. We can also visualise this on a log-log scale, which helps to reveal the diffusive regime.

:

plt.errorbar(diff.dt, diff.mstd, diff.mstd_std)
plt.axvline(4, color='g')
plt.ylabel('MSD/Å$^2$')
plt.xlabel('$\Delta t$/ps')
plt.yscale('log')
plt.xscale('log')
plt.show() We will set the start of the diffusive regime to be 4 ps using the argument of 4.

:

diff.jump_diffusion(4, jump_diffusion_params={'random_state': rng})

Likelihood Sampling: 100%|██████████| 1500/1500 [00:01<00:00, 1267.54it/s]


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). Now we can probe these objects.

We can get the median and 95 % confidence interval using,

:

diff.D_J.n, diff.D_J.con_int()

:

(1.306269499961833e-05, array([1.33259472e-06, 3.09890243e-05]))

:

diff.intercept.n, diff.intercept.con_int()

:

(-153.82649129180245, array([-617.42572524,  209.28224073]))


Or we can get all of the sampled values from one of these objects.

:

diff.D_J.samples

:

array([1.89830541e-05, 9.34711468e-06, 1.95322989e-05, ...,
1.76636510e-05, 2.82372625e-05, 2.14902449e-06])


Having completed the analysis, we can save the object for use later (such as downstream by a plotting script).

:

diff.save('example.hdf')


The data is saved in a HDF5 file format which helps us efficiently store our data. We can then load the data with the following class method.

:

loaded_diff = JumpDiffusionAnalyzer.load('example.hdf')


We can plot the data with the credible intervals from the $$D$$ and $$D_{\text{offset}}$$ distribution.

:

credible_intervals = [[16, 84], [2.5, 97.5], [0.15, 99.85]]
alpha = [0.6, 0.4, 0.2]

for i, ci in enumerate(credible_intervals):
alpha=alpha[i],
color='#0173B2',
lw=0)
plt.ylabel('MSD/Å$^2$')
plt.xlabel('$\Delta t$/ps')
plt.show() Additionally, we can visualise the distribution of the diffusion coefficient that has been determined.

:

plt.hist(loaded_diff.D_J.samples, density=True)
plt.xlabel('$D$/cm$^2$s$^{-1}$')
plt.ylabel('$p(D$/cm$^2$s$^{-1})$')
plt.show() Or the joint probability distribution for the diffusion coefficient and intercept.

:

from corner import corner

:

corner(loaded_diff.flatchain, labels=['$D$/cm$^2$s$^{-1}$', 'intercept/Å$^2$'])
plt.show() [ ]: