In previous works we have reported the general concept of the TAPED approach24–26. Here, for the first time, we provide a detailed description of the underlying algorithms and the overall methodology. And most importantly, we expand the approach for the study of ion transport phenomena in complex systems such as polycrystalline or disordered ion conductors. Specifically, we illustrate the capabilities and reach of this newly implemented approach by examining two case studies. First, we focus on the grain boundary Li-ion transport phenomena considering a polycrystalline Li3ClO solid electrolyte. And second, we build on a kinetic Monte Carlo (KMC) model to compute diffusion coefficients as a function of charge carrier concentration in spinel LiTiS2, an electrode material that presents Li-vacancy disorder and two symmetrically non-equivalent sites.

**TAPED** vs. **DFT**

We validate this approach by thoroughly comparing ion trajectories and diffusion barriers obtained using TAPED with the DFT-based climbing image nudged elastic band (DFT-CI-NEB) method27. In total, we consider five Li-ion migration pathways in the anti-perovskite bulk Li3OCl structure (Fig. 1a-e) and three Li-ion hops in the spinel LiTiS2, with one, two, and three Li vacancies along the pathway (Fig. 1f-h). We find that the minimal electron density paths (MEDP) predicted by the TAPED approach are in very good agreement with the corresponding minimum energy paths (MEPs) for Li-ion migration yielded by the DFT-CI-NEB method. Specifically, the root mean square deviation (RMSD) between the MEDP and MEP for Li3OCl and LiTiS2 are 0.17 Å and 0.05 Å, respectively, with the greatest RMSD of 0.23 Å observed for the 3.9 Å Li hop in Li3OCl (Fig. 1c). Interestingly, when comparing the MEDP and MEP barriers (i.e., *Δρ* = *ρ**max*-*ρ**min* and *ΔE* = *E**max*-*E**min*, respectively) we find an excellent linear correlation (Figs. 2a and 2d). Furthermore, the electron density profile reproduces accurately the shape of the energy profile (Figs. 2b, 2c 2e and 2f). There are some divergences in the high energy barriers for Li3OCl, but the overall correlation is high, with coefficients of determination (*R**2*) of 0.953 and 0.996 for Li3OCl and LiTiS2, respectively.

To understand the origin of the observed MEDP deviations from the corresponding MEP trajectories and profiles, we need to consider that the TAPED approach neglects geometrical changes in the local atomic environment of the diffusing Li-ion. The approximation leads to an intrinsic error in the predicted geometry and barriers of the migration pathways. And this error is expected to become larger as the trajectory of the diffusing Li-ion is closer to nearby atoms of the structure, especially when there are anisotropic geometrical changes in the local atomic environment along the pathway. To illustrate this, consider a diffusing ion trajectory that goes through a narrow window that is smaller than the size of the ion (Fig. S1a, b). To squeeze through the window the ion should push the atoms out of their equilibrium positions. In the case of isotropic deformation, the bottleneck positions of the equilibrium and the deformed window coincide, while for anisotropic deformation, the positions differ. The same happens with the whole Li-ion migration trajectory. For example, the DFT-CI-NEB calculation of the 3.9 Å Li-ion diffusion (Fig. 1c) with frozen atomic coordinates yields a trajectory equivalent to the MEDP (Fig. S1d). In contrast, gradient paths calculated for the transition state geometry coincide with the DFT-CI-NEB trajectory (Fig. S1c). Note that this limitation is shared by all methods that neglect geometrical relaxation in the modeling of the diffusion process, like bond valence10,28,29, pinball8, pathfinder30, charge-based analysis31,32, or Voronoi-based methods33,34. However, even with this limitation, the results yielded by the TAPED approach are in good quantitative agreement with those obtained using the DFT-CI-NEB method. Furthermore, the TAPED-computed ion migration pathways, barriers, and shape of the energy profiles can readily be used to improve the effectiveness of DFT-CI-NEB simulations by guiding the setting of the initial position of each image along the diffusion paths; this is especially useful when tackling complex multi-step diffusion processes.

## Progress beyond the state of the art

One of the main disadvantages of NEB-based methods is a critical dependency on the initial guess of the diffusion pathway. Generally, the pathway geometry is approximated by linear interpolation between stable initial and final atomic positions. However, for some pathways (e.g., Fig. 1b, e) such crude interpolation can cause atom overlapping that leads to slow convergence or simulations to fail. To overcome this issue various algorithms exist, such as the pathfinder30 which preoptimizes the initial trajectory using an electrostatic force field. Yet, the problem of searching for the global MEP remains. In other words, a NEB optimization yields a local MEP, which can differ from the global MEP. This occurs when the initial pathway and global MEP are separated hills on the potential energy surfaces (PES) as illustrated in Fig. S2a. For example, when the NEB optimization, starts from the linear interpolation of one of the K-ion hops in K2SbPO6, the energy barrier is overestimated by about 7.5 eV due to this issue24. Thus, finding the global MEP requires scanning the entire PES, which is unfeasible at the DFT level.

One alternative strategy is to use simple interatomic force fields to approximate the complete PES. In this way, the MEPs can be found applying topological analysis of the field critical points (minima and adjacent saddle points) and the gradient lines connecting them. Follow-up DFT-CI-NEB calculations can then be used to refine such approximated MEPs, achieve better accuracies and maximize computational resources. For example, Adams *et al.* have implemented the bond valence pathway analyzer (BVPA) program for automatically finding ion migration paths based on scanning bond valence energy landscapes (BVELs).10 In this method, the critical points of a BVEL are calculated using a grid-based algorithm, and the path geometries determined by means of a modified Dewar, Healy, and Stewart algorithm35. To reduce the number of potential paths under consideration, they introduced geometrical restrictions, a cut-off distance of 3.5 Å between the minima and the saddle point, and an angle range of 90–270 degrees for the two minima and adjacent saddle points. These criteria are justified by the absence in available literature of reported ion jumps with longer distances between the site and the saddle point; and the assumption that straight-like hops are the most relevant, even if such pathways imply slightly higher energy barriers than other competing paths. However, the lack of reports is not solid proof that long ion hops can be systematically ruled out. And ignoring possible low-energy transition states contradicts the transition state theory36, according to which the hop rate depends only on the activation energy, temperature, and an effective frequency associated with vibration of the migrating ion in the direction of the transition state (Eq. 3 in the Methods section). Furthermore, the BVPA approach neglects the interaction between mobile ions; this can cause significant underestimation of migration energy barriers in non-dilute systems. For example, we show that there is almost twice difference in the Li-ion migration energy barriers in LiTiS2 with one and three vacancies along the pathway (Fig. 2f). This limitation applies to similar implementations of BV methods37,38. Replacing BEVL by electron density and avoiding geometrical constraints enable eliminating these disadvantages. Shen and co-workers proposed an alternative method based on topological charge density analysis that considers the electrostatic repulsion of mobile ions31. The method provides information about possible sites of mobile ions, but it does not allow to approximate pathways geometry and shape of the energy profile without involving expensive NEB calculations. Our TAPED approach eliminates the bottlenecks of the methods mentioned above and provides a straightforward way for finding the global MEP.

Another potential issue with the NEB method is the proper identification of transition states within multi-step, MEPs because accurate representations of energy profiles require many images. However, increasing the number of images drastically increases the computation time of NEB calculations. In this regard, the CI-NEB approach helps mitigate this issue and guarantees to find the transition state and corresponding activation energy using a minimal number of images. The CI-NEB algorithm forces the intermediate image that has maximum energy to climb the hill of PES along the MEP until it reaches the top. Specifically, the algorithm maximizes the energy of the climbing image towards the elastic band and minimizes it perpendicular to the band27. Nevertheless, in the case of an energy profile with more than one maximum (i.e., at least one intermediate local minimum between the initial and final states), there is a probability that the method will converge to a low-energy transition state (Fig. S2b). The best practical solution is to split the multi-step diffusion path into a series of individual steps, ensuring that each CI-NEB calculation only involves one transition state. In principle, the CI-NEB optimization of such elementary pathways only needs one intermediate image and, therefore, the spring forces are not applied, which enhances convergence. Yet, the prerequisite to split the diffusion pathways into suitable single-transition-state steps implies having a good approximation for the energy profile which is generally unknown. In this regard, the TAPED approach can provide all required information to split the pathways (i.e., the positions of the migrating ion at local minima and transition states), which can be obtained from the MEDP and the corresponding electron density profiles.

Overall, the TAPED approach efficiently solves the critical issues of NEB-based calculations, namely the choice of starting diffusion geometry and image sampling. Since the approach is much less computationally-intensive than DFT-CI-NEB, it can be used as an individual tool for screening for potential ion conductors. Unlike the methods based on visual analysis of maps and isosurfaces of electron density32,39,40 or BVEL28, the TAPED approach is fully automated. The only required input information are the crystallographic data and the type of migrating ion. In addition, the fact that the migration energy barriers and the maxima of electron density for the MEDPs are linearly correlated (Fig. 2a, d) enables the implementation of quantitative simulations of complex systems, as demonstrated in the following.

## Grain boundary Li-ion transport in polycrystalline Li3ClO solid electrolytes

We have considered four low-energy periodic grain boundary (GB) models of anti-perovskite Li3OCl generated using the aimsgb online service41 and geometrically optimized using DFT. The Li-ion migration pathways were then calculated between all neighboring Li sites, with a total of 1536, 1956, 1408, and 2988 computed pathways for Σ3(111), Σ3(112), Σ5(210), and Σ5(310) GB models, respectively. The diffusion barriers through the GBs were found by analysis of the Li-ion migration maps as described in the Methods section.

Figure 3 shows the Li-ion migration map of the Σ3(111) GB model, where we see high-energy migration barriers (0.63 eV) perpendicular to the interface (red pathways) and low-energy barriers within the grains. Notice how the lowest-energy barrier pathways (0.28 eV) form 1-periodic maps (blue); and these maps are connected to the 2-periodic map parallel to the GB by 0.29 eV pathways (violet). Since the grains have finite size this 2-periodic map corresponds to a thin high conductive layer close to the surface of the grains. However, polycrystalline Li-ion diffusion is limited by the high-energy barrier (0.63 eV) through GBs. Figs. S3-S5 show the Li-ion migration maps of the Σ3(112), Σ5(210), and Σ5(310) GB models. In general, we find that the Li-ion migration energy barriers through and along the GBs are consistently higher than that of the bulk structure (Table 1). This is in good agreement with previous MD calculations of Dawson *et. al.* 42 and available experimental data43,44. The limiting energy barriers for Σ3(111), 0.63 eV, and Σ3(112), 0.57 eV, are very close to the experimental value of 0.59 eV for polycrystalline Li3OCl. Note that the concentration of Σ5 GBs and, therefore, their contribution to the total Li-ion diffusion characteristics should be much less than Σ3 GBs because Σ5 GBs are less energetically favorable than Σ3 GBs42. Thus, in the following we focus the discussion on the results for the Σ3(111) and Σ3(112) models.

Table 1

Computed Li-ion migration energy barriers for different GBs models and bulk Li3OCl. Experimental values are also provided.

| Grain boundary | Bulk |

| **Σ3(111)** | **Σ3(112)** | **Σ5(210)** | **Σ5(310)** |

**E (eV) this work** | 0.63 | 0.57 | 0.46 | 0.47 | 0.38 |

**E (eV)**42 | 0.52 | 0.56 | 0.40 | 0.49 | 0.29 |

**E****exp.** **(eV)**43 | 0.59 | 0.35 |

Further analysis of the Σ3(112) GB model leads to a more fundamental understanding of the migration mechanism. Its migration map reveals pathways with relatively low-energy barriers within the interface region (Fig. 4a, S3): 0.23 eV (blue) and 0.25 eV (purple). Therefore, one would expect an increase in ionic conductivity at the GB with respect to the bulk. However, in practice these pathways cannot be realized due to a special structural feature of the Σ3(112) GB. Notice the electron density profile of the pathway that crosses the interface (Fig. 4a) has a concave shape (Fig. 4b). This is a consequence of the repulsive forces between neighboring Li ions located in the same cage of the framework. Thus, the formation of a vacancy on one side of the pathway forces a Li-ion on the opposite side to fell into the potential energy well, which is in the center of the cage (site II in Fig. 4b, c). The geometrical DFT optimization of the Σ3(112) GB model with the vacant Li site leads to the displacement of the opposite Li to the interface (Fig. 4c). Since this position is energetically more favorable than the initial site, the effective energy barrier for a diffusing Li-ion to cross the GB is 0.79 eV, as depicted in Fig. 4b. Among the alternative diffusion routes through the GB (Fig. S3), we find the lowest barrier is 0.57 eV. Thus, as in the case of Σ3(111) GBs, the Σ3(112) GBs also exhibits higher ionic resistance than the bulk.

It is important to highlight that the identified redistribution of Li ions and vacancies at the Σ3(112) GBs should also be observed in MD simulations. However, the detection of this effect might be cumbersome because it requires a detailed analysis of MD trajectories. For example, the MD investigation of Li3OCl GBs by Dawson *et al.* 42 mainly focused on Li-ion diffusion properties but omitted a complementary analysis capable of revealing such subtle dynamics. In contrast to MD, the TAPED approach provides comprehensive information on all possible Li-ion diffusion trajectories inside and between the grains, as well as values of migration barriers for each pathway. This approach can systematically examine the effect of the structural features of GB on ionic conductivity, which could help gain a better understanding of ion transport phenomena in polycrystalline ionic conductors. Moreover, from an efficiency viewpoint, it is important to point out that we explicitly accessed here the Li-ion migration energy barriers of over a thousand of non-equivalent trajectories for a range of GB models − a problem that would be intractable using a plain DFT-CI-NEB approach. And for this, we only needed to compute five DFT-CI-NEB pathways in a 2×2×2 supercell of the bulk Li3OCl to obtain the linear correlation shown in Fig. 2a for converting density barriers into energy barriers. The good agreement between our results and alternative MD simulations and excellent agreement with experiment (Table 1) indicates that the approach can be quantitative and scalable.

**Coupling TAPED with KMC to compute Li-ion diffusion coefficients in LiTiS** **2** **electrodes as a function of Li content**

Now we consider the application of the TAPED approach for computing Li-ion diffusion coefficients in LiTiS2 using KMC simulations. For this task, we first need to explore the complete configurational space of the Li*x*TiS2 system as a function of Li content (0 ≤ *x* ≤ 1), which was performed by discretizing the Li content using Δ*x* steps of 0.62. This corresponds to a total of 531 Li*x*TiS2 structures with symmetrically unique Li-vacancy arrangements, which were relaxed using DFT. And, to build the KMC model, we have applied the TAPED approach to compute migration trajectories and migration energy barriers for 12504 Li-ion hops between the closest pairs of occupied and vacant Li sites. The diffusion coefficients were averaged over 100 independent KMC trajectories, including 107 steps.

The distribution of Li-ion migration energy barriers shows a diverse dependency with the Li content (Fig. 5a). This is due to a strong sensitivity of the migration barriers to the local atomic environment. The distribution has three distinct energy ranges corresponding to 1-, 2- and 3-vacancy Li-ion migration mechanisms45. Regardless of the Li content, the highest and lowest barrier pathways correspond to the 1- and 3-vacancies ones, respectively. Obviously, the number of possible low barrier pathways increases as the Li content decreases. Consequently, diluted systems should yield the highest ionic conductivities. This is confirmed by the calculated diffusion coefficients (Fig. 5b). Starting from ~ 10− 12.4 cm2/s at a Li content of 6.25%, the trace- and jump-diffusion coefficients reached their maxima (~ 10− 11.9 cm2/s and ~ 10− 11.4 cm2/s, respectively) at a Li content of 25%. Within this Li range (less than 25%), there are significant changes in the unit cell volume (Fig. 5c). The lowering of the diffusion coefficient as the Li content approaches to zero is related to the contraction of the unit cell. A further increase in the Li content from 25–93.75% leads to a decrease in the diffusion coefficient to ~ 10− 18.4 cm2/s and ~ 10− 17.3 cm2/s for trace and jump diffusion, respectively. The swift decline of the diffusion coefficients at 50% and 80% Li are associated with the significant reduction of available 3- and 2-vacancies Li hops, respectively (Fig. 5a).

These results are consistent with the previous cluster-expansion-based KMC simulation of the spinel Li*x*TiS2 by Bhattacharya and Van der Ven45. However, we find a systematic difference of ~ 0.2–0.3 eV between the energy barriers calculated here (Fig. 2f) and in the previous study which we attribute to the underestimation of the energy barriers by the standard NEB method used in the previous work. This also leads to a deviation in the values and shape of the curves of the diffusion coefficients as a function of the Li content (Fig. 5b, d). Specifically, the diffusion coefficients calculated by Bhattacharya and Van der Ven are several orders of magnitude higher than ours because of the mismatch in the energy barriers. And their dependencies with the Li content are smoother than ours due to the smaller energy gaps between 1-, 2-, and 3-vacancies pathways. Nevertheless, the underlying trend of the dependency and following conclusions agree qualitatively. Thus, the TAPED approach can be an alternative to the cluster-expansion-based method for exploring Li-ion migration energy barriers in the configuration space of disordered structures, requiring a significantly lower number of NEB-based computationally-intensive calculations. Furthermore, we put forward that combining TAPED with KMC can allow the computation of Li-ion diffusion coefficients in complex systems. This opens new perspectives in assessing the impact varying the charge carrier concentration has not only in pristine bulk materials, but also in GBs, doped materials, and other simulations of practical interest.