Investigation of an Arrhenius relationship

In addition to being able to determine the mean-squared displacement and diffusion coefficient from a given simulation, kinisi also includes tools to investigate Arrhenius relationships. In this tutorial, we will look at how we can take advantage of these tools to study short, approximately 50 ps, simulations of lithium lanthanum zirconium oxide (LLZO).

Warning

The warnings that are being ignored are related to the parsing of the files by MDAnalysis and lead to unnecessary print out to the screen that we want to avoid in the web documentation.

import numpy as np
import MDAnalysis as mda
import matplotlib.pyplot as plt
from kinisi.analyze import DiffusionAnalyzer
from kinisi.arrhenius import StandardArrhenius
import warnings
np.random.seed(42)
warnings.filterwarnings("ignore", category=UserWarning)

/home/docs/checkouts/readthedocs.org/user_builds/kinisi/envs/latest/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
Matplotlib is building the font cache; this may take a moment.

To investigate this we will loop through a series of four temperatures and append each diffusion coefficient to a list.

temperatures = np.array([500, 600, 700, 800])
D = []
analyzers = []

To read these simulations we will use MDAnalysis (however, it is also possible to use data from a VASP simulation). The parser, bootstrap, and diffusion parameters are all defined for all simulations, here we only consider the diffusive regime to begin after 5 ps. Additionally, we include in the p_params a sub_sample_atoms key, this defines the sampling frequency of atoms to be used in the analysis and the sub_sample_traj key, which defined the sampling frequency for the trajectory. This facility can be particularly useful for large simulations where kinisi might encounter issues related to out-of-memory problems.

p_params = {'specie': 'LI',
            'time_step': 5.079648e-4,
            'step_skip': 100,
            'min_dt': 0.001,
            'sub_sample_atoms': 2,
            'sub_sample_traj': 2,
            'progress': False}
u_params = {'progress': False}
d_params = {'progress': False}

File parsing and diffusion determination is then performed in a loop here.

for t in temperatures:
    u = mda.Universe(f'_static/traj_{t}.gro', f'_static/traj_{t}.xtc')
    d = DiffusionAnalyzer.from_universe(u, p_params, uncertainty_params=u_params)
    d.diffusion(10, d_params)
    D.append(d.D)
    analyzers.append(d)

The list of diffusion coefficient objects (which are uravu.distribution.Distribution type objects) and array of temperatures can then be passed to the kinisi.arrhenius.StandardArrhenius class, where we use the bounds keyword argument to give a minimum and maximum value for the activation energy and preexponential factor.

s = StandardArrhenius(temperatures, D, bounds=((0.01, 0.2), (1e-5, 1e-1)))

Having created the object, we can determine the maximum likelihood values for the parameters of activation energy and the preexponential factor.

s.max_likelihood('mini')

After determining the maximum likelihood values, we can use Markov chain Monte Carlo (MCMC) to sample the probability distributions for these.

s.mcmc()

100%|██████████| 1000/1000 [00:58<00:00, 16.98it/s]

We can then visualise the probability distributions for the parameters as histograms.

from corner import corner

corner(s.flatchain, labels=['$E_a$/eV', '$A$/eV'])
plt.show()

It is also possible to plot these probability distributions as Arrhenius relations on the data measured values.

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

plt.errorbar(1000/s.x, s.y.n, s.y.ci(), marker='o', ls='', color='k', zorder=10)
for i, ci in enumerate(credible_intervals):
    plt.fill_between(1000/s.x,
                     *np.percentile(s.distribution, ci, axis=1),
                     alpha=alpha[i],
                     color='#0173B2',
                     lw=0)
plt.yscale('log')
plt.xlabel('$1000T^{-1}$/K$^{-1}$')
plt.ylabel('$D$/cm$^2$s$^{-1}$')
plt.show()