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.

[1]:
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).

[2]:
p_params = {'specie': 'Li',
            'time_step': 2.0,
            'step_skip': 50,
            'progress': False
            }

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.

[3]:
u_params = {'dimension': 'xy',
            'progress': False}
[4]:
xd = Xdatcar('./example_XDATCAR.gz')
diff = JumpDiffusionAnalyzer.from_Xdatcar(xd, parser_params=p_params, uncertainty_params=u_params)
/home/docs/checkouts/readthedocs.org/user_builds/kinisi/envs/latest/lib/python3.12/site-packages/kinisi/parser.py:116: UserWarning: Converting triclinic cell to orthorhombic this may have unexpected results. Triclinic a, b, c are not equilivalent to orthorhombic x, y, z.
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/kinisi/envs/latest/lib/python3.12/site-packages/kinisi/parser.py:438: UserWarning: UserWarning: Be aware that the expected input unit for 'time_step' is femtoseconds, not picoseconds.
  warnings.warn("UserWarning: Be aware that the expected input unit for 'time_step' is femtoseconds, not picoseconds.")

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.

[5]:
import matplotlib.pyplot as plt
[6]:
plt.errorbar(diff.dt, diff.mstd, diff.mstd_std)
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\Delta t$/ps')
plt.show()
_images/vasp_dj_9_0.png

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.

[7]:
plt.errorbar(diff.dt, diff.mstd, diff.mstd_std)
plt.axvline(4, color='g')
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\Delta t$/ps')
plt.yscale('log')
plt.xscale('log')
plt.show()
_images/vasp_dj_11_0.png

We will set the start of the diffusive regime to be 4 ps using the argument of 4.

[8]:
diff.jump_diffusion(4, jump_diffusion_params={'random_state': rng, 'progress': False})
/home/docs/checkouts/readthedocs.org/user_builds/kinisi/envs/latest/lib/python3.12/site-packages/scipy/optimize/_numdiff.py:596: RuntimeWarning: invalid value encountered in subtract
  df = fun(x1) - f0
/home/docs/checkouts/readthedocs.org/user_builds/kinisi/envs/latest/lib/python3.12/site-packages/scipy/optimize/_numdiff.py:596: RuntimeWarning: invalid value encountered in subtract
  df = fun(x1) - f0

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,

[9]:
diff.D_J.n, diff.D_J.con_int()
[9]:
(np.float64(2.1617667582483697e-06), array([8.68747246e-08, 8.47800806e-06]))
[10]:
diff.intercept.n, diff.intercept.con_int()
[10]:
(np.float64(378.11306480852204), array([197.90257266, 557.02370067]))

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

[11]:
diff.D_J.samples
[11]:
array([1.32836868e-06, 4.34868171e-06, 1.64058755e-06, ...,
       7.26396559e-07, 8.91074296e-07, 1.13748772e-06], shape=(3200,))

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

[13]:
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.

[14]:
loaded_diff = JumpDiffusionAnalyzer.load('example.hdf')

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

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

plt.plot(loaded_diff.dt, loaded_diff.mstd, 'k-')
for i, ci in enumerate(credible_intervals):
    plt.fill_between(loaded_diff.dt,
                     *np.percentile(loaded_diff.distribution, ci, axis=1),
                     alpha=alpha[i],
                     color='#0173B2',
                     lw=0)
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\Delta t$/ps')
plt.show()
_images/vasp_dj_26_0.png

Additionally, we can visualise the distribution of the diffusion coefficient that has been determined.

[16]:
plt.hist(loaded_diff.D_J.samples, density=True)
plt.axvline(loaded_diff.D_J.n, c='k')
plt.xlabel('$D$/cm$^2$s$^{-1}$')
plt.ylabel('$p(D$/cm$^2$s$^{-1})$')
plt.show()
_images/vasp_dj_28_0.png

Or the joint probability distribution for the diffusion coefficient and intercept.

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