ABSTRACT Ab initio molecular dynamics (AIMD) simulation is widely employed in studying diffusion mechanisms and in quantifying diffusional properties of materials. However, AIMD simulations

are often limited to a few hundred atoms and a short, sub-nanosecond physical timescale, which leads to models that include only a limited number of diffusion events. As a result, the

diffusional properties obtained from AIMD simulations are often plagued by poor statistics. In this paper, we re-examine the process to estimate diffusivity and ionic conductivity from the

AIMD simulations and establish the procedure to minimize the fitting errors. In addition, we propose methods for quantifying the statistical variance of the diffusivity and ionic

diffusional properties. Our work provides the foundation for quantifying the statistical confidence levels of diffusion results from AIMD simulations and for correctly employing this

powerful technique. SIMILAR CONTENT BEING VIEWED BY OTHERS COMPUTING INELASTIC NEUTRON SCATTERING SPECTRA FROM MOLECULAR DYNAMICS TRAJECTORIES Article Open access 12 April 2021 CONVERGENCE

AND EQUILIBRIUM IN MOLECULAR DYNAMICS SIMULATIONS Article Open access 07 February 2024 APPLYING BAYESIAN INFERENCE AND DETERMINISTIC ANISOTROPY TO RETRIEVE THE MOLECULAR STRUCTURE ∣Ψ(R)∣2

DISTRIBUTION FROM GAS-PHASE DIFFRACTION EXPERIMENTS Article Open access 13 November 2023 INTRODUCTION As a powerful modeling technique, ab initio molecular dynamics (AIMD) simulation has

been recently applied to a wide range of research topics in chemistry and materials science.1,2,3,4 AIMD simulations are carried out using the potential energy surface and atomistic forces

calculated from ab initio methods, such as density functional theory (DFT), and model the dynamics of atomistic systems with ab initio level of accuracy and chemical versatility. This

accuracy is lacking in classical molecular dynamics (MD) simulations, which are based on interatomic potentials (a.k.a. force fields). In addition, the limited availability of reliable

interatomic potentials for studying new materials greatly limits the applicability and chemical versatility of classical MD simulations. Therefore, AIMD simulation is the method of choice

for studying the dynamics of atoms with complex chemistry changes and simulating materials that cannot be described by available interatomic potentials. The wide applicability of AIMD

simulations has been successfully demonstrated in studies of diffusional properties,5,6,7,8,9 reaction processes,10,11 vibrational frequency,12,13,14 amorphous materials,15,16,17 phase

transition,18 etc. Mo et al.5 pioneered the application of AIMD simulations to the study of ionic diffusion in lithium ionic conductor materials, such as Li10GeP2S12 (LGPS). Since then, AIMD

simulation has been widely adopted as a standard technique for studying fast ionic conductor materials with a wide range of mobile species (e.g., Li+, Na+, Mg2+, O2−).6,7,19,20,21,22,23

Both ab initio and classical MD simulations model the real-time dynamics of atoms with a femtosecond-scale time resolution, which is difficult to directly obtain from experimental

characterizations or from other non-MD computational methods. Due to their ab initio level of accuracy and chemical versatility, AIMD simulations provide more accurate quantification of

diffusion properties than classical MD and have a unique ability to predict new diffusion mechanisms. AIMD simulation technique has achieved great successes in identifying diffusion

mechanism at the atomistic scale and in quantifying the diffusional properties. Major advantages of AIMD simulations (as well as classical MD simulations) are that the diffusion events,

i.e., atom/ion hops, in the materials are directly observed from the dynamics and trajectory of atoms and that specific diffusion mechanism is not presumed as calculation input. In contrast,

the widely used nudged-elastic-band (NEB) calculations only obtain the energy barrier of a specific static migration pathway and require the pathway as input. That NEB input pathway is

often derived from the guess of the researchers. The prior guess of diffusion mechanism may be difficult for complex materials such as super-ionic conductors, where the mobile-ion sublattice

is highly disordered. For example, recent AIMD simulations identified that the diffusion mechanism in lithium super-ionic conductors, such as thio-LISICON9,24 and lithium garnet,9,25,26 is

a concerted migration of multiple ions rather than isolated ion hopping. AIMD simulations played unique roles in uncovering and identifying these concerted migration mechanisms, whereas

isolated ion hopping was often assumed for previous NEB calculations. In addition, AIMD simulations sample the statistical contributions of all diffusional events and obtain diffusional

properties such as diffusivity and ionic conductivity from the collective contributions of all different diffusion modes. In addition, performing AIMD simulations at different temperatures

obtains the Arrhenius relation of diffusivity and ionic conductivity, including the pre-factor and overall activation energy, which is similar to experimental measurements. However, the high

computational expenses of ab initio calculations limit the accessible length scale and timescale of AIMD simulations. The small system size and short physical time scale are major

shortcomings of AIMD simulations when compared to classical MD simulations.27,28 Even on state-of-the-art supercomputing facilities, typical AIMD simulations are limited to a system size of

a few hundred atoms and a total physical time duration of tens to thousands of picoseconds. As a result, a limited number of diffusion events, i.e., ion hops, are observed during an AIMD

simulation. The limited sampling of diffusion events may lead to significant statistical variances in calculated diffusional properties or even, if only a few ion hops are sampled, to

erroneous values. As AIMD simulation has become a widely used technique in studying diffusion mechanisms and in quantifying diffusional properties of materials, the procedures for applying

AIMD simulations to extract diffusional properties need to be systematically established. First, the fitting procedure of AIMD simulation results requires a number of steps and parameters.

Currently, every study chooses their own tested parameters, but the choice of these parameters has not been discussed. As we show in this study, the choice may significantly impact obtained

diffusional properties. Second, most previous studies correlated the error bar and statistical uncertainty of diffusivity to the goodness of linear fitting to the Einstein relation. However,

as we show in this study, the goodness of linear fitting to the Einstein relation does not necessarily capture the true statistical variance of fitted diffusivity. Instead, the statistical

variance and uncertainty of fitted diffusivity should be a direct consequence of the limited observation of diffusion events (i.e., ion hops) during the MD simulation. In the literature,

several studies estimated the errors in the fitted diffusivity from single-particle tracking,29,30,31,32 classical MD simulations33 and kinetic Monte Carlo simulations.34 However, these

studies, which are developed and tested for systems and measurements with a significantly larger number of diffusion events than typical AIMD simulations, are not ideal for analyzing errors

from AIMD simulations, where the total number of diffusion events are small. Third, a significant number of diffusion events should be captured to quantify the diffusion with statistical

validity. However, in many studies, the AIMD simulations only capture and sample ion diffusion from a few ion hops with a total mean-squared displacement (MSD) of a few Angstrom.2 The

diffusion properties from such poor sampling of ion hops may have low statistical significance. The limitation of AIMD simulations should be quantified and charted out. In summary, the

following three problems of using AIMD simulations to study diffusional properties need to be addressed: (1) How shall one properly extract diffusional properties from the dynamics of atoms

with minimal error? (2) How shall one quantify the statistical variances of diffusional properties extracted from AIMD simulations? (3) What are the accessible ranges of materials properties

and physical conditions for typical AIMD simulations in studying diffusion? The aim of this paper is to resolve the aforementioned problems and to establish a systematic procedure for

quantifying diffusional properties and their statistical variances from AIMD simulations. In the following sections, we examine and establish the proper procedure to extract diffusivity _D_

from the atomic trajectory of AIMD simulations. We illustrate the origin and the functional dependence of the statistical variance of the diffusivity _D_. We estimate the proper temperature

range for AIMD simulations to obtain _D_ with reasonable accuracy. We estimate the errors of activation energy and ionic conductivity from the Arrhenius relation for typical super-ionic

conductor materials. RESULTS QUANTIFYING DIFFUSIVITY, IONIC CONDUCTIVITY, AND ACTIVATION ENERGY FROM AIMD SIMULATIONS The diffusional properties are calculated from the trajectory of ions

(or atoms) \({\boldsymbol{r}}_i(t)\) from AIMD simulations. The displacement \(\Delta {\boldsymbol{r}}_i\) of ion _i_ from time _t_1 to _t_2 can be calculated as $$\Delta

{\boldsymbol{r}}_i\left( {\Delta t} \right) = {\boldsymbol{r}}_i\left( {t_2} \right) - {\boldsymbol{r}}_i\left( {t_1} \right),\,{\mathrm{where}}\,\Delta t = t_2 - t_1.$$ (1) The total

squared displacement sums up the squared displacement of all mobile ions, \(\mathop {\sum}\nolimits_{i = 1}^N {\left( {\left| {\Delta r_i\left( {\Delta t} \right)} \right|^2} \right)}\), and

describes the movement of all _N_ mobile ions over a time interval Δ_t_. During AIMD simulations over a total time duration _t_tot, there are many (_N_Δ_t_) time intervals with the same

duration Δ_t_ (Δ_t_ _<_ _t_tot) but with different starting time _t_. Since the displacement of ions over Δ_t_ reflects the mobility of ions, the total mean squared displacement (TMSD)

for all diffusional ions over time interval Δ_t_ is calculated as: $${\mathrm {TMSD}}\left( {\Delta t} \right) = \mathop {\sum}\nolimits_{i = 1}^N {\left\langle {\left|

{{\boldsymbol{r}}_i\left( {\Delta t} \right) - {\boldsymbol{r}}_i(0)} \right|^2} \right\rangle } = \mathop {\sum}\nolimits_{i = 1}^N {\frac{1}{{N_{\Delta t}}}} \mathop {\sum}\nolimits_{t =

0}^{t_{tot} - \Delta t} {\left| {{\boldsymbol{r}}_i\left( {t + \Delta t} \right) - {\boldsymbol{r}}_i(t)} \right|^2},$$ (2) by averaging over a total of _N_Δ_t_ time intervals with the same

duration Δ_t_. This averaging over different _N_Δ_t_ time intervals provides essential ensemble sampling to obtain more accurate diffusional properties. To estimate the diffusivity of the

mobile-ion species, the MSD over time interval Δ_t_ is calculated as the TMSD per mobile ion: $${\mathrm {MSD}}\left( {\Delta t} \right) = \frac{1}{N}{\mathrm {TMSD}}\left( {\Delta t}

\right).$$ (3) _N_ is the number of ions that are assumed to be the mobile carriers contributing to diffusion. In general, the dependence of MSD over time interval Δ_t_ follows a linear

relationship if a large amount of diffusional displacement is captured during the MD simulation. The diffusivity of these ions is calculated as the slope of the MSD over time interval Δ_t_

according to the Einstein relation: $$D = \frac{{{\mathrm {MSD}}\left( {\Delta t} \right)}}{{2d\Delta t}} + D_{{\mathrm{offset}}},$$ (4) where _d_ = 3 is the dimension of the system, and the

offset _D_offset of this linear dependence is discussed in later sections. This calculated diffusivity _D_ is the tracer diffusivity of the mobile-ion species and is an intrinsic property

of the material at the given condition. From the diffusivity _D_, the ionic conductivity is calculated based on the Nernst–Einstein relation: $$\sigma = \frac{{Nq^2}}{{VkT}}D,$$ (5) where

_V_ is the total volume of the model system, _q_ is the charge of the mobile-ion species, _T_ is the temperature, and _k_ is the Boltzmann constant. By combining Eqs. (3–5), the ionic

conductivity is directly determined by TMSD(Δ_t_)/Δ_t_: $$\sigma = \frac{{q^2}}{{VkT}}\frac{{{\mathrm {TMSD}}\left( {\Delta t} \right)}}{{2d\Delta t}}.$$ (6) The ionic conductivity and TMSD

calculated in Eq. (6) are independent of the specific choice of the diffusion carrier, such as vacancy or interstitial. In comparison, the diffusivity _D_ can be calculated for specific

carriers, e.g., Li+ or vacancy, by normalizing the specific number of carriers in Eq. (3). The choice and counting of mobile carriers do not affect the calculated ionic conductivity from Eq.

(6). Therefore, the use of conductivity and TMSD is more straightforward in describing the overall diffusion that occurred during the AIMD simulations. In addition, the Nernst–Einstein

relation (Eq. (5)) assumes dilute, non-interacting mobile ions in the materials systems. However, recent AIMD simulation studies have shown the strong correlation of the ionic diffusion in

fast ionic conductors.9,24,25,26,35,36,37 The Haven ratio is often used to describe the correlation factor, and can be quantified by the ratio of the jump diffusion coefficient _D_J against

the tracer diffusivity _D_ calculated from Eq. (4).38,39,40 The jump diffusion coefficient _D_J, which describes the collective migration of all ions, is estimated from the MSD of the mass

center of all mobile ions. However, our analysis found that tracking the displacement of the single center point exhibits higher statistical variance. In this study, the scheme for

estimating statistical variance is developed for tracer diffusivity _D_ in Eqs. (5, 6), but a similar scheme can be applied and developed for analyzing the variance of _D_J. As in the

experiments,41,42,43,44 AIMD simulations can be performed at multiple temperatures to obtain the Arrhenius relation of _D_ as a function of _T_: $$D = D_0{\mathrm{exp}}\left( { -

\frac{{E_{\mathrm a}}}{{kT}}} \right).$$ (7) The activation energy _E_a of ion diffusion can be obtained through fitting the data of log _D_ vs. 1/_T_ to the Arrhenius relationship (Eq.

(7)). The fitted Arrhenius relationship can be used to extrapolate the diffusivity _D_ and conductivity _σ_ to other temperatures. It should be noted that, by extrapolating this Arrhenius

relation to other temperatures, the identical diffusion mechanism is assumed at those extrapolated temperatures. REGIONS OF MSD–Δ_T_ DEPENDENCE Figure 1 shows a typical MSD–Δ_t_ curve from

AIMD simulations of the Li-ion superionic conductor Li1.33Ti1.67Al0.33(PO4)3 (LATP). The linear MSD–Δ_t_ dependence as described in the Einstein relationship (Eq. (4)) only holds within a

certain range of time intervals Δ_t_, and a notable fraction of this dependence is not linear. The MSD–Δ_t_ curve at short time interval Δ_t_ _<_ 0.1 ps follows \({\mathrm {MSD}} \propto

\Delta t^{1.42}\), which is consistent with the local harmonic vibration motion model as shown in the Supplementary Note 1 and Supplementary Figure 1. This portion of the MSD–Δ_t_ curve is

named the ballistic region, corresponding to the ballistic and vibrational motion of Li ions around their local equilibrium sites rather than Li-ion hopping to new sites.33,45 More details

of this ballistic region are discussed in the next section. In Δ_t_ range of 10–100 ps, the MSD–Δ_t_ curve follows \({\mathrm {MSD}} \propto \Delta t\), the linear Einstein relationship for

diffusional displacement. However, when Δ_t_ reaches a large fraction of _t_tot (>100 ps in the case of Fig. 1), the MSD_–_Δ_t_ curve shows notable variance and deviation from the linear

relationship. This deviation is a result of poor statistics for large Δ_t_ values, as values that are comparable to _t_tot result in averaging a smaller number of time intervals _N_Δ_t_ in

Eq. (2) than would be averaged for small Δ_t_ values that are a small fraction of _t_tot. The effect of this deviation is quantified in the next section. Given that the linear MSD–Δ_t_

dependence only holds at certain Δ_t_ range, one should not utilize the entire MSD–Δ_t_ curve to fit the Einstein relation (Eq. (4)) and to deduce the diffusivity. The ballistic region at

short times intervals and the poor statistical region at large time intervals should be excluded from fitting diffusivity. The linear fitting of MSD–Δ_t_ curve should be performed in a range

[Δ_t_ _low_ , Δ_t_ _up_ ], where Δ_t_ _low_ and Δ_t_ _up_ are the lower and upper bound, respectively. In last sections, we establish the procedures to determine the bounds Δ_t_ _low_ and

Δ_t_ _up_ for the linear fitting of MSD–Δ_t_ in order to minimize the errors caused by improper fitting to the Einstein relation. LOWER BOUND OF LINEAR DIFFUSION REGION In Fig. 2, the local

derivative dMSD/dΔ_t_ calculated using a finite difference method shows the transition from the ballistic region to the linear region and is used to determine the lower fitting bound

Δ_t_low. The derivative dMSD/dΔ_t_ has a large value at small Δ_t_ values (_<_0.2 ps), and decreases significantly as Δ_t_ increases. The MSD–Δ_t_ curve becomes linear at Δ_t_ ≳ ~ 1 ps,

and the value of dMSD/dΔ_t_ reaches a plateau value of ~3 Å2/ps. In the case of Fig. 2, the cut-off of the ballistic region (black dotted line in Fig. 2) is at an MSD value of ~5 Å2. An MSD

of 5 Å2 is approximately 0.5_a_2, where _a_ is the distance between two neighboring Li sites and is 3.2 Å in LATP, which is typical of Li ionic conductors. The ballistic region has the MSD

cut-off at a fraction of _a_2 because the local vibration of Li ions on their equilibrium sites is confined within the potential well between two nearby sites. This local vibrational

displacement at small Δ_t_ _<_ Δ_t_low does not represent the ionic diffusion from site to site, and should be excluded from the linear fitting for _D_. One may identify the specific

range of the ballistic region and the precise lower fitting bound Δ_t_low for each individual AIMD simulation of different materials systems at different temperatures using the same

procedure shown in Figs. 1 and 2. The cut-off values of Δ_t_ (i.e., Δ_t_low) and MSD for the ballistic region are dependent on the temperatures and the materials, which yield different

potential energy surfaces near the equilibrium sites. In general, the displacement of local vibration is within a fraction of site distance _a_, so the ballistic region corresponds to the

MSD ranging from 0 to a fraction of _a_2 (Supplementary Note 1 and Supplementary Figure 1). Therefore, we propose to quantify the cut-off values of Δ_t_, i.e., Δ_t_low, by defining the

cut-off in MSD using the value of a fraction of _a_2 (e.g., 0.5 _a_2). As a safer measure for lower fitting bound (i.e., Δ_t_low), the entire region with MSD less than this cut-off based on

_a_2 may be excluded. The linear fitting to Einstein relation should only be performed on the linear region corresponding to diffusional displacement. Otherwise, a part of ballistic motion

would be included into the fitting for _D_, leading to an over-estimation of _D_. This over-estimation of _D_ would be pronounced for the MSD–Δ_t_ curves with small values of the maximum MSD

(~ a few Å2), due to a significant fraction of ballistic motion mixed into the MSD–Δ_t_ curve. Therefore, AIMD simulation should be long enough to have the MSD per mobile ion larger than a

few _a_2, so that the ballistic region (and Δ_t_low) can be distinguished and separated from the linear fitting for _D_. UPPER BOUND OF LINEAR DIFFUSION REGION The upper fitting bound Δ_t_up

is determined by the transition from the linear diffusion region to the region at large Δ_t_ with large variance and deviation. Ten different MSD–Δ_t_ curves from AIMD simulations of the

same LATP structure model at 1200 K over 50 ps (Fig. 3a) were obtained by dividing a total AIMD simulation over 500 ps into ten non-overlapping parts. The significant deviations from the

linear dependence of these ten MSD–Δ_t_ curves are typically observed at large values of Δ_t_, i.e., >25 ps in Fig. 3 or ≳50% of _t_tot. At large Δ_t_, a smaller number of time intervals

_N_Δ_t_ is averaged in Eq. (2) and many of these Δ_t_ intervals overlap other intervals that contain physically identical trajectories of ions, thus leading to larger deviation from the

linear MSD–Δ_t_ relation. Therefore, the linear fitting for _D_ should be performed below the upper bound Δ_t_up of the linear diffusion region. To determine the upper fitting bound Δ_t_up,

we fitted these MSD–Δ_t_ curves to the Einstein relation (Eq. (4)) with different Δ_t_up values ranging from 10% to 100% of _t_tot. For the fitting of each curve, the value of _R_2 in the

linear regression was calculated to evaluate the goodness of fitting. The average and standard deviation of _R_2 values over these ten curves are shown in Fig. 3b. At Δ_t_up ≤ 0.3_t_tot, all

values of _R_2 are very close to 1, showing good linearity of MSD–Δ_t_ curves at small Δ_t_. At Δ_t_up > 0.7_t_tot, _R_2 values decrease significantly from 1, indicating poor linearity

of MSD–Δ_t_ curves at large Δ_t_, and the standard deviation of _R_2 also increases. Therefore, one should only fit the linear region of MSD–Δ_t_ below an upper fitting bound of Δ_t_up < 

0.7_t_tot. By performing the same test on other fast ion conductor materials, such as LGPS and Li7La3Zr2O12 (LLZO) at a few different temperatures, we found this Δ_t_up < 0.7_t_tot to be

generally applicable (Supplementary Figure 3). As the optimal Δ_t_up values may depend on the material (the mobile ions) and temperatures, for unique systems one may use this scheme to

determine the specific Δ_t_up values. Within properly determined lower and upper fitting bounds, the goodness of fitting to the MSD–Δ_t_ curve would always be good. Therefore, the goodness

of fitting itself does not reflect the statistical variance in the fitted _D_, the slope of the MSD–Δ_t_ curve, and is different among different AIMD simulations for the same materials model

(Fig. 3a). The changes in the slopes of MSD–Δ_t_ curves below Δ_t_up reflect the statistical variances in the fitted diffusivity _D_ from different runs of AIMD simulations, which are

quantified and analyzed in the next section. STATISTICAL VARIANCE OF DIFFUSIVITY AND CONDUCTIVITY In Fig. 3a, the diffusivities, i.e., the slopes of different MSD–Δ_t_ curves, exhibit

significant variance among different AIMD simulations of the exact same material model. The variance among the fitted diffusivities _D_ are a result of the stochastic nature of the diffusion

process, which causes different numbers of ion hops during AIMD simulations with different initial conditions. In this section, we quantify the statistical variance of the fitted

diffusivity _D_ from AIMD simulations. The statistics were calculated based on a set of MD simulations that were created by dividing a long MD simulation into several non-overlapping shorter

MD simulations, each of which was treated as an individual MD simulation. Following our established fitting procedure, the value of _D_ was extracted from each MD simulation. The standard

deviation of _D_, _s_ _D_ , was calculated from the set of fitted _D_ values (Supplementary Figure 2). The relative standard deviation (RSD) of _D_ was calculated as _s_ _D_ _/D_true, where

the true value of the diffusivity _D_true is calculated from the longest available MD simulation. Since more ion hops improve the sampling of the diffusional property, we found that the RSD

of _D_ decreases with the total effective ion hops _N_eff as $$\frac{{S_D}}{{D_{{\mathrm{true}}}}} = \frac{A}{{\sqrt {N_{{\mathrm{eff}}}} }} + B.$$ (8) _N_eff is calculated as

$$N_{{\mathrm{eff}}} = \frac{{\max _{\Delta t}\left[ {TMSD\left( {\Delta t} \right)} \right]}}{{a^2}},$$ (9) where \({\max_{\Delta t}\left[ {{\mathrm{TMSD}}\left( {\Delta t} \right)}

\right]},\) is the maximum value of TMSD over the entire MSD–∆_t_ curve. _N_eff can be considered as the effective number of ion hops that contributed to the TMSD of all mobile ions in the

entire duration of the MD simulation. The value of _a_ is 2.4, 3.2, 2.8, 3.4 Å for LLZO, LATP, LGPS, RbAg4I5, respectively, as averaged distance between neighboring mobile-ion sites. For all

materials in Fig. 4, the values of _A_ and _B_ are fitted as 3.43 and 0.04, respectively. As shown in Fig. 4, Eq. (8) and its parameters are general for ionic conductor materials and also

hold for a classical MD simulation of LLZO over much longer time duration with a large number of ion hops (Fig. 4b). Equation (8) can be used to estimate the statistical variance of _D_

calculated from any given AIMD simulation, where the number of total effective ion hops _N_eff can be first estimated from the maximum TMSD of all mobile ions using Eq. (9). For example, in

a material with site distance _a_ = 3 Å (assumed for all following estimations), an AIMD simulation that reaches a maximum TMSD of 1800 Å2, corresponding to _N_eff = ~ 200, results in an RSD

of _D_ of ~28%. In addition, since the diffusivity _D_ and ionic conductivity _σ_ are linearly dependent (Eq. (5)), the RSD of _σ_ is equivalent to the RSD of _D_. Thus, Eq. (8) is also

valid for estimating the _RSD_ of ionic conductivity _σ_ from AIMD simulations. This result suggests that in order to obtain more accurate diffusivity from AIMD simulations, one should run

longer MD simulations that sample more diffusion events, i.e., have larger _N_eff. For AIMD simulations with an estimated _N_eff of 50, 100, and 150, corresponding to max (TMSD) of 450, 900,

and 1350 Å2, respectively, RSD of _D_ is 52%, 38%, and 32%, respectively. These levels of RSD are reasonable, but still allow for a notable statistical error of the fitted _D_. To obtain

_D_ with <20% RSD, the AIMD simulation should observe more than max (TMSD) of ~4150 Å2 (_N_eff = ~ 460). The RSD is reduced to ~ 10% for _N_eff = 3200 and max (_TMSD_) = ~ 30,000 Å2. It

should be noted that _N_eff and TMSD are from all mobile ions, while only MSD per ion is often presented in the literature. A max(TMSD) of a few thousand Å2 is quite significant for typical

AIMD simulations within 1 ns and can only be reached in fast ion conductors at relatively high temperatures. In addition to running longer MD simulations, running MD simulations on larger

systems with more mobile ions, for example as is done in classical MD simulations, can achieve more effective ion hops _N_eff and hence less statistical variance in _D_ and _σ_. ACCESSIBLE

RANGE OF DIFFUSIVITY, ACTIVATION ENERGY, AND TEMPERATURE Given the short physical time duration of 100 ps to 1 ns in typical AIMD simulations, AIMD simulations of materials with large

activation energy _E_a or at low temperature _T_ may not observe enough number of ion hops. In order to achieve a reasonable accuracy of the fitted _D_, AIMD simulations may only be

applicable to materials with relatively high ionic conductivity. Given the high computational costs of AIMD simulations, it is often desired to pre-estimate the range of temperatures that

can be performed for a given material (with an estimated _E_a). The accessible ranges of _T_ at a given _E_a can be estimated as follows. The estimation here assumes an AIMD simulation over

a total physical time duration of 1 ns and a supercell with a volume _V_ of 1000 Å3 (i.e., 10 Å × 10 Å × 10 Å) with _N_ = 20 mobile ions and a site distance _a_ = 3 Å. In order to achieve

<50% RSD of _D_ and _σ_, _N_eff should be >55 according to Eq. (8). For this RSD limit, the ionic conductivity should be >0.025 S/cm at 600 K to achieve _N_eff > 55 over 1 ns

(Eqs. (4), (5)). This ionic conductivity corresponds to a Li+ diffusivity of ~4.2 × 10–7 cm2/s (Eq. (5)) for the assumed supercell model. For the materials supercell with different size _V_

or different number of mobile carriers _N_, the accessible range of diffusivity can be estimated in the same way using Eqs. (5–8). In general, a minimum diffusivity of ~10−7 cm2/s is

accessible by the AIMD simulations with typical size (~1 nm3) and time duration (~1 ns), in order to have a reasonable number of effective ion hops (_N_eff > 50). Figure 5 shows a plot of

_D_ as a function of _E_a and _T_ estimated from the transition state theory $$D = a^2v^ \ast \exp \left( { - \frac{{E_a}}{{kT}}} \right),$$ (10) where \(v^ \ast\) is the attempting jump

frequency and is chosen to be 1012 Hz, _a_ is the neighboring-site distance, and _k_ is the Boltzmann constant. The geometric factor and correlation factor, which also affect the value of

_D_, were neglected in this back-of-the-envelope estimation. This plot can be used to estimate the range of _E_a and _T_ accessible to typical AIMD simulations. In the same example model

above, a minimum diffusivity of ~5 × 10–7 cm2/s is needed to achieve an RSD of ~50% within 1 ns. The highest _E_a and lowest temperatures to have the desired _D_ value can be read from this

plot (Fig. 5). For materials with an _E_a of 0.6 eV, as in some cathode materials, AIMD simulations below 900 K may not observe enough ion hops and the fitted _D_ would have a large error.

For super ionic conductors with an _E_a of 0.20 eV or lower, the temperature accessible by AIMD simulations may go as low as 300 K. In addition, this plot can be used to determine the

appropriate temperatures for AIMD simulations to obtain _D_ with reasonable accuracy. For example, in order to achieve a desired RSD of ~20%, a maximum TMSD of 4150 Å2 over 1 ns is needed

(Eqs. (7–8)), and this maximum TMSD corresponds to a minimum diffusivity of ~3 × 10–6 cm2/s (Eq. (5)) in the same example material model. To achieve that diffusivity, AIMD simulations of a

material with _E_a of 0.2 eV need to be at above 450 K (Fig. 5). An ionic conductor with _E_a of 0.3 eV should have AIMD simulations above 700 K to achieve similar accuracy (RSD = ~20%). An

AIMD simulation at >1150 K is needed for an ionic conductor with _E_a of 0.5 eV. Therefore, for materials with slow diffusion, long AIMD simulations at high temperatures are essential to

obtain _D_ with low statistical uncertainty. In general, AIMD simulations containing a significant number of ion hops are crucial for achieving small error bounds and a high confidence level

in the fitted _D_. ESTIMATING ERRORS OF DIFFUSIONAL PROPERTIES FROM ARRHENIUS RELATION Given the statistical uncertainty of fitted _D_, it is crucial to include the statistical uncertainty

of every _D_ data point into the fitting of the Arrhenius relation. In particular, the data points from the simulations at lower temperatures or over shorter durations may have a

significantly smaller number of ion hops and thus a higher uncertainty, and this uncertainty should be taken into account during the fitting. The data points in Fig. 6 indeed show

significantly larger variance at lower temperatures. Therefore, these lower-temperature data points should have less weight in the fitting. As a common practice in linear regression, the

inverse square of the standard deviation of log(_D_) is factored into the fitting as the weight of each data point to account for the statistical uncertainty.46 Given that the fitting of the

Arrhenius relationship is usually performed as a linear fitting of log(_D_) and 1/_T_, the derived standard deviation of log(_D_) shall be used as the weight on each point (as shown as

error bars in Fig. 6) for the linear fitting (Supplementary Note 2). Therefore, the error bounds and statistical intervals of _E_a and _D_0 in the Arrhenius relation (Eq. (7)) can be

obtained using standard error analysis for linear regression, and so are also the errors for diffusivity and conductivity when extrapolated to other temperatures. Table 1 shows the error

bounds of the activation energy _E_a and ionic conductivity _σ_ extrapolated to 300 K for LATP, LGPS and LLZO Li-ion conductors. _E_a and extrapolated _σ_ at 300 K agree well with

experimental values, showing the capability of AIMD simulations for quantifying diffusional properties.41,42,47,48 For all three materials, the _E_a has a standard deviation of ~0.02 eV, and

the extrapolated _σ_ at 300 K has the error bound within 1–2 orders of magnitude. It should be noted that extrapolated conductivity (i.e., room temperature conductivity) assumes that there

is no change in the slope of the Arrhenius relation due to a phase change. While excellent agreement of AIMD simulations and experiments is widely reported, one should note that the

statistical uncertainty of extrapolated _σ_ at 300 K can be as large as an order of magnitude for typical AIMD simulations. The errors in _E_a and _σ_ may be larger for AIMD simulations over

shorter duration or on slower ion conductors, which sample fewer diffusion events. As shown in Supplementary Note 3 and Supplementary Figure 4, using short AIMD simulations, especially at

low temperatures, would lead to significant overestimation of diffusivity and, in particular, the RT conductivity. In addition, the errors of _E_a and extrapolated _σ_ also depend on the

number of data points and the linearity of the relationship according to the error analyses of linear regression.46 To minimize the error of _E_a and extrapolated _σ_, more data points are

helpful. For example, the error of _E_a and _σ_ may increase significantly if some data points are omitted for the LGPS data in Fig. 6, leading to more deviations from the experimental

values (Supplementary Figure 5 and Supplementary Table 1). Therefore, having more data points of _D_ (from sufficiently long AIMD simulations) would lead to smaller error and less

statistical uncertainty for fitted diffusional properties. DISCUSSION In this work, we systematically test and establish the procedures to obtain the diffusional properties and their

statistical variances from AIMD simulations. Only the linear region of the MSD–Δ_t_ curve within the lower and upper fitting bound Δ_t_low and Δ_t_up, which corresponds to diffusional

displacements, should be used to fit the diffusivity. The cut-off of the ballistic region, Δ_t_low, is at the MSD per mobile ion of less than a fraction of _a_2 (_a_ is the distance between

neighboring mobile-ion sites), since the local ballistic vibration is limited within the potential well of the equilibrium site. The AIMD simulation should be long enough to have the MSD per

mobile ion larger than a few times _a_2 so that the distinction between the linear diffusion region and the ballistic region can be clearly observed. Otherwise, a part of ballistic motion

would be mixed into the fitting, leading to the over-estimation of _D_. Moreover, the upper fitting bound Δ_t_up should be set as <30–70% of _t_tot to exclude the region with poor

linearity at large Δ_t_. An adequately long AIMD simulation with properly determined lower and upper fitting bounds is crucial to have a linear MSD–Δ_t_ curve for fitting _D_. If the proper

procedure is followed to determine the linear diffusion region of the MSD–Δ_t_ curve, the goodness of fitting to the MSD–Δ_t_ curve would be nearly perfect and would not reflect the

statistical variance of the fitted diffusivity. The statistical variance of the fitted diffusivity and conductivity is instead a direct result of the total number of ion hops sampled during

the MD simulations. The statistical variances of _D_ from fitting Einstein relation should be quantified by the total displacement TMSD using the empirically determined relation in Eq. (8).

On the basis of these results, AIMD simulations should be adequately long to collect a sufficient number of ion hopping events. As shown in typical alkali-ion conductor materials, a simple

rule of thumb is that a maximum TMSD of a few hundred Å2 is necessary to obtain a reasonable estimation of _D_ (with ~50% RSD), and a maximum TMSD of a few thousand Å2 is necessary to obtain

a more accurate _D_ (~20–30% RSD). Given the high uncertainty of _D_ from typical AIMD simulations, the error bounds of diffusivity and conductivity should be reported. Given the need to

collect a large number of ion hops, AIMD simulation is more suitable for studying fast-ion-conducting materials. Likewise, it is necessary to perform AIMD simulations at high temperature for

most materials. A minimum diffusivity of ~10−7 cm2/s is necessary for a material to be accessible by AIMD simulations and a diffusivity of ~10−5 cm2/s may allow for the calculation of a

more accurate _D_. On the basis of these estimated bounds of minimum diffusivity, the approximate accessible temperature ranges of AIMD simulations may be estimated using the plot in Fig. 5

or using a back-of-the-envelope estimation based on the Arrhenius relation. For materials with a high activation energy, running AIMD simulations at sufficiently high temperatures is crucial

for simulating an adequate number of ion hops. In addition, it is also crucial to include the statistical error bounds of each diffusivity data point when fitting to the Arrhenius relation

(Eq. (7)). The error bounds of _D_ should be accounted for as a weight into the linear fitting of log(_D_) over 1/_T_. The error of _E_a and extrapolated _σ_ can be established through

standard regression error analysis. Such error analyses indicate that the typical statistical error bounds of extrapolated _σ_ at 300 K may be as large as an order of magnitude, even for

fast ionic conductors with many data points at different temperatures and with long AIMD simulations. One should consider the statistical variances when interpreting the diffusional

properties calculated from AIMD simulations. While the results of this study are largely based on AIMD simulations, our analyses, schemes, and conclusions are applicable for diffusional

studies using classical MD simulations. While AIMD simulations provide ab initio potential energy surfaces, leading to more accurate diffusion properties, classical MD simulations can be

performed on significantly larger model systems with more mobile ions and over longer time scales, so more diffusional events can be observed. When the TMSD and total effective ion hops

_N_eff from classical MD simulations are large, the statistical variances of diffusional properties are expected to be small. However, in classical MD simulations in which TMSD or _N_eff is

small, which is common for systems with slow diffusion or at low temperatures, our scheme and analyses should be carried out to properly quantify diffusional properties and their statistical

variances. In summary, the major conclusions for calculating diffusional properties from AIMD simulations are as follows. These conclusions also serve as specific guidelines for future

practices. * 1. In the calculation of diffusivity, fitting to the Einstein relation (Eq. (4)) should be limited to the linear diffusion region of the MSD–Δ_t_ curve. This linear diffusion

region excludes the ballistic region at MSD ≲ _a_2 and the poor linearity region at large Δ_t_ > Δ_t_up (Δ_t_up < 0.7_t_tot). * 2. The AIMD simulations should be long enough and should

only be performed on materials with reasonably fast diffusion so that a large number of diffusion events can be captured to obtain accurate diffusivity. In addition, the MSD per mobile ion

should be larger than a few times _a_2 to distinguish the ballistic region (i.e., Δ_t_low). * 3. The statistical variance of the diffusivity should be derived from the total diffusional

displacements (such as maximum TMSD using Eqs. (8–9)) rather than the goodness of fitting to the Einstein relation. The statistical variance of fitted diffusivity values from AIMD

simulations should be reported. * 4. The statistical uncertainties of each data point of diffusivity should be included when fitting to the Arrhenius relation. The error bounds of activation

energy and extrapolated ionic conductivity are non-negligible and should be evaluated using standard regression analysis. Our study establishes the proper calculation procedures and

statistical error analyses for the correct application of AIMD simulations in estimating diffusional properties. The obtained knowledge and established procedures are the basis for

quantifying diffusional properties, drawing proper conclusions from the AIMD simulation results, and further developing AIMD simulations. METHODS In this study, all DFT calculations were

performed using the Vienna ab initio Simulation Package49 within the projector augmented-wave approach with the Perdew−Burke−Ernzerhof generalized-gradient approximation.50 The initial

structures for AIMD simulations are statically relaxed in DFT using the standard parameters from Materials Project.51,52,53 In AIMD simulations, the DFT-based force evaluations were

non-spin-polarized with a single _Γ_-center _k_-point grid. The AIMD simulations performed in this study were based on the Born–Oppenheimer approximation with a time step of 2 fs. Classical

MD simulations of garnet were performed using Large-scale Atomic/Molecular Massively Parallel Simulator54 in a supercell model of 8 formula units of Li7La3Zr2O12 and the interatomic

potential from ref. 55. All AIMD and classical MD simulations used a time step of 2 fs, an _NVT_ ensemble using a Nose–Hoover thermostat, and a velocity–Verlet algorithm for integrating the

equation of motion. The crystal structures of LGPS, cubic phase LLZO, and RbAg4I5 investigated were obtained from Inorganic Crystal Structure Database (ICSD)56 and Materials

Project_._51,52,53 We used the Li super-ionic conductor LATP based on the NASICON structure to establish the fitting scheme, and then tested the same scheme in other Li super-ionic

conductors. The structure of LATP was derived from the LiTi2(PO4)3 structure by partially substituting Ti with Al and by inserting extra Li atoms. The structures with disordered site

occupancy were ordered using the same method in the previous studies.5,7 DATA AVAILABILITY The computation data to support the findings of this study is available from the corresponding

author on reasonable request. The code to analyze AIMD simulation is available at: https://github.com/mogroupumd/aimd REFERENCES * Carloni, P., Rothlisberger, U. & Parrinello, M. The

role and perspective of ab initio molecular dynamics in the study of biological systems. _Acc. Chem. Res._ 35, 455–464 (2002). Article  Google Scholar  * Iftimie, R., Minary, P. &

Tuckerman, M. E. Ab initio molecular dynamics: concepts, recent developments, and future trends. _Proc. Natl. Acad. Sci. USA_ 102, 6654–6659 (2005). Article  Google Scholar  * Kirchner, B.,

di Dio, P. J. & Hutter, J. _Real-World Predictions from Ab Initio Molecular Dynamics Simulations_. Vol. 307 (Springer, Berlin, Heidelberg, 2011). * Hassanali, A. A., Cuny, J., Verdolino,

V. & Parrinello, M. Aqueous solutions: state of the art in ab initio molecular dynamics. _Philos. Trans. R. Soc. A_ 372, 20120482 (2014). Article  Google Scholar  * Mo, Y., Ong, S. P.

& Ceder, G. First principles study of the Li10GeP2S12 lithium super ionic conductor material. _Chem. Mater._ 24, 15–17 (2012). Article  Google Scholar  * Mo, Y., Ong, S. P. & Ceder,

G. Insights into diffusion mechanisms in P2 layered oxide materials by first-principles calculations. _Chem. Mater._ 26, 5208–5214 (2014). Article  Google Scholar  * He, X. & Mo, Y.

Accelerated materials design of Na0.5Bi0.5TiO3 oxygen ionic conductors based on first principles calculations. _Phys. Chem. Chem. Phys._ 17, 18035–18044 (2015). Article  Google Scholar  *

Deng, Z., Mo, Y. & Ong, S. P. Computational studies of solid-state alkali conduction in rechargeable alkali-ion batteries. _NPG Asia Mater._ 8, e254 (2016). Article  Google Scholar  *

He, X., Zhu, Y. & Mo, Y. Origin of fast ion diffusion in super-ionic conductors. _Nat. Commun._ 8, 15893 (2017). Article  Google Scholar  * He, K. et al. Sodiation kinetics of metal

oxide conversion electrodes: a comparative study with lithiation. _Nano Lett._ 15, 5755–5763 (2015). Article  Google Scholar  * Mosconi, E., Azpiroz, J. M. & De Angelis, F. Ab initio

molecular dynamics simulations of methylammonium lead iodide perovskite degradation by water. _Chem. Mater._ 27, 4885–4892 (2015). Article  Google Scholar  * Zhang, C., Wu, J., Galli, G.

& Gygi, F. Structural and vibrational properties of liquid water from van der Waals density functionals. _J. Chem. Theory Comput._ 7, 3054–3061 (2011). Article  Google Scholar  * Lin, I.

C., Seitsonen, A. P., Tavernelli, I. & Rothlisberger, U. Structure and dynamics of liquid water from ab initio molecular dynamics-comparison of BLYP, PBE, and revPBE density functionals

with and without van der Waals corrections. _J. Chem. Theory Comput._ 8, 3902–3910 (2012). Article  Google Scholar  * Thomas, M., Brehm, M., Fligg, R., Vöhringer, P. & Kirchner, B.

Computing vibrational spectra from ab initio molecular dynamics. _Phys. Chem. Chem. Phys._ 15, 6608–6622 (2013). Article  Google Scholar  * Car, R. & Parrinello, M. Structural,

dynamical, and electronic properties of amorphous silicon: an ab initio molecular dynamics study. _Phys. Rev. Lett._ 60, 204–207 (1988). Article  Google Scholar  * Kresse, G. & Hafner,

J. Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium. _Phys. Rev. B_ 49, 14251–14269 (1994). Article  Google Scholar  * Johari, P.,

Qi, Y. & Shenoy, V. B. The mixing mechanism during lithiation of Si negative electrode in Li-ion batteries: an ab initio molecular dynamics study. _Nano Lett._ 11, 5494–5500 (2011).

Article  Google Scholar  * Lee, T. & Elliott, S. Ab initio computer simulation of the early stages of crystallization: application to Ge2Sb2Te5 phase-change materials. _Phys. Rev. Lett._

107, 145702 (2011). Article  Google Scholar  * Ong, S. P. et al. Phase stability, electrochemical stability and ionic conductivity of the Li10±1MP2X12 (M = Ge, Si, Sn, Al or P, and X = O, S

or Se) family of superionic conductors. _Energy Environ. Sci._ 6, 148–156 (2013). Article  Google Scholar  * Deng, Z., Radhakrishnan, B. & Ong, S. P. Rational composition optimization

of the lithium-rich Li3OCl1-_x_Br _x_ anti-perovskite superionic conductors. _Chem. Mater._ 27, 3749–3755 (2015). Article  Google Scholar  * Zhu, Z., Chu, I.-H., Deng, Z. & Ong, S. P.

Role of Na+ interstitials and dopants in enhancing the Na+ conductivity of the cubic Na3PS4 superionic conductor. _Chem. Mater._ 27, 8318–8325 (2015). Article  Google Scholar  * Ling, C.,

Zhang, R., Arthur, T. S. & Mizuno, F. How general is the conversion reaction in Mg battery cathode: a case study of the magnesiation of α-MnO2. _Chem. Mater._ 27, 5799–5807 (2015).

Article  Google Scholar  * Ling, C. & Suto, K. Thermodynamic origin of irreversible magnesium trapping in chevrel phase Mo6S8: importance of magnesium and vacancy ordering. _Chem.

Mater._ 29, 3731–3739 (2017). Article  Google Scholar  * Xu, M., Ding, J. & Ma, E. One-dimensional stringlike cooperative migration of lithium ions in an ultrafast ionic conductor.

_Appl. Phys. Lett._ 101, 031901 (2012). Article  Google Scholar  * Jalem, R. et al. Concerted migration mechanism in the Li ion dynamics of garnet-type Li7La3Zr2O12. _Chem. Mater._ 25,

425–430 (2013). Article  Google Scholar  * Meier, K., Laino, T. & Curioni, A. Solid-state electrolytes: revealing the mechanisms of Li-ion conduction in tetragonal and cubic LLZO by

first-principles calculations. _J. Phys. Chem. C_ 118, 6668–6679 (2014). Article  Google Scholar  * Panchmatia, P. M. et al. Oxygen defects and novel transport mechanisms in apatite ionic

conductors: combined 17O NMR and modeling studies. _Angew. Chem. Int. Ed._ 50, 9328–9333 (2011). Article  Google Scholar  * Deng, Y. et al. Structural and mechanistic insights into fast

lithium-ion conduction in Li4SiO4-Li3PO4 solid electrolytes. _J. Am. Chem. Soc._ 137, 9136–9145 (2015). Article  Google Scholar  * Saxton, M. J. Single-particle tracking: the distribution of

diffusion coefficients. _Biophys. J._ 72, 1744–1753 (1997). Article  Google Scholar  * Berglund, A. J. Statistics of camera-based single-particle tracking. _Phys. Rev. E_ 82, 011917 (2010).

Article  Google Scholar  * Michalet, X. Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium. _Phys. Rev. E_ 82,

041914 (2010). Article  Google Scholar  * Michalet, X. & Berglund, A. J. Optimal diffusion coefficient estimation in single-particle tracking. _Phys. Rev. E_ 85, 061916 (2012). Article

  Google Scholar  * Chitra, R. & Yashonath, S. Estimation of error in the diffusion coefficient from molecular dynamics simulations. _J. Phys. Chem. B_ 101, 5437–5445 (1997). Article 

Google Scholar  * Leetmaa, M. & Skorodumova, N. V. Mean square displacements with error estimates from non-equidistant time-step kinetic Monte Carlo simulations. _Comput. Phys. Commun._

191, 119–124 (2015). Article  Google Scholar  * Catti, M. Short-range order and Li+ ion diffusion mechanisms in Li5La9□2(TiO3)16 (LLTO. _Solid State Ion._ 183, 1–6 (2011). Article  Google

Scholar  * Lang, B., Ziebarth, B. & Elsässer, C. Lithium ion conduction in LiTi2(PO4)3 and related compounds based on the NASICON structure: a first-principles study. _Chem. Mater._ 27,

5040–5048 (2015). Article  Google Scholar  * Burbano, M., Carlier, D., Boucher, F., Morgan, B. J. & Salanne, M. Sparse cyclic excitations explain the low ionic conductivity of

stoichiometric Li7La3Zr2O12. _Phys. Rev. Lett._ 116, 135901 (2016). Article  Google Scholar  * Uebing, C. & Gomer, R. Determination of surface diffusion coefficients by Monte Carlo

methods: comparison of fluctuation and Kubo–Green methods. _J. Chem. Phys._ 100, 7759–7766 (1994). Article  Google Scholar  * Van der Ven, A., Ceder, G., Asta, M. & Tepesch, P.

First-principles theory of ionic diffusion with nondilute carriers. _Phys. Rev. B_ 64, 184307 (2001). Article  Google Scholar  * Murch, G. E. The Haven ratio in fast ionic conductors. _Solid

State Ion._ 7, 177–198 (1982). Article  Google Scholar  * Murugan, R., Thangadurai, V. & Weppner, W. Fast lithium ion conduction in garnet‐type Li7La3Zr2O12. _Angew. Chem. Int. Ed._ 46,

7778–7781 (2007). Article  Google Scholar  * Kamaya, N. et al. A lithium superionic conductor. _Nat. Mater._ 10, 682–686 (2011). Article  Google Scholar  * Arbi, K., Mandal, S., Rojo, J.

& Sanz, J. Dependence of ionic conductivity on composition of fast ionic conductors Li1+_x_Ti2-_x_Al _x_ (PO4)3, 0≤ _x_ ≤ 0.7. a parallel NMR and electric impedance study. _Chem. Mater._

14, 1091–1097 (2002). Article  Google Scholar  * Kosova, N., Devyatkina, E., Stepanov, A. & Buzlukov, A. Lithium conductivity and lithium diffusion in NASICON-type Li1+_x_Ti2–_x_Al _x_

(PO4)3 (_x_=0; 0.3) prepared by mechanical activation. _Ionics_ 14, 303–311 (2008). Article  Google Scholar  * Huang, R. et al. Direct observation of the full transition from ballistic to

diffusive Brownian motion in a liquid. _Nat. Phys._ 7, 576–580 (2011). Article  Google Scholar  * Ruppert, D. & Wand, M. P. Multivariate locally weighted least squares regression. _Ann.

Stat._ 22, 1346–1370 (1994). Article  Google Scholar  * Aono, H., Sugimoto, E., Sadaoka, Y., Imanaka, N. & Adachi, G.-y. Ionic-conductivity and sinterability of lithium titanium

phosphate system. _Solid State Ion._ 40/41, 38–42 (1990). Article  Google Scholar  * Li, Y. et al. Ionic distribution and conductivity in lithium garnet Li7La3Zr2O12. _J. Power Sources_ 209,

278–281 (2012). Article  Google Scholar  * Kresse, G. & Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_

54, 11169–11186 (1996). Article  Google Scholar  * Perdew, J. P., Ernzerhof, M. & Burke, K. Rationale for mixing exact exchange with density functional approximations. _J. Chem. Phys._

105, 9982–9985 (1996). Article  Google Scholar  * Jain, A. et al. A high-throughput infrastructure for density functional theory calculations. _Comput. Mater. Sci._ 50, 2295–2310 (2011).

Article  Google Scholar  * Jain, A. et al. Formation enthalpies by mixing GGA and GGA+U calculations. _Phys. Rev. B_ 84, 045115 (2011). Article  Google Scholar  * Jain, A. et al. Commentary:

the materials project: a materials genome approach to accelerating materials innovation. _Apl. Mater._ 1, 011002 (2013). Article  Google Scholar  * Plimpton, S. Fast parallel algorithms for

short-range molecular dynamics. _J. Comput. Phys._ 117, 1–19 (1995). Article  Google Scholar  * Adams, S. & Rao, R. P. Ion transport and phase transition in Li7−_x_La3(Zr2−_x_M _x_

)O12(M = Ta5+, Nb5+, _x_ = 0, 0.25). _J. Mater. Chem._ 22, 1426–1434 (2012). Article  Google Scholar  * Belsky, A., Hellenbrandt, M., Karen, V. L. & Luksch, P. New developments in the

Inorganic Crystal StructureDatabase (ICSD): accessibility in support of materials research and design. _Acta Crystallogr. Sect. B Struct. Sci._ 58, 364–369 (2002). Article  Google Scholar 

Download references ACKNOWLEDGEMENTS The authors acknowledge the support from Office of Naval Research (ONR) and from National Science Foundation under award No. 1550423. We thank Shuo Zhang

for testing the initial parameters of classical molecular dynamics simulations. This research used computational facilities from the University of Maryland supercomputing resources, the

Maryland Advanced Research Computing Center (MARCC), and the Extreme Science and Engineering Discovery Environment (XSEDE) supported by National Science Foundation Award No. DMR150038.

AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Department of Materials Science and Engineering, University of Maryland, College Park, MD, 20742, USA Xingfeng He, Yizhou Zhu, Alexander Epstein

 & Yifei Mo * Maryland Energy Innovation Institute, University of Maryland, College Park, MD, 20742, USA Yifei Mo Authors * Xingfeng He View author publications You can also search for

this author inPubMed Google Scholar * Yizhou Zhu View author publications You can also search for this author inPubMed Google Scholar * Alexander Epstein View author publications You can

also search for this author inPubMed Google Scholar * Yifei Mo View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS Y.M. conceived the project.

Y.M. and X.H. developed the computation methods and analyses, and X.H. performed the computation and analyses. Y.M. and X.H. wrote the manuscript. All authors contributed to the discussions

and revisions of the manuscript. CORRESPONDING AUTHOR Correspondence to Yifei Mo. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION

PUBLISHER'S NOTE: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. ELECTRONIC SUPPLEMENTARY MATERIAL SUPPLEMENTARY

MATERIAL RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and

reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes

were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If

material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain

permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS

ARTICLE He, X., Zhu, Y., Epstein, A. _et al._ Statistical variances of diffusional properties from ab initio molecular dynamics simulations. _npj Comput Mater_ 4, 18 (2018).

https://doi.org/10.1038/s41524-018-0074-y Download citation * Received: 05 November 2017 * Revised: 09 March 2018 * Accepted: 16 March 2018 * Published: 03 April 2018 * DOI:

https://doi.org/10.1038/s41524-018-0074-y SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative