Enhancing first-principles simulations of complex solid-state ion conductors using topological analysis of procrystal electron density

doi:10.21203/rs.3.rs-1742415/v1

Download PDF

Article

Enhancing first-principles simulations of complex solid-state ion conductors using topological analysis of procrystal electron density

https://doi.org/10.21203/rs.3.rs-1742415/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 30 Aug, 2022

Read the published version in npj Computational Materials →

You are reading this latest preprint version

Atomistic-level understanding of ion migration mechanisms holds the key to design high-performance solid-state ion conductors for a breadth of electrochemical devices. First-principles simulations play an important role in this quest. Yet, these methods are generally computationally-intensive, with limited access to complex, low-symmetry structures, such as interfaces. Here we show how topological analysis of the procrystal electron density can be applied to efficiently mitigate this issue. We discuss how this methodology goes beyond current state of the art capabilities and demonstrate this with two examples. In the first, we examine Li-ion transport across grain boundaries in Li₃ClO electrolyte. Then we compute diffusion coefficients as a function of charge carrier concentration in spinel LiTiS₂ electrode material. These two case studies do not exhaust the opportunities and might constitute motivations for still more complex applied materials.

Solid-state ionic conductors find application in many materials science domains. For example, they are deemed to play a key role in the swiftly evolving energy storage and conversion technologies¹, as well as the development of new functional materials for sensing^2,3. Understanding and controlling ionic transport in these systems is essential for their successful application in high-performance electrochemical devices. This typically involves identifying ion occupancies, ion migration paths, charge carrier types, point defects and ion-doping sites at both bulk and interfacial structures such as grain boundaries and heterophases^4,5.

Quantum mechanics is a well-established theoretical approach for accessing these aspects at the atomistic level. In particular, density functional theory (DFT) offers a good trade-off in terms of accuracy and efficiency, with a successful track record^6,7. Yet, gathering in-depth information about ion migration processes and their underlying mechanisms in complex systems requires exhaustive explorations of tortuous potential energy landscapes or the evaluation of atomic dynamics, which are computationally expensive tasks. This limits the use of DFT-based calculations in high-throughput screening campaigns to discover new ionic conductors or the assessment of the role of charge carrier concentration and interfacial structures in performance. And less computationally intensive alternatives to DFT-based computations^8–10 often come at the expense of structural detail, quantitative predictability, and transferability.

Workflows for the computation of kinetic properties (i.e. ion diffusivity) require not only finding the lowest energy geometries of initial and final states, but also considering all plausible migration paths and computing their transition states. And, from a general perspective, such workflows should consider multi-step diffusion mechanisms and the role of possible metastable sties. Current DFT-based methods are well suited to perform such tasks in relatively simple bulk materials and provide accurate predictions^7,11. The two most successful approaches in this context are ab initio molecular dynamics (AIMD) simulations^12–15 and nudged elastic band (NEB) methods^16–18. Despite important advances to enhance these approaches by, for example, exploiting crystal symmetry^19,20 or using machine learning techniques^21–23, the main bottleneck associated to their high computational time remains. And this limitation is largely amplified when dealing with complex, low-symmetry systems (such as interfaces) and the systematic exploration of large compositional and configurational spaces, all those key aspects for truly applied materials development.

Here, we implement a topological analysis of procrystal electron density (TAPED) to compute the kinetic properties of ion conductors effectively and efficiently. Compared with conventional AIMD and NEB techniques, this approach requires much less computational time, while allowing for the consideration of interface models and large supercells containing thousands of symmetrically nonequivalent ion migration pathways. Furthermore, the implemented algorithm is fully automated and requires minimal user involvement. We show that the approach can be calibrated to yield quantitative agreement with DFT-based NEB calculations by modeling two representative case systems (one solid electrolyte and one electrode material) for which high quality theoretical and experimental data is available. Accurate estimates of the ionic kinetic properties can be achieved in both cases at a relatively low computational cost.

In previous works we have reported the general concept of the TAPED approach^24–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 Li₃ClO 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 LiTiS₂, 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) method²⁷. In total, we consider five Li-ion migration pathways in the anti-perovskite bulk Li₃OCl structure (Fig. 1a-e) and three Li-ion hops in the spinel LiTiS₂, 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 Li₃OCl and LiTiS₂ are 0.17 Å and 0.05 Å, respectively, with the greatest RMSD of 0.23 Å observed for the 3.9 Å Li hop in Li₃OCl (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 Li₃OCl, but the overall correlation is high, with coefficients of determination (R²) of 0.953 and 0.996 for Li₃OCl and LiTiS₂, 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 valence^10,28,29, pinball⁸, pathfinder³⁰, charge-based analysis^31,32, or Voronoi-based methods^33,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 pathfinder³⁰ 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 K₂SbPO₆, the energy barrier is overestimated by about 7.5 eV due to this issue²⁴. 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).¹⁰ 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 algorithm³⁵. 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 theory³⁶, 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 LiTiS₂ with one and three vacancies along the pathway (Fig. 2f). This limitation applies to similar implementations of BV methods^37,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 ions³¹. 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 band²⁷. 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 density^32,39,40 or BVEL²⁸, 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 Li₃ClO solid electrolytes

We have considered four low-energy periodic grain boundary (GB) models of anti-perovskite Li₃OCl generated using the aimsgb online service⁴¹ 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. ⁴² and available experimental data^43,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 Li₃OCl. 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 GBs⁴². 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 Li₃OCl. Experimental values are also provided.
	Grain boundary				Bulk
	Σ3(111)	Σ3(112)	Σ5(210)	Σ5(310)	Bulk
E (eV) this work	0.63	0.57	0.46	0.47	0.38
E (eV)⁴²	0.52	0.56	0.40	0.49	0.29
E_exp. (eV)⁴³	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 Li₃OCl GBs by Dawson et al. ⁴² 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 Li₃OCl 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 ₂ electrodes as a function of Li content

Now we consider the application of the TAPED approach for computing Li-ion diffusion coefficients in LiTiS₂ using KMC simulations. For this task, we first need to explore the complete configurational space of the Li_xTiS₂ 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_xTiS₂ 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 10⁷ 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 mechanisms⁴⁵. 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 cm²/s at a Li content of 6.25%, the trace- and jump-diffusion coefficients reached their maxima (~ 10^− 11.9 cm²/s and ~ 10^− 11.4 cm²/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 cm²/s and ~ 10^− 17.3 cm²/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_xTiS₂ by Bhattacharya and Van der Ven⁴⁵. 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.

Calculation of ionic diffusion trajectories and migration barriers using TAPED

The TAPED method approximates the PES by a field to obtain an analytical or numerical function of the potential energy from coordinates of the migrating ion. Based on this representation of the PES, it is possible to identify all the MEPs between the initial and final positions of the migrating ion and find the global MEP. In particular, the equilibrium and metastable sites of the ion are associated with the minima of the field, so the MEPs correspond to gradient paths connecting these minima. These gradient paths starting from minima end in saddle points representing the transition state between two adjacent minima. The value of the field in the saddle point characterizes the migration barrier. The MEP can include one or more transition states. Given an equilibrium structure, the ratio between the number (N) of the saddle points, minima, and gradient paths within the MEP is N : N + 1 : 2N, with the only exception of post-transition state bifurcations⁴⁶. Finding the global MEP requires the iteration over all MEPs to identify the one with the lowest values of the barrier. Due to the periodicity of the crystal structure, the number of MEPs is infinite, but the number of unique transition states is finite. Thus, this task can be solved by a simple graph-based algorithm, as illustrated in Fig. 6. Initially, the topology of the field is represented as an edge weighted graph in which vertices, edges, and weights correspond to the critical points, gradient paths, and field value at the saddle points. Next, the ends of the MEPs are associated with the closest minima, saddle points, or points of the gradient paths. Accounting for the last two involves considering non-equilibrium structures, without available geometry optimization. An example is I-I´ pathways in Fig. 4. In the case of equilibrium structures, the initial and final position of the migrating ion matches with the minima of the PES. Further steps of the algorithm search for the route with the lowest value of the barrier (Fig. 6c), which is solved using a Dijkstra-like algorithm^47,48. Starting from one of the two given vertices, the algorithm iterates over all possible routes, discarding variants with a higher weight than the most optimal route found in the previous steps. As a result, we obtain the ion migration trajectory, the field profile of the pathway, as well as the positions of local minima and transition states (Fig. 6d). The detailed information about the implementation of the approach can be found in the IonExplorer2 code (https://github.com/mcs-cice/IonExplorer2). The calculation of the gradient paths and the critical points was carried out using the Critic2 program⁴⁹. Visualization of ion diffusion trajectories and crystal structures was performed by means of Vesta⁵⁰.

To approximate the PES, we use procrystal electron density³⁹, which yields good agreement with experimental⁵¹ and ab initio derived electron density distributions⁵². The process of ion migration is considered as the movement of a single negative point charge in the constant electron field of the structure (Fig. 6b). The approximation assumes that the main forces that affect the ion trajectory are the electrostatic repulsion of electron shells of the ion and atoms of the structure. In this case, the equilibrium and metastable ion positions are defined as minima of electron density: the center of intrastructural cages, namely cage critical points (CCP). The positions of the migrating ion at transition states are in the maxima of electron density along the gradient path: the center of the atomic window, namely ring critical points (RCP). Notice that the algorithm to calculate the procrystal electron density distribution only requires information about the atomic types and coordinates.

Choice of the initial and final positions for calculating the diffusion trajectory

One of the input data required for calculating the lowest-barrier ion migration pathway is the initial and final positions of the ion. In ion conductors, mobile ions hop between nearest vacant atomic positions. In other words, the possible final positions of any diffusion trajectory should correspond to the neighbor sites of the ions. A common practice to choose the final positions of the pathways is by simple visual analysis of the structure or subjective distance cutoff criteria. However, such approach cannot be automated for calculation of a large data set and complex geometries. Alternatively, potential ion hops can be found using information about neighbor sites derived from the Voronoi partition of the ionic sublattice³³. For a set of M geometrical objects in Euclidean space, the Voronoi partition divides the space into M regions (domains). Each domain surrounds one of the objects so that the distance from any point of the domain to the enclosed object is no greater than the distance to any other object of the set⁵³. Thus, the domain allocates a region of space that belongs to the object. Depending on the geometry and arrangement of the objects, Voronoi domains take various shapes. In the case of ion sublattice approximated by a set of points or spheres with an equal radius, the atomic domains are convex polyhedra (Fig. S6a). Atoms whose domains share common faces, edges, or vertices are neighbors (Fig. S6b). Nearest (direct) neighbors are adjacent by faces, while more distant (indirect neighbors) by edges or vertices of the polyhedra. In most cases, ion hops between direct neighbor positions are the most energetically favorable ones, but there are exceptions. For the high-throughput screening of ion conductors, considering only hops between direct neighbors can significantly decrease the computational cost. However, the detailed analysis of the diffusion requires the calculation of all possible unique trajectories. The whole set of neighbor pairs forms a graph that reflects the map of potential ion hops (Fig. S6c). Due to the periodicity of the crystal structure, this graph is always 3-periodic. Depending on the distribution of ions, graph vertices can have a quite high degree. For example, each atom of the Li sublattice in the Li₃OCl structure has 8 direct and 26 indirect neighbors (Fig. S6b). This means that in the unit cell with 3 Li atoms there are about 100 possible hops for neighbor sites. However, among all of them, there are only 5 symmetrically nonequivalent ones. Thus, taking symmetry into account can drastically decrease the numbers of calculations of the lowest-barrier ion migration trajectories.

In summary, the implemented searching scheme of initial and final position migration pathways include 3 main steps:

Calculation of Voronoi partition for the sublattice of ions involved in diffusion process.

Search for the pairs of ions whose Voronoi polyhedra are adjacent by faces, edges, or vertices and build the graph of neighbors.

Identification of symmetrically unique edges of the graph.

In the case of studying ion conductors at various levels of ion content (e.g., Li_xTiS₂ system), the graph of possible hops is calculated for a structure with fully occupied positions. The graph for a partially occupied structure is then generated by removing all edges connecting two occupied or two vacant positions.

Analysis of ion migration map

The presence of low barrier migration pathways is not a sufficient criterion for high ionic conductivity. Since ion transport is limited by the highest barrier along the diffusion path, low barrier pathways should form a periodic migration map allowing ions to diffuse through the structure. The preference for certain diffusion crystallographic directions and their associated barriers can be found based on analysis of the ion migration map. For this purpose, we iteratively removed pathways in decreasing order of the barriers until the map loses periodicity. At each iteration, nets of pathways and their periodicity are determined, and all non-periodic parts of the migration map are removed. The maximal barrier at the step before the map loses the periodicity in a certain direction corresponds to the diffusion barrier in this direction. For example, in the case of LiTiS₂ (Fig. 7a, b), the migration map of Li hops between neighboring sites includes three symmetrical nonequivalent pathways. The density barriers calculated by the TAPED approach equal to 18.7 a.u., 7.7 a.u., and 3.0 a.u. Removing first two pathways does not change the periodicity of the map (Fig. 7a); thus, the lowest barrier Li migration map of the LiTiS₂ structure is 3-periodic and formed by one type of the pathway. And the diffusion barrier in all directions is equal to 3.0 a.u. In the case of structures with lower symmetry the algorithm can pass through 2- and 1-periodic maps. Moreover, there is a possibility of splitting the migration map into several independent subnets of pathways, which may have equal or lesser periodicity relative to the parent net. The migration maps at the step before changing the periodicity or splitting the net of pathway are milestones, which provide comprehensive information about diffusion barriers in all crystallographic directions. Furthermore, this iterative procedure shows the evolution of the topology of the ion migration map depending on the amount of energy transferred to the system to overcome the diffusion barriers, for example, by heating or applying an electric potential.

The periodicity of the net is determined by the number of independent translation vectors. These vectors can be found using the breadth-first search algorithm. For this purpose, we represent the net using a periodic graph^54,55. In this notation, each vertex (P_i) of the graph is associated with its translation image by the (r,s,t) symbol, which indicates the ra + sb + tc translation from P_i(0,0,0), where a, b, and c are basis translation vectors (Fig. 7c). Starting from some arbitrary vertex, the algorithm traverses graph vertices by the neighbors. Before moving on the next vertex, the algorithm explores all vertices at the current depth level. Once a vertex P_i is traversed the algorithm does not explore the translational equivalent images P_i(r,s,t), but keeps the information about the vectors between these vertices. The procedure stops after traversing all translationally non-equivalent vertices.

Diffusion coefficient and KMC simulations

The KMC methods are based on stochastic simulations of dynamic processes within a predefined set of states of the system and transition rates between them⁵⁶. Detailed information of the methodology for calculating ion diffusion coefficients in undiluted systems can be found in Refs. ^45,57 by Van der Ven and co-workers. Specifically, we use here the Bortz-Kalos-Lebowitz algorithm^58,59 to estimate ionic diffusion coefficients and ion hop frequencies, where the selection of a transition to a new state and time increment is computed using the following two formulas:

$$\frac{1}{{{\Gamma }}_{tot}}\sum _{m=1}^{k-1}{{\Gamma }}_{m}<p\le \frac{1}{{{\Gamma }}_{tot}}\sum _{m=1}^{k}{{\Gamma }}_{m}$$

where Γ_tot is the sum of individual mth hop frequencies, Γ_m, and p is a random number from 0 to 1; and

$$\varDelta t=-\frac{1}{{{\Gamma }}_{tot}}\text{ln}\xi$$

where ∆t and ξ are the time increment and a random number from 0 to 1, respectively.

The hop frequencies are calculated as:

$${{\Gamma }}_{m}={{\nu }}^{*}exp\left(\frac{-\varDelta E}{{k}_{B}T}\right)$$

where ΔE is the energy barrier for the hop of the ion and ν* is the vibrational pre-exponential factor. Here we assumed ν*=10¹³ Hz as a commonly accepted approximate value.

We compute ΔE using the TAPED approach. To convert TAPED derived density barriers into actual energy barriers, we use a linear equation of Δρ vs. ΔE (Fig. 2d) obtained through a combination of DFT calculation and the CI-NEB method, considering a few different diffusion pathways. For the calculation of ion self-diffusion from KMC trajectories, we use the following two formulas for jump- and trace-diffusion coefficients:

$${D}_{J}=\underset{t\to \infty }{\text{lim}}\left[\frac{1}{2dt}⟨\frac{1}{N}{\left(\sum _{i=1}^{N}{\overrightarrow{r}}_{i}\left(t\right)\right)}^{2}⟩\right]$$

and

$${D}^{*}=\underset{t\to \infty }{\text{lim}}\left[\frac{1}{2dt}\left(\frac{1}{N}\sum _{i=1}^{N}⟨{\left[{\overrightarrow{r}}_{i}\left(t\right)\right]}^{2}⟩\right)\right]$$

where t stands for time, d is the dimension of the ion migration map, N is the number of ions involved in the diffusion process, and ${\overrightarrow{r}}_{i}$ is the vector that connects the position of the ith ion at initial (t₀) and final (t) times.

DFT calculations

The DFT calculations were performed using the Vienna ab initio simulation package (VASP) 5.4 code^60–62. We use the Perdew − Burke − Ernzerhof (PBE) exchange − correlation functional⁶³ to calculate total energies and electron densities, together with a plane-wave kinetic energy cutoff of 400 eV and a Γ-point k-point sampling. The Li (1s, 2s, 2p), S (3s, 3p), Cl (3s, 3p), O (2s, 2p), Ti (3p, 4s, 3d) electrons were treated explicitly. The remaining electrons were replaced by PBE-based projector-augmented wave potentials^64,65. Structural optimization of ground states was performed by allowing the lattice parameters, cell shape, and atomic positions to relax until the maximum force on each atom was less than 0.05 eV/Å. The CI-NEB calculations were carried out employing the transition state tools (TST) for VASP (https://theory.cm.utexas.edu/vtsttools/references.html).

DATA AVAILABILITY

All the data supporting the findings of this work are available from the authors upon reasonable request. The DFT-relaxed structures of the Li₃OCl (upload Id: YgSjkRuvTduzU1HSpX4Ymg) and LiTiS₂ (upload Id: 9exhy71TSwKctviNKmYAYA) systems can be found in the NOMAD public repository https://nomad-lab.eu/.

CODE AVAILABILITY

The program implementation of algorithms reported in this work is available at https://github.com/mcs-cice/IonExplorer2.

ACKNOWLEDGEMENTS

This work was supported by the European Commission through the H2020 program under Grant agreement number 875028 (SUBLIME Project). The authors gratefully acknowledge the technical support provided by the Barcelona Supercomputing Center and the computer resources from Calendula (RES-QS-2021-1-0003 and RES-QS-2021-2-0004) and SGI/IZO-SGIker UPV/EHU.

AUTHOR CONTRIBUTIONS

Andrey Golov: Conceptualization, methodology, investigation, writing – original draft, writing – review and editing. Javier Carrasco: Conceptualization, funding acquisition, supervision, writing – review and editing.

COMPETING INTERESTS

The authors declare no competing interests.

ADDITIONAL INFORMATION

No additional information.

Zhao, Q., Stalin, S., Zhao, C.-Z. & Archer, L. A. Designing solid-state electrolytes for safe, energy-dense batteries. Nature Reviews Materials 5, 229–252 (2020).
Liu, T., Wang, X., Li, L. & Yu, J. Review—Electrochemical NO_x Gas Sensors Based on Stabilized Zirconia. Journal of The Electrochemical Society 164, B610–B619 (2017).
Gorbova, E. et al. Fundamentals and Principles of Solid-State Electrochemical Sensors for High Temperature Gas Detection. Catalysts 12, 1 (2022).
Famprikis, T., Canepa, P., Dawson, J. A., Islam, M. S. & Masquelier, C. Fundamentals of inorganic solid-state electrolytes for batteries. Nature Materials 18, 1278–1291 (2019).
Augustyn, V., McDowell, M. T. & Vojvodic, A. Toward an Atomistic Understanding of Solid-State Electrochemical Interfaces for Energy Storage. Joule 2, 2189–2193 (2018).
Jain, A., Shin, Y. & Persson, K. A. Computational predictions of energy materials using density functional theory. Nature Reviews Materials 1, 15004 (2016).
Van der Ven, A., Deng, Z., Banerjee, S. & Ong, S. P. Rechargeable Alkali-Ion Battery Materials: Theory and Computation. Chemical Reviews 120, 6977–7019 (2020).
Kahle, L., Marcolongo, A. & Marzari, N. Modeling lithium-ion solid-state electrolytes with a pinball model. Physical Review Materials 2, 065405 (2018).
Katcho, N. A. et al. An investigation of the structural properties of Li and Na fast ion conductors using high-throughput bond-valence calculations and machine learning. Journal of Applied Crystallography 52, 148–157 (2019).
Wong, L. L. et al. Bond Valence Pathway Analyzer—An Automatic Rapid Screening Tool for Fast Ion Conductors within softBV. Chemistry of Materials 33, 625–641 (2021).
Franco, A. A. et al. Boosting Rechargeable Batteries R&D by Multiscale Modeling: Myth or Reality? Chemical Reviews 119, 4569–4627 (2019).
Jalem, R. et al. Concerted Migration Mechanism in the Li Ion Dynamics of Garnet-Type Li₇La₃Zr₂O₁₂. Chemistry of Materials 25, 425–430 (2013).
Wang, Y. et al. Design principles for solid-state lithium superionic conductors. Nature Materials 14, 1026–1031 (2015).
He, X., Zhu, Y., Epstein, A. & Mo, Y. Statistical variances of diffusional properties from ab initio molecular dynamics simulations. npj Computational Materials 4, 18 (2018).
Kahle, L., Marcolongo, A. & Marzari, N. High-throughput computational screening for solid-state Li-ion conductors. Energy & Environmental Science 13, 928–948 (2020).
Carrasco, J. Role of van der Waals Forces in Thermodynamics and Kinetics of Layered Transition Metal Oxide Electrodes: Alkali and Alkaline-Earth Ion Insertion into V₂O₅. The Journal of Physical Chemistry C 118, 19599–19607 (2014).
Deng, Z., Radhakrishnan, B. & Ong, S. P. Rational Composition Optimization of the Lithium-Rich Li₃OCl_1–xBr_xAnti-Perovskite Superionic Conductors. Chemistry of Materials 27, 3749–3755 (2015).
Liu, M. et al. Spinel compounds as multivalent battery cathodes: a systematic evaluation based on ab initio calculations. Energy & Environmental Science 8, 964–974 (2015).
Bölle, F. T. et al. Autonomous Discovery of Materials for Intercalation Electrodes. Batteries & Supercaps 3, 488–498 (2020).
Bölle, F. T., Bhowmik, A., Vegge, T., Maria García Lastra, J. & Castelli, I. E. Automatic Migration Path Exploration for Multivalent Battery Cathodes using Geometrical Descriptors. Batteries & Supercaps 4, 1516–1524 (2021).
Lombardo, T. et al. Artificial Intelligence Applied to Battery Research: Hype or Reality? Chemical Reviews (2021) doi:10.1021/acs.chemrev.1c00108.
Jalem, R. et al. Bayesian-Driven First-Principles Calculations for Accelerating Exploration of Fast Ion Conductors for Rechargeable Battery Application. Scientific Reports 8, 5845 (2018).
Fujimura, K. et al. Accelerated Materials Design of Lithium Superionic Conductors Based on First-Principles Calculations and Machine Learning Algorithms. Advanced Energy Materials 3, 980–985 (2013).
Zolotarev, P. N., Golov, A. A., Nekrasova, N. A. & Eremin, R. A. Topological analysis of procrystal electron densities as a tool for computational modeling of solid electrolytes: A case study of known and promising potassium conductors. AIP Conference Proceedings 2163, 020007 (2019).
Golov, A. & Carrasco, J. Molecular-Level Insight into the Interfacial Reactivity and Ionic Conductivity of a Li-Argyrodite Li₆PS₅Cl Solid Electrolyte at Bare and Coated Li-Metal Anodes. ACS Applied Materials & Interfaces 13, 43734–43745 (2021).
Borisov, A. S. et al. Expanding the family of mineral-like anhydrous alkali copper sulfate framework structures: new phases, topological analysis and evaluation of ion migration potentialities. Journal of Applied Crystallography 54, 237–250 (2021).
Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. Journal of Chemical Physics 113, 9901–9904 (2000).
Adams, S. From bond valence maps to energy landscapes for mobile ions in ion-conducting solids. Solid State Ionics 177, 1625–1630 (2006).
Chen, H., Wong, L. L. & Adams, S. SoftBV – a software tool for screening the materials genome of inorganic fast ion conductors. Acta Crystallographica Section B 75, 18–33 (2019).
Rong, Z., Kitchaev, D., Canepa, P., Huang, W. & Ceder, G. An efficient algorithm for finding the minimum energy path for cation migration in ionic materials. The Journal of Chemical Physics 145, 074112 (2016).
Shen, J.-X., Li, H. H., Rutt, A. C., Horton, M. K. & Persson, K. A. Rapid discovery of cathodes, ionic conductors and solid-state electrolytes through topological migration analysis. (2022) arXiv:2202.00222.
Filsø, M. Ø. et al. Visualizing Lithium-Ion Migration Pathways in Battery Materials. Chemistry - A European Journal 19, 15535–15544 (2013).
Gulino, V. et al. Combined DFT and geometrical–topological analysis of Li-ion conductivity in complex hydrides. Inorganic Chemistry Frontiers 7, 3115–3125 (2020).
He, B. et al. CAVD, towards better characterization of void space for ionic transport analysis. Scientific Data 7, 153 (2020).
Dewar, M. J. S., Healy, E. F. & Stewart, J. J. P. Location of transition states in reaction mechanisms. Journal of the Chemical Society, Faraday Transactions 2 80, 227-233 (1984).
Vineyard, G. H. Frequency factors and isotope effects in solid state rate processes. Journal of Physics and Chemistry of Solids 3, 121–127 (1957).
He, B. et al. High-throughput screening platform for solid electrolytes combining hierarchical ion-transport prediction algorithms. Scientific Data 7, 151, (2020).
He, B. et al. A highly efficient and informative method to identify ion transport networks in fast ion conductors. Acta Materialia 203, 116490, (2021).
Filsø, M. Ø., Eikeland, E. & Iversen, B. B. Procrystal analysis as a tool for the visualization of ion migration pathways. AIP Conference Proceedings 1765, 020010 (2016).
Liu, Y., Jiang, X., Zhao, J. & Hu, M. Electronic charge density as a fast approach for predicting Li-ion migration pathways in superionic conductors with first-principles level precision. Computational Materials Science 192, 110380, (2021).
Cheng, J., Luo, J. & Yang, K. Aimsgb: An algorithm and open-source python library to generate periodic grain boundary structures. Computational Materials Science 155, 92–103 (2018).
Dawson, J. A., Canepa, P., Famprikis, T., Masquelier, C. & Islam, M. S. Atomic-Scale Influence of Grain Boundaries on Li-Ion Conduction in Solid Electrolytes for All-Solid-State Batteries. J Am Chem Soc 140, 362–368 (2018).
Lü, X. et al. Antiperovskite Li₃OCl Superionic Conductor Films for Solid-State Li-Ion Batteries. Advanced Science 3, 1500359 (2016).
Zhao, Y. & Daemen, L. L. Superionic Conductivity in Lithium-Rich Anti-Perovskites. J Am Chem Soc 134, 15042–15047 (2012).
Bhattacharya, J. & van der Ven, A. First-principles study of competing mechanisms of nondilute Li diffusion in spinel Li_xTiS₂. Physical Review B 83, 144302 (2011).
Hare, S. R. & Tantillo, D. J. Post-transition state bifurcations gain momentum – current state of the field. Pure and Applied Chemistry 89, 679–698 (2017).
Knuth, D. E. A generalization of Dijkstra’s algorithm. Information Processing Letters 6, 1–5 (1977).
Dijkstra, E. W. A note on two problems in connexion with graphs. Numerische Mathematik 1, 269–271 (1959).
Otero-de-la-Roza, A., Johnson, E. R. & Luaña, V. Critic2: A program for real-space analysis of quantum chemical interactions in solids. Computer Physics Communications 185, 1007–1018 (2014).
Momma, K. & Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography 44, 1272–1276 (2011).
Seiler, P. & Dunitz, J. D. Are ionic solids really built of ions? New evidence from x-ray diffraction. Helvetica Chimica Acta 69, 1107–1112 (1986).
Downs, R. T., Gibbs, G. V., Boisen Jr, M. B. & Rosso, K. M. A comparison of procrystal and ab initio model representations of the electron-density distributions of minerals. Physics and Chemistry of Minerals 29, 369–385 (2002).
Blatov, V. A. Voronoi–dirichlet polyhedra in crystal chemistry: theory and applications. Crystallography Reviews 10, 249–318 (2004).
Chung, S. J., Hahn, Th. & Klee, W. E. Nomenclature and generation of three-periodic nets: the vector method. Acta Crystallographica Section A 40, 42–50 (1984).
Delgado-Friedrichs, O. & O’Keeffe, M. Identification of and symmetry computation for crystal nets. Acta Crystallographica Section A 59, 351–360 (2003).
Andersen, M., Panosetti, C. & Reuter, K. A Practical Guide to Surface Kinetic Monte Carlo Simulations. Frontiers in Chemistry 7, (2019).
Van der Ven, A., Ceder, G., Asta, M. & Tepesch, P. D. First-principles theory of ionic diffusion with nondilute carriers. Physical Review B 64, 184307 (2001).
Bortz, A. B., Kalos, M. H. & Lebowitz, J. L. A new algorithm for Monte Carlo simulation of Ising spin systems. Journal of Computational Physics 17, 10–18 (1975).
Bulnes, F. M., Pereyra, V. D. & Riccardo, J. L. Collective surface diffusion: η-fold way kinetic Monte Carlo simulation. Physical Review E 58, 86–92 (1998).
Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science 6, 15–50 (1996).
Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals. Physical Review B 47, 558–561 (1993).
Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B 54, 11169–11186 (1996).
Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Physical Review Letters 77, 3865–3868 (1996).
Blöchl, P. E. Projector augmented-wave method. Physical Review B 50, 17953–17979 (1994).
Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B 59, 1758–1775 (1999).

(Not answered)

SI.pdf
Supplementary information

Download PDF

Journal Publication

published 30 Aug, 2022

Read the published version in npj Computational Materials →

Editorial decision: revise
20 Jul, 2022
Review #2 received at journal
19 Jul, 2022
Review #3 received at journal
05 Jul, 2022
Reviewer #3 agreed at journal
05 Jul, 2022
Reviewer #2 agreed at journal
02 Jul, 2022
Review #1 received at journal
28 Jun, 2022
Reviewer #1 agreed at journal
22 Jun, 2022
Reviewers invited by journal
18 Jun, 2022
Submission checks completed at journal
10 Jun, 2022
Editor assigned by journal
09 Jun, 2022
First submitted to journal
09 Jun, 2022

You are reading this latest preprint version

Enhancing first-principles simulations of complex solid-state ion conductors using topological analysis of procrystal electron density

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results And Discussions

Progress beyond the state of the art

Grain boundary Li-ion transport in polycrystalline Li₃ClO solid electrolytes

Methods

Calculation of ionic diffusion trajectories and migration barriers using TAPED

Choice of the initial and final positions for calculating the diffusion trajectory

Analysis of ion migration map

Diffusion coefficient and KMC simulations

DFT calculations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1

Enhancing first-principles simulations of complex solid-state ion conductors using topological analysis of procrystal electron density

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results And Discussions

Progress beyond the state of the art

Grain boundary Li-ion transport in polycrystalline Li3ClO solid electrolytes

Methods

Calculation of ionic diffusion trajectories and migration barriers using TAPED

Choice of the initial and final positions for calculating the diffusion trajectory

Analysis of ion migration map

Diffusion coefficient and KMC simulations

DFT calculations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1

Grain boundary Li-ion transport in polycrystalline Li₃ClO solid electrolytes