Finding the MSD uncertainties from a VASP file#
Using kinisi to obtain the mean-squared displacement and uncertainty in a VASP Xdatcar type file is straightforward and involves using the DiffusionAnalyzer class.
[1]:
import scipp as sc
import numpy as np
import matplotlib.pyplot as plt
from pymatgen.io.vasp import Xdatcar
from kinisi.analyze import DiffusionAnalyzer
Below, the params dictionary describes some details about our simulation, this is unpacked when passed to the relevant class method. You can see other keyword arguments for a given class method in the relevant documentation. We outline the params included in this example below:
specie: the atomic/ionic species to calculate the MSD fortime_step: the simulation timestep in your molecular dynamics simulation, for this example, this XDATCAR is from a VASP ab-initio molecular dynamics simulation and the simulation timestep is 2 femtoseconds. Note, that we use the fs unit from thescipplibrary.step_skip: the frequency with which the data was written in your simulation trajectory. This XDATCAR contains every 50th molecular dynamics step (NBLOCK=50in VASP parlance).progress: is you want a progress meter to be printed to the screen, this is not necessary for this documentation example.
[2]:
params = {'specie': 'Li',
'time_step': 2.0 * sc.Unit('fs'),
'step_skip': 50 * sc.Unit('dimensionless'),
'progress': False
}
Therefore, for this example, we have a simulation that had a timestep of 2 femtoseconds but we only stored the trajectory information every 50 steps and we want to investigate only the lithium motion.
The next step is to read the data file in to the appropriate analyzer. Below, the XDATCAR file is read into a pymatgen Xdatcar object, called xd. A DiffusionAnalyzer object is then created using the from_xdatcar class method, passing also the params dictionary.
[3]:
xd = Xdatcar('./example_XDATCAR.gz')
diff = DiffusionAnalyzer.from_xdatcar(xd, **params)
The DiffusionAnalyzer will determine the mean-squared displacements and variances (using the variance rescaling approach detailed in the methodology). The diff.msd object is a scipp.DataArray, which can be plotted as shown below.
[4]:
fig, ax = plt.subplots()
ax.errorbar(diff.dt.values, diff.msd.values, np.sqrt(diff.msd.variances))
ax.set_xlabel(f'Time / {diff.dt.unit}')
ax.set_ylabel(f'MSD / {diff.msd.unit}')
ax.set_xlim(0, None)
ax.set_ylim(0, None)
plt.show()
If you want to delve deeper into the data, this can be achieved by investigating the scipp.DataArray directly. This can be interpreted as a HTML object, as shown below.
[5]:
diff.da
[5]:
- time interval: 140
- n_samples(time interval)float64𝟙2.688e+04, 1.344e+04, ..., 193.381, 192.0
Values:
array([26880. , 13440. , 8960. , 6720. , 5376. , 4480. , 3840. , 3360. , 2986.66666667, 2688. , 2443.63636364, 2240. , 2067.69230769, 1920. , 1792. , 1680. , 1581.17647059, 1493.33333333, 1414.73684211, 1344. , 1280. , 1221.81818182, 1168.69565217, 1120. , 1075.2 , 1033.84615385, 995.55555556, 960. , 926.89655172, 896. , 867.09677419, 840. , 814.54545455, 790.58823529, 768. , 746.66666667, 726.48648649, 707.36842105, 689.23076923, 672. , 655.6097561 , 640. , 625.11627907, 610.90909091, 597.33333333, 584.34782609, 571.91489362, 560. , 548.57142857, 537.6 , 527.05882353, 516.92307692, 507.16981132, 497.77777778, 488.72727273, 480. , 471.57894737, 463.44827586, 455.59322034, 448. , 440.6557377 , 433.5483871 , 426.66666667, 420. , 413.53846154, 407.27272727, 401.19402985, 395.29411765, 389.56521739, 384. , 378.5915493 , 373.33333333, 368.21917808, 363.24324324, 358.4 , 353.68421053, 349.09090909, 344.61538462, 340.25316456, 336. , 331.85185185, 327.80487805, 323.85542169, 320. , 316.23529412, 312.55813953, 308.96551724, 305.45454545, 302.02247191, 298.66666667, 295.38461538, 292.17391304, 289.03225806, 285.95744681, 282.94736842, 280. , 277.11340206, 274.28571429, 271.51515152, 268.8 , 266.13861386, 263.52941176, 260.97087379, 258.46153846, 256. , 253.58490566, 251.21495327, 248.88888889, 246.60550459, 244.36363636, 242.16216216, 240. , 237.87610619, 235.78947368, 233.73913043, 231.72413793, 229.74358974, 227.79661017, 225.88235294, 224. , 222.14876033, 220.32786885, 218.53658537, 216.77419355, 215.04 , 213.33333333, 211.65354331, 210. , 208.37209302, 206.76923077, 205.19083969, 203.63636364, 202.10526316, 200.59701493, 199.11111111, 197.64705882, 196.20437956, 194.7826087 , 193.38129496, 192. ]) - time interval(time interval)float64fs100.0, 200.0, ..., 1.390e+04, 1.400e+04
Values:
array([ 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000., 1100., 1200., 1300., 1400., 1500., 1600., 1700., 1800., 1900., 2000., 2100., 2200., 2300., 2400., 2500., 2600., 2700., 2800., 2900., 3000., 3100., 3200., 3300., 3400., 3500., 3600., 3700., 3800., 3900., 4000., 4100., 4200., 4300., 4400., 4500., 4600., 4700., 4800., 4900., 5000., 5100., 5200., 5300., 5400., 5500., 5600., 5700., 5800., 5900., 6000., 6100., 6200., 6300., 6400., 6500., 6600., 6700., 6800., 6900., 7000., 7100., 7200., 7300., 7400., 7500., 7600., 7700., 7800., 7900., 8000., 8100., 8200., 8300., 8400., 8500., 8600., 8700., 8800., 8900., 9000., 9100., 9200., 9300., 9400., 9500., 9600., 9700., 9800., 9900., 10000., 10100., 10200., 10300., 10400., 10500., 10600., 10700., 10800., 10900., 11000., 11100., 11200., 11300., 11400., 11500., 11600., 11700., 11800., 11900., 12000., 12100., 12200., 12300., 12400., 12500., 12600., 12700., 12800., 12900., 13000., 13100., 13200., 13300., 13400., 13500., 13600., 13700., 13800., 13900., 14000.])
- (time interval)float64Å^20.443, 0.707, ..., 11.620, 11.814σ = 0.003, 0.007, ..., 0.902, 0.928
Values:
array([ 0.44302929, 0.70655629, 0.89235476, 1.02097822, 1.12603226, 1.22171929, 1.31557299, 1.41057845, 1.50508755, 1.59926923, 1.69197674, 1.78497179, 1.87238976, 1.95812694, 2.04634983, 2.13427355, 2.2251105 , 2.30746871, 2.3851957 , 2.46455221, 2.54733011, 2.63140368, 2.71515786, 2.80328499, 2.8827326 , 2.96221317, 3.0406781 , 3.12485007, 3.21228671, 3.29702384, 3.38412521, 3.47213267, 3.55604937, 3.64162357, 3.73038097, 3.81936626, 3.90970167, 4.00250713, 4.09432587, 4.17815965, 4.26496592, 4.35031936, 4.43909235, 4.53301713, 4.62667761, 4.71980438, 4.81947935, 4.92633342, 5.02302782, 5.11082792, 5.1927592 , 5.27814133, 5.36412577, 5.45414657, 5.54038982, 5.626678 , 5.70500244, 5.78362404, 5.85729732, 5.93835067, 6.02156549, 6.10097306, 6.18919942, 6.2750682 , 6.35895016, 6.44117195, 6.51954637, 6.59801615, 6.67261818, 6.7428452 , 6.80631239, 6.87210483, 6.95200991, 7.03272087, 7.12144909, 7.20390133, 7.28768648, 7.37263007, 7.44975712, 7.51565205, 7.59379822, 7.67400039, 7.76187187, 7.83711758, 7.91030433, 7.97908311, 8.05480952, 8.13274041, 8.21094614, 8.28441732, 8.35504301, 8.42901659, 8.50298369, 8.56961437, 8.63473134, 8.69366385, 8.74888618, 8.80389912, 8.86750111, 8.9372199 , 9.00375859, 9.07589728, 9.13559701, 9.19586384, 9.23473129, 9.28727536, 9.34760169, 9.41883469, 9.49605429, 9.57657648, 9.67211549, 9.75220801, 9.83102276, 9.92704782, 10.01037201, 10.09124807, 10.15626185, 10.20614993, 10.26117425, 10.33047751, 10.41641562, 10.492222 , 10.54600808, 10.58822323, 10.62851195, 10.6755299 , 10.71922037, 10.77233732, 10.81188068, 10.8736523 , 10.92249725, 10.97909647, 11.04875533, 11.12978573, 11.2050379 , 11.30116001, 11.40843418, 11.52798086, 11.62017274, 11.81370595])
Variances (σ²):
array([8.90150277e-06, 4.63900568e-05, 1.13768803e-04, 2.16162673e-04, 3.53135139e-04, 5.17277896e-04, 7.15683361e-04, 9.50374943e-04, 1.21305914e-03, 1.50009701e-03, 1.85092706e-03, 2.25634059e-03, 2.68901692e-03, 3.17336420e-03, 3.66280313e-03, 4.18565539e-03, 4.76390191e-03, 5.41871625e-03, 6.06143294e-03, 6.78574713e-03, 7.48607571e-03, 8.36842284e-03, 9.34213859e-03, 1.04614844e-02, 1.15256737e-02, 1.26034363e-02, 1.36873181e-02, 1.49494204e-02, 1.62699794e-02, 1.76091563e-02, 1.89552391e-02, 2.05283836e-02, 2.21962498e-02, 2.38804827e-02, 2.57441817e-02, 2.77449334e-02, 2.95460057e-02, 3.14082420e-02, 3.35550237e-02, 3.58050377e-02, 3.82286634e-02, 4.08094523e-02, 4.37185473e-02, 4.65075561e-02, 4.92390271e-02, 5.21393157e-02, 5.53024417e-02, 5.84933260e-02, 6.16680813e-02, 6.51010451e-02, 6.82277730e-02, 7.14265932e-02, 7.52295651e-02, 7.94570080e-02, 8.33941040e-02, 8.72656678e-02, 9.12273513e-02, 9.53928743e-02, 9.97905486e-02, 1.04166775e-01, 1.09020517e-01, 1.13382208e-01, 1.17856569e-01, 1.22505754e-01, 1.26987391e-01, 1.30920965e-01, 1.35144169e-01, 1.39971315e-01, 1.45066849e-01, 1.49710134e-01, 1.54114282e-01, 1.57360116e-01, 1.61702772e-01, 1.66660926e-01, 1.72217467e-01, 1.78043498e-01, 1.83662707e-01, 1.89906426e-01, 1.96385141e-01, 2.02736063e-01, 2.09474056e-01, 2.16081432e-01, 2.23423636e-01, 2.30947414e-01, 2.38458255e-01, 2.46904820e-01, 2.53367311e-01, 2.61291562e-01, 2.70380852e-01, 2.78732917e-01, 2.85797648e-01, 2.93762950e-01, 3.01953014e-01, 3.08914732e-01, 3.15357198e-01, 3.20846442e-01, 3.28050344e-01, 3.35173326e-01, 3.42604832e-01, 3.51409942e-01, 3.57011867e-01, 3.65253185e-01, 3.73396747e-01, 3.80500885e-01, 3.86201400e-01, 3.92079861e-01, 4.00116269e-01, 4.09737609e-01, 4.19204979e-01, 4.28488560e-01, 4.40247714e-01, 4.51290163e-01, 4.61367290e-01, 4.73440769e-01, 4.83825372e-01, 4.93876751e-01, 5.00471355e-01, 5.07809782e-01, 5.18132601e-01, 5.29874755e-01, 5.44973605e-01, 5.60501672e-01, 5.72259712e-01, 5.82141950e-01, 5.94166072e-01, 6.05495879e-01, 6.14831202e-01, 6.20662972e-01, 6.28017009e-01, 6.39987858e-01, 6.54597840e-01, 6.75216284e-01, 6.87267214e-01, 6.99100178e-01, 7.17464270e-01, 7.30560165e-01, 7.56139515e-01, 7.77811256e-01, 8.13531325e-01, 8.60436389e-01])
Note that by giving the appropriate units to the time_step input, the units of the time interval (the duration in which the particle travelled) can be correctly calculated.