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,
            '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.

u_params = {'dimension': 'xy',
            'progress': False}

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.

import matplotlib.pyplot as plt

plt.errorbar(diff.dt, diff.mstd, diff.mstd_std)
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\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(r'$\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, '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,

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

(np.float64(2.1617667582483697e-06), array([8.68747246e-08, 8.47800806e-06]))

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

(np.float64(378.11306480852204), array([197.90257266, 557.02370067]))

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

diff.D_J.samples

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).

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]

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()

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

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()

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()