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 DiffusionAnalyzer
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 b_params dictionary describes the details of the bootstrapping 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 = DiffusionAnalyzer.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:437: 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.msd, diff.msd_std)
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\Delta t$/ps')
plt.show()

We can visualise this on a log-log scale, which helps to reveal the diffusive regime (the region where the gradient stops changing).

plt.errorbar(diff.dt, diff.msd, diff.msd_std)
plt.axvline(2, color='g')
plt.ylabel('MSD/Å$^2$')
plt.xlabel(r'$\Delta t$/ps')
plt.yscale('log')
plt.xscale('log')
plt.show()

The green line at 2 ps appears to be a reasonable estimate of the start of the diffusive regime. Therefore, we want to pass 2 as the argument to the diffusion analysis below.

diff.diffusion(2, 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:592: 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.n, diff.D.ci()

(1.3328599142011123e-05, array([1.13998021e-05, 1.52771231e-05]))

diff.intercept.n, diff.intercept.ci()

(0.5643281864864429, array([0.46956222, 0.66207265]))

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

diff.D.samples

array([1.33939630e-05, 1.32179803e-05, 1.42601657e-05, ...,
       1.52350661e-05, 1.48898418e-05, 1.30407890e-05])

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 = DiffusionAnalyzer.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.msd, '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.samples, density=True)
plt.axvline(loaded_diff.D.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._diff.flatchain, labels=['$D$/cm$^2$s$^{-1}$', 'intercept/Å$^2$'])
plt.show()

Above, the posterior distribution is plotted as a function of \(\Delta t\). kinisi can also compute the 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. 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. 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.

ppd = loaded_diff.posterior_predictive({'n_posterior_samples': 128,
                                        'progress': False})

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

<>:9: SyntaxWarning: invalid escape sequence '\D'
<>:9: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_1202/137998712.py:9: SyntaxWarning: invalid escape sequence '\D'
  plt.xlabel('$\Delta t$/ps')