Scripts to calculate ionic diffusivity

Question

I found several AIMD analyzers to estimate mean square displacement and ionic diffusivity such as pymatgen, mo, vaspkit. I tried some tests and found that the same data set might result in different MSD slops by using different analyzers. The ways of calculating MSD are not similar (I do not see vaspkit source code). I wonder which way/script is reliable and its obtained results are accepted?

Answer 1

Short Version

The differences have to do with possible discrepancies in the order of operations when computing MSD's, or with lengths of trajectories in different studies.

For software packages, in addition to the ones listed in the question, there is also kinisi. There is also CHAMPION, which is not published but possible to reverse engineer, and is nicer for complex fluids. See especially Figure 7 in the CHAMPION paper for an example of where differences between approaches come from.

Personally, I would bias towards verifying the underlying physics over trusting any code which does not show it to me.

Long Version: Computing MSD Properly

Nomenclature:

• $\bar{x} = \frac{1}{\Delta} \int_0^\Delta x(\delta) d\delta$ for time-averaging single-ion quantities over multiple starting points in the simulation trajectory of length $\Delta$. We are interested in the displacements over time intervals $t$, $x_1(t,\delta) = r_i(t+\delta) - r_i(\delta)$, and squared displacements, $x_2(t,\delta) = (r_i(t+\delta) - r_i(\delta))^2$. Here, little $\delta$ is the integrand of the time averaging.
• $\langle x \rangle = \frac{1}{N} \sum_1^N x$ for the ensemble averaging over all ions $i$ in the simulation.

There are three operations to keep track of: taking the square of the displacements for a given interval $t$ and starting point $\delta$, taking the time average, and taking the ensemble average.

The proper way of computing the self diffusivity is to first compute single-ion MSD's, then average them over the ensemble of the ions in a simulation, i.e. $\Big\langle \bar{r^2} \Big\rangle$ - bar inside the brackets. I'll call this quantity tMSD. When this quantity is linear with time interval (not simulation length), Fickian diffusion can be assumed and the self diffusion coefficient extracted.

Computing tMSD properly requires custom code because simulation engines pack trajectories by time point ($\delta$), not by ion. Hence all the packages. I have tried to understand what exactly pymatgen and others compute, but in finite time was not able to for layers of abstraction. Hopefully some of the developers can comment in more detail.

Sometimes the computation is done either in reverse order, i.e. take the square, then the ensemble average, then time-average, or skipping the time-averaging. This makes built-in methods such as the computes msd and msd/nongauss in LAMMPS not give reliable answers: they only sample one starting point. I would not recommend using either.

A third order of operations, i.e. first computing the ensemble center-of-mass displacement, then its square, $\langle r \rangle^2$, and then averaging over times, gives the "jump" diffusion coefficient, which yields conductivity rather than diffusivity. Both for real physics and for having effectively one particle, this takes a very long trajectory to converge. I am bringing this up because often from AIMD a unity Haven ratio is assumed, and the self diffusion coefficient used in place of the jump one - which is not always a good assumption.

Long Version : How Long is Long Enough?

Besides possible errors in the order of operations, sometimes the (AI)MD trajectories are simply not long enough, and the results will vary with changing the length of the trajectory. So how does one tell whether a trajectory is "long enough" and a long-time limit is reached? There are a few ways.

• The one used most often for self diffusivity is to verify linearity of tMSD (if done correctly) on a log-log plot vs the time interval $t$. That involves human judgement over when the long-time linear regime starts, i.e. how close to 1 is "good enough".
• The maximum-entropy distribution with the constraints of a mean and a variance is a Gaussian. This means that distributions of tMSD's over longer and longer time intervals will asymptotically tend towards Gaussians. So an indirect test for "long enough" is whether the distributions of tMSD's are Gaussian. This is why the non-Gaussian parameter (NGP) is sometimes used. It has the same drawback: it does not go to zero at finite interval $t$, but approaches it asymptotically.
• Next, there is the ergodicity breaking parameter, $EB(t,\Delta)$, which is the relative variance of the tMSD, very similar to Mandel's Q parameter. This is typically computed for a fixed small-ish $t$ over increasing trajectory lengths $\Delta$. The benefit of using $EB$ is that it can have a well-defined cutoff for a long-time limit, i.e. you can tell when a trajectory is "long enough".
• Fourth, there is the correlation of the diffusion kernel $C_D$, which is another test loosely based on when nonlinear path-dependent effects, which are represented by higher-exponent averages of displacements, become negligible relative to tMSD. Like the $EB$, it has a well-defined long-time limit, and reaches that limit at finite time intervals. The drawback is that it involves a Laplace transform and an inverse Laplace transform, which may need to be tuned numerically.

The python code I wrote for a recent paper is far from great, and is not a package. But as of the time of this writing (July 31, 2022), it is to my knowledge the only one that implements the three metrics above (NGP, $EB$, $C_D$) for verifying that the "infinite-time" limit has been reached in the simulation.

References

• Song et al. Transport dynamics of complex fluids. PNAS 2019, 116 (26), 12733.
• He, Burov, Metzler, Barkai. Random Time-Scale Invariant Diffusion and Transport Coefficients. PRL 2008, 101, 058101.
• Jeon et al. Protein Crowding in Lipid Bilayers Gives Rise to Non-Gaussian Anomalous Lateral Diffusion of Phospholipids and Proteins. PRX 2016, 6, 021006.
• Lanoiselée, Grebenkov. A model of non-Gaussian diffusion in heterogeneous media. J. Phys. A. 2018, 51, 145602.
• Poletayev, Dawson, Islam, Lindenberg. Defect-driven anomalous transport in fast-ion conducting solid electrolytes. Nat. Mater. 2022.