Unsupervised landmark analysis for jump detection in molecular dynamics simulations

TL;DR

This paper introduces an unsupervised method for detecting ion jumps in molecular dynamics simulations, enabling analysis without prior knowledge of the diffusive process, and applies it to various lithium-ion conductors.

Contribution

The novel approach combines landmark-based projection, clustering, and atomic environment descriptors to identify and analyze jump events in ionic diffusion simulations.

Findings

Successfully applied to ten lithium-ion systems

Revealed detailed diffusive behaviors in specific materials

Compared favorably with existing methods

Abstract

Molecular dynamics is a versatile and powerful method to study diffusion in solid-state ionic conductors, requiring minimal prior knowledge of equilibrium or transition states of the system's free energy surface. However, the analysis of trajectories for relevant but rare events, such as a jump of the diffusing mobile ion, is still rather cumbersome, requiring prior knowledge of the diffusive process in order to get meaningful results. In this work, we present a novel approach to detect the relevant events in a diffusive system without assuming prior information regarding the underlying process. We start from a projection of the atomic coordinates into a landmark basis to identify the dominant features in a mobile ion's environment. Subsequent clustering in landmark space enables a discretization of any trajectory into a sequence of distinct states. As a final step, the use of the smooth overlap of atomic positions descriptor allows distinguishing between different environments in a straightforward way. We apply this algorithm to ten Li-ionic systems and conduct in-depth analyses of cubic Li$_{7}$La$_{3}$Zr$_{2}$O$_{12}$, tetragonal Li$_{10}$GeP$_{2}$S$_{12}$, and the $\beta$-eucryptite LiAlSiO$_{4}$. We compare our results to existing methods, underscoring strong points, weaknesses, and insights into the diffusive behavior of the ionic conduction in the materials investigated.

Tables

Table 1. Table 1: Comparison between the presented landmark analysis and the sites listed in CIF files taken from structural databases. A checkmark indicates structures where the sites in the CIF file and the results of the unsupervised analysis agree, both in number of site types and number of sites of each type. This is the case in all but two structures, \ce Li20Re4N16 and \ce Li24Zn4O16, which are discussed in the text.

Structure	Sites	Sites
	Wyckoff (CIF)	(Landmark + SOAP)
\ceLi32Al16B16O64	16e + 16e	16 + 16 ✓
\ceLi24Sc8B16O48	N/A	8 + 16
\ceLi24Ba16Ta8N32	8e + 16f	8 + 16 ✓
\ceLi20Re4N16	4a + 16g	8 + 16
\ceLi12Rb8B4P16O56	4d + 8g	4 + 8 ✓
\ceLi6Zn6As6O24	3b + 3b	3 + 3 ✓
\ceLi24Zn4O16	16f + 8d	20 + 8

Table 2. Table 2: For every structure analyzed in Sec.  III.4 , we list the simulation length T s ​ i ​ m subscript 𝑇 𝑠 𝑖 𝑚 T_{sim} in picoseconds, the database (DB) from which the structure was retrieved, and the structure’s database ID.

Structure	$T_{sim}$ (ps)	DB	DB-ID
\ceLi32Al16B16O64	72	ICSD	50612
\ceLi24Sc8B16O48	218	COD	2218562
\ceLi24Ba16Ta8N32	58	ICSD	75031
\ceLi20Re4N16	159	ICSD	92468
\ceLi12Rb8B4P16O56	116	ICSD	424352
\ceLi6Zn6As6O24	226	ICSD	86184
\ceLi24Zn4O16	407	ICSD	62137

Equations

$$R_i=\{x\in {\mathbb{R}}^n:|\mathbf{x}-{\mathbf{s}}_i|\le |\mathbf{x}-{\mathbf{s}}_j|\text{ for all }j\ne i\}$$

$$d_{m,h}^{\mathrm{A}}(t)=\frac{\left|{\mathbf{r}}_m(t)-{\mathbf{r}}_h(t)\right|}{\left|{\bar{\mathbf{r}}}_{\text{VN}}^{\mathrm{A}}-{\bar{\mathbf{r}}}_h\right|},h\in \mathrm{A},$$

$$l_{\mathrm{A}}^m=\prod _h^{N_{\text{coord}}}{\left[f\left(d_{m,h}^{\mathrm{A}}\right)\right]}^{\frac{1}{N_{\text{coord}}}},$$

$$\sigma (d;d_0,k)=\frac{1}{1+e^{-k\left(d-d_0\right)}}.$$

$$f(d;d_0,k)=$$

$$S({\mathbf{l}}_A,{\mathbf{l}}_B)=\frac{{\mathbf{l}}_A\cdot {\mathbf{l}}_B}{|{\mathbf{l}}_A||{\mathbf{l}}_B|},$$

$$[m_{AB}]=\{\begin{array}{ll}0& \text{if }|{\mathbf{r}}_A-{\mathbf{r}}_B|>\text{cutoff}\\ p_{A\to B}& \text{otherwise}\end{array}$$

$$n_M(\mathbf{r})={⟨\sum _m^M\delta (\mathbf{r}-{\mathbf{r}}_m(t))⟩}_t,$$

$$D_{\mathrm{tr}}^M=\underset{\tau \to \infty }{\lim }\frac{1}{6\tau }\frac{1}{N_M}\sum _m^M{⟨|{\mathbf{r}}_m(\tau +t)-{\mathbf{r}}_m(t){|}^2⟩}_t,$$

$$g_{M-S}(r)=$$

$$=$$

$$\sigma =\frac{Z^2{e}^{2}N}{k_{B}T}\frac{D_{tr}}{H},$$

$$j=\underset{j\in [0,N^{\prime })}{\mathrm{argmax}}S({\mathbf{c}}_i,{\mathbf{c}}_j^{\prime })$$

$$S({\mathbf{c}}_i,{\mathbf{c}}_j^{\prime })=\frac{{\mathbf{c}}_i\cdot {\mathbf{c}}_j^{\prime }}{|{\mathbf{c}}_i||{\mathbf{c}}_j^{\prime }|}$$

$$S({\mathbf{c}}_i,{\mathbf{c}}_j^{\prime })>\text{clustering threshold}$$

$${\mathbf{c}}_j^{\prime }=\frac{n{\mathbf{c}}_j^{\prime }+{\mathbf{c}}_i}{n+1}$$

$${\mathbf{c}}_{N^{\prime }}^{\prime }$$

$$N^{\prime }$$

$$s=\underset{j\in [0,N)}{\max }S(\mathbf{l},{\mathbf{c}}_j)$$

Full text

I Introduction

Lithium-ion batteries power an increasingly broad and critical set of technologies Armand and Tarascon (2008). Commercially available batteries use organic electrolytes that impose constraints on their safety, power and energy density Schaefer et al. (2012) and can introduce chemical instabilities that require the incorporation of fuses and safety vents Balakrishnan et al. (2006). Solid-state electrolytes are widely considered to be a promising alternative for next-generation batteries: many structural families of candidate solid-state ionic conductors have been identified and are under investigation Bachman et al. (2016). While a good solid-state electrolyte must meet several criteria, such as low electronic mobility, easy device integration, and electrochemical and mechanical stability Manthiram et al. (2017), it must first be a fast Li-ion conductor, and consequently optimization of conductivity and analysis of the mechanisms of Li-ion diffusion has been the focus of a large body of literature Bachman et al. (2016).

Atomistic modeling techniques, and in particular molecular dynamics (MD), have been used to study a wide variety of candidates for solid-state electrolytes and the factors that influence their ionic conductivity. Classical/empirical force fields were chosen in several studies Adams and Rao (2012a, b); Xu et al. (2012a); Deng et al. (2015); Kozinsky et al. (2016); Klenk and Lai (2016); Burbano et al. (2016); Dawson et al. (2018) due to their computational efficiency and access to the time and length scales required to characterize ionic transport. Accurate, yet expensive, first-principles simulations have also been employed for selected systems Wood and Marzari (2006); Ong et al. (2013); Mo et al. (2012); Xu et al. (2012b); Mo et al. (2014); Meier et al. (2014); Wang et al. (2015); Chu et al. (2016); Zhu et al. (2017); Marcolongo and Marzari (2017); Sagotra et al. (2019). The necessary compromise between the transferability of first-principles potential energy surfaces and the computational efficiency of force fields has also motivated the development of novel hybrid quantum/classical approaches Kahle et al. (2018) to model diffusion. The estimate of transport coefficients from the Green-Kubo or Einstein relations using molecular dynamics can be done in a straightforward yet expensive way, though improved methods for obtaining accurate estimates from short trajectories are being developed Ercole et al. (2017). In addition to computing ionic conductivity, design and characterization of new materials requires detailed understanding of the atomistic mechanisms of ionic transport. The central challenge is to develop automated methods for accurately analyzing the structure and dynamics of lithium’s local atomic environments and for detecting rare transitions and subtle correlations in large amounts of data.

In many solid-state Li-ion conductors, Li ions form a mobile, often disordered, sublattice within a non–diffusive sublattice of the other species, which we refer to as the host lattice hereafter. In the jump-diffusion model, the mobile ions spend the majority of their time in the local minima of the potential energy surface and vibrate within such sites for a sufficiently long time to lose memory of their previous locations while intermittently acquiring sufficient kinetic energy to overcome the barrier separating them from a different potential well. This formulation of Li-ion diffusion as occupation of and exchange between well-defined crystallographic sites can be used to model diffusion as a Markov-chain model using kinetic Monte Carlo Van der Ven et al. (2001). Also, using this discrete formulation to understand the microscopic origin of diffusion is a common theme in the literature, and site analysis tools have been used to explore the effects of site volume Wang et al. (2015); Kweon et al. (2017) and anion sublattice structure Wang et al. (2015) on ionic conductivity, to identify conduction pathways and rate-limiting steps Kozinsky et al. (2016); de Klerk et al. (2018), and correlated diffusion events Chen et al. (2017), and to design new descriptors for conductivity Kweon et al. (2017). Ideally, an automated site analysis (see Fig. 1) approach should: (1) automatically identify relevant Li sites, (2) accurately track migration of mobile ions through those sites, (3) require no prior knowledge of the material, and (4) work with the same parameters over a broad range of materials.

Existing approaches fall mostly into three classes: distance-based, topology-based, and density-based methods. Distance-based methods Varley et al. (2017); Kweon et al. (2017); de Klerk et al. (2018) use preexisting knowledge of the equilibrium positions of all Li-ion sites and consider a Li ion to be resident at a site when it comes within a given cutoff distance from the site’s position. Cutoffs can be smooth Kweon et al. (2017) or discrete de Klerk et al. (2018), but in both cases, they need to be tailored to the structure at hand and are uniform across all sites within it. The positions of sites can also be coupled to the instantaneous positions of nearby host-lattice atoms Kweon et al. (2017) to decrease sensitivity to thermal noise. Nevertheless, such methods rely on the crystallographic information they are given and also do not account for the varied or non-spherical geometry of sites Deng et al. (2018). Starting from a prior knowledge of the host structure and possible Li sites, mobile ions can also be automatically assigned to sites based on convex-hull analysis of site polyhedra Kozinsky et al. (2016); Kozinsky (2018). This topology-based method deals with arbitrary site geometries, eliminates thermal noise and does not require arbitrary distance cutoffs, but does require the site polyhedra to be specified. Density-based methods Chen et al. (2017) identify regions of high Li-ion density separated by areas of low Li-ion density, as determined by a threshold, and define each high-density region as a site. These methods thus do not require prior knowledge about the material and can resolve sites with different geometries. In materials with nearby or rapidly exchanging sites, however, choosing a density threshold that can distinguish such sites from one another can be difficult. Richards et al. Richards et al. (2016) used a k-means clustering of Na-ion positions in \ceNa10GeP2S12, initialized with known ionic positions for the similar ionic conductor \ceLi10GeP2S12, combining prior information with a density based method.

In order to overcome some of these challenges, this work introduces an algorithm for accurately and automatically analyzing molecular dynamics trajectories and detecting jumps of the mobile ion through the host lattice with minimal human supervision and no prior knowledge. This algorithm can be combined with the automatic detection of important structural motifs Gasparotto and Ceriotti (2014); Gasparotto et al. (2018), leading to a versatile tool for the unsupervised analysis of trajectories and detection of diffusion events. The algorithm will be discussed in Sec. II; in Sec. III we apply it to three known ionic conductors and to seven non-diffusing materials and discuss the results; some details of the implementation are given in Sec. IV; and our final conclusions are presented in Sec. V.

II Algorithm

Landmarks are persistent local features in an environment and therefore can be used to describe positions in the absence of global information (i.e. real-space coordinates). Landmark-based navigation explains the homing of social insects Wehner et al. (1996) and has been applied in the field of autonomous navigation and artificial intelligence Möller (2000). Landmark models employ a vector-based description of the environment via a landmark vector $\bm{l}$. Such a vector representation is useful for navigation if the distances between the landmark vectors corresponding to two states or positions $A$ and $B$ decrease with reducing distances in real space: $|\bm{r}_{A}-\bm{r}_{B}|=f(|\bm{l}_{A}-\bm{l}_{B}|$), where $f$ is a monotonically increasing function of its input.

When analyzing trajectories, the real-space positions are obviously known beforehand. However, atomic coordinates are inefficient descriptors for most properties since they are not invariant under rigid translation or rotation of the structure. We will describe the positions of mobile ions through landmarks that encode all the information necessary to detect changes in the ions’ environments and are invariant under these transformations.

First, we deduce that the descriptors should only encode local information since the local environment mostly defines the potential energy landscape for the mobile ion, a principle reminiscent of the nearsightedness of electronic matter Prodan and Kohn (2005). In addition, we know from Pauling’s rules in crystal structures Pauling (1929) that ionic systems minimize their energy by packing into coordination numbers that are determined, among other factors, by the ratio of the radii of the cations and anions. Therefore, possible coordination polyhedra in the local environment are meaningful features. Checking all the possible polyhedra in a crystal is not feasible because of the combinatorial complexity this induces, so we need to restrict the description via a meaningful subset of convex hulls or polyhedra formed by the host lattice. A site description and trajectory discretization via pre-selected convex hulls has been previously developed and applied Kozinsky et al. (2016) to study Li-ion diffusion in garnets.

Using polyhedra defined by host-lattice atoms as landmarks relies on the assumption that these atoms fluctuate around equilibrium positions, such that well-defined coordination polyhedra persist throughout the simulation. Equivalently, the host lattice is not changing in a way that causes sites to appear or disappear. Due to this assumption, the present landmark analysis cannot be applied to systems with liquid-like host structure such as polymers, where inter-site hopping happens on a longer time scale than the host motion, and dynamic coordination tracking must be used Molinari et al. (2018). It also cannot be applied to systems with a “paddle-wheel” diffusion mechanism, where the slow rotation of polyatomic anions creates a constantly changing set of local potential minima, such as shown for proton diffusion in \ceCsHSO4 Wood and Marzari (2007) or lithium- and sodium-ionic diffusion in the closoborate structures Kweon et al. (2017). Similar to the site analysis presented in the literature, our method does not assume that the occupation a given site are Markovian, i.e. whether a mobile ion completely loses memory of its past at any site: We define and find a site based on stable and persistent features in the environment of mobile ions, described by the landmark vectors, without considering information in the time domain. Whether the underlying process is Markovian can be determined by analyzing the resulting statistics Chen et al. (2017); Morgan and Madden (2014).

The basic algorithm has three steps: (1) definition of suitable landmarks, (2) expression of the coordinates of the mobile ions during their trajectory in the landmark basis, and (3) clustering of the landmark vectors to reveal sites and discretize the trajectory of each mobile atom. We also implemented, as an option, the possibility to: (4) merge nearby sites that have high exchange rates and that fulfill some distance criteria and (5) determine site types based on the geometry and chemistry of the local environment.

While the two last steps are optional and independent from each other, we always apply them in the analysis that we show in Sec. III. Step 4 reduces significantly the noise in the data, and step 5 supplies information on the local geometry and chemistry. In the following we explain each of the steps in greater detail and finish with a discussion on why certain design choices were taken.

II.1 Step 1: Define Landmarks

The landmark analysis we introduce here is based on the Voronoi tessellation of the equilibrium configuration of the host lattice and its geometric dual, the Delaunay triangulation. Given a set of points in space, termed seeds, a Voronoi tessellation divides space into regions such that all points in each region are closer to the region’s seed than to any other seed Okabe et al. (2009). Formally, the Voronoi region determined by the seed point $s_{i}\in\mathbb{R}^{n}$ is given by:

[TABLE]

Voronoi regions connect at Voronoi facets, as shown in the schematic in Fig. 2. Any point on such a facet is equidistant to the seeds of the adjacent Voronoi regions. Voronoi nodes are, in a space of $D$ dimensions, points where at least $D$ facets intersect and therefore are equidistant to at least $D+1$ seed points. It follows that each Voronoi node locally maximizes the distance to its adjacent seed points. The geometric dual of the Voronoi tessellation is the Delaunay triangulation. Such a triangulation or simplicial 111A simplex in $\mathbb{R}^{D}$ is the convex hull of $D+1$ points that do not lie on a hyperplane. decomposition is obtained by connecting seed points that share a Voronoi facet. The Delaunay triangulation has the useful property that the circumcircles of all formed triangles have empty interiors, i.e. there are no seed points inside any circumcircle. It follows from the duality between the two tessellations that every Voronoi node is associated to exactly one Delaunay simplex. The Voronoi node lies at the center of the circumcircle of the associated Delaunay simplex. In the remainder of the text we will work in three dimensions unless otherwise specified; in three dimensions a Voronoi node is equidistant to at least four coordinating seeds.

While a Voronoi node is a reasonable guess for a low-energy position since it maximizes the distance to its coordinating seeds, the associated Delaunay simplex corresponds to the coordinating polyhedron of the site or a subset thereof. A Voronoi node and its coordinating host-lattice atoms are together referred to as a landmark. The coordinating host-lattice atoms of a landmark are the host-lattice atoms that are vertices of the Delaunay simplex – dual to the Voronoi node – in the equilibrium configuration.

II.2 Step 2: Landmark Vectors

We start from a molecular dynamics trajectory that gives the real-space positions $\mathbf{r}_{i}(t)\in\mathbb{R}^{3}$ of each atom $i$ at time $t=N\Delta t_{s}$, an integer multiple of the sampling timestep $\Delta t_{s}$. In the remainder, we will use the index $h$ for host-lattice atoms and $m$ for mobile ions. First, we calculate the time-averaged positions for host-lattice atoms $\bar{\mathbf{r}}_{h}=\langle\mathbf{r}_{h}\rangle_{t}$ and use these as seed points for a Voronoi decomposition, resulting in Voronoi nodes $\bar{\mathbf{r}}_{\text{VN}}^{\mathrm{A}}$. The instantaneous position of a mobile particle, $\mathbf{r}_{m}(t)$, is expressed in terms of a proximity or similarity to each landmark in the system. That is to say, any real space position $\mathbf{r}_{m}$ can be transformed into a vector in the $N$-dimensional landmark space, where $N$ is the number of landmarks in the system, equal to the number of Voronoi nodes and also equal to the number of Delaunay simplices, due to the duality discussed in Sec. II.1. We index landmarks with capital latin characters. For a landmark $\mathrm{A}$, we first define the normalized instantaneous distance between a mobile particle $m$ and a host lattice atom $h$, where atom $h$ is one of the coordinating seed atoms of the landmark $\mathrm{A}$:

[TABLE]

where $\mathbf{r}_{h}(t)$ and $\mathbf{r}_{m}(t)$ are the instantaneous real-space positions of host-lattice atom $h$ and mobile ion $m$, respectively, $\bar{\mathbf{r}}_{h}$ is the time-averaged position of host-lattice atom $h$, and $\bar{\mathbf{r}}_{\text{VN}}^{\mathrm{A}}$ is the position of the landmark’s Voronoi node. The corresponding component of the landmark vector is then computed as:

[TABLE]

where $h$ ranges over the set of $N_{\text{coord}}$ coordinating host-lattice atoms and $f(d)$ is a cutoff function that smoothly goes from 1 to 0. We base the cutoff function on the logistic function $\sigma(d;d_{0},k)$, a sigmoid curve that goes from 0 to 1 around a midpoint $d_{0}$, with a steepness $k$:

[TABLE]

To obtain a cutoff function suitable for Eq. (2), we subtract the logistic function from 1:

[TABLE]

The function $f(d;d_{0},k)$ varies smoothly from 1 to 0, reaching $\frac{1}{2}$ at the set midpoint $d_{0}$. How fast it varies is tuned by the hyperparameter $k$. As can be seen from Eq. (2), we normalize $f\left(d_{m,h}^{\mathrm{A}}\right)$ for a varying $N_{\text{coord}}$. In three dimensions and in the present framework, $N_{\text{coord}}$ is always $4$, since we use a simplicial decomposition to determine the landmarks. However, the framework could be changed to include a varying number of coordinating host-lattice atoms, motivating this normalization. The cutoff function in Eq. (4) was preferred due to its continuity and simplicity. Because distances are normalized to the equilibrium distance between the Voronoi node and the host atoms, the magnitude of each landmark vector component depends on neither the volume nor shape of the corresponding landmark’s polyhedron. This allows landmark analysis to distinguish between sites whose coordination polyhedra have very different volumes, as well as accurately tracking mobile particles through highly distorted sites.

II.3 Step 3: Landmark Clustering

The magnitude of each component of the landmark vector indicates the extent to which a mobile atom’s position is dominated by that landmark. If, for example, a mobile atom occupies a tetrahedral site, its landmark vector would have one large value at the corresponding landmark’s component and some low-magnitude noise for neighboring landmarks. During a transition between sites, there are no dominant contributions, as shown schematically in two dimensions in Fig. 3.

If an atom occupies an octahedral site, however, the landmark vector will have four major contributions, corresponding to the four tetrahedrons resulting from the Delaunay triangulation of the octahedron. We show this schematically for two dimensions in Fig. 4.

Because we have chosen a smooth function of position for the landmark vector components, the landmark vectors are a continuous function of trajectory time. By definition, landmark vectors are invariant under rigid translations or rotations of the system and as such are ideally suited as descriptors for dominant recurring features. A clustering of the landmark vectors can be used to group similar landmark vectors and therefore discretize our trajectory in landmark space. We use density-based clusters of landmark vectors, where each cluster is described by a high-density region in landmark space, corresponding to a frequent feature in the local environment of the mobile ion. Therefore, we define sites as clusters in landmark space.

We use a custom hierarchical clustering algorithm (described in more detail in Appendix A) with a simple cosine similarity metric:

[TABLE]

where $\mathbf{l}_{A/B}$ are landmark vectors. The clustering algorithm scales linearly with the number of landmark vectors.

The clustering algorithm is run on the landmark vectors computed from the real-space positions of all mobile atoms every $n$ frames, where $n$ is sufficiently small and corresponds to a time span that is below the jump rate. A mobile atom is said to be occupying site $i$ at time $t$ if its corresponding landmark vector at that time is a member of the $i$-th landmark cluster. If the mobile atom’s landmark vector is not a member of any cluster, the atom is said to be unassigned at that time. The time sequence of such site assignments for a given mobile atom is its discretized trajectory; every change of site in that discretized trajectory is defined as a jump event.

The center of each site is defined as the spatial average of all real-space positions of mobile ions assigned to it.

II.4 Step 4 (optional): Merge Sites

While one of the main strengths of the landmark analysis is its ability to distinguish between very close sites, that level of resolution often identifies multiple sites where only one should exist. This is mainly due to a lack of data for the clustering. This issue is particularly prominent in host lattices containing sites with greater than four-fold coordination whose coordination polyhedra are highly distorted from the corresponding regular polyhedra. To merge such split sites, a post-processing clustering of the sites themselves can be applied, taking into account information from the time domain. We define $M$ as the stochastic matrix observed from the exchanges of ions between sites:

[TABLE]

where $\mathbf{r}_{A}$ is the center of site $A$ and $p_{A\rightarrow B}$ is the probability that an ion occupies site $B$, conditional on the ion’s having occupied site $A$ in the previous frame (for $A\neq B$). For $A=B$ it is the probability that an ion remains at site $A$ until the next frame. We apply Markov Clustering Van Dongen (2008) to the weighted graph defined by the stochastic matrix $M$, resulting in clusters of highly-connected subgraphs. Sites belonging to the same subgraph are merged (additional details are given in Appendix B).

II.5 Step 5 (optional): Site Type Analysis

Sites are commonly defined by their Wyckoff points, and symmetry-equivalent sites can be interpreted as one site type. Such analysis depends on preexisting crystallographic data and also neglects that the energetics of a site are defined by the local geometry and chemistry. In line with our goal of making unsupervised site analysis possible, we developed a method for determining the type of the sites identified by the steps described from Sec. II.1 to Sec. II.4. Different sites whose environments cannot be distinguished are said to be of the same site type.

We describe local atomic environments using the smooth overlap of atomic positions (SOAP) Bartók et al. (2013) as implemented in the QUIP molecular dynamics framework noa (2018). Briefly, a SOAP descriptor is a vector that describes the local geometry around a point in a rotation-, translation-, and permutation-invariant way. The descriptor changes smoothly with the Cartesian coordinates of the structure. For these reasons, SOAP descriptors have become a powerful tool to express local geometry for machine-learning applications De et al. (2016) and the detection of structural motifs Gasparotto and Ceriotti (2014); Gasparotto et al. (2018).

Multiple SOAP vectors must be computed for each site to provide sufficient data density for subsequent clustering. Computing these vectors for a site requires some procedure for sampling the real-space positions of both the site and its surrounding host-lattice atoms. We implemented two sampling schemes. (1) Real-space averaging: the real-space positions of all mobile atoms when they occupy the site are collected, and $n$ average real-space positions are computed for the site, where $n$ is a parameter chosen by the user. SOAP is computed on the averaged sites. (2) SOAP-space averaging: SOAP vectors are computed for all real-space positions with the host-lattice atoms at their corresponding instantaneous positions. Then, $n$ average descriptor vectors are computed in SOAP space.

After reducing the dimensionality of the SOAP vectors with Principal Component Analysis, we cluster them using density-peak clustering Rodriguez and Laio (2014) with a Euclidean distance metric. A simple parameter estimation scheme is used to determine the number of clusters (see Appendix C). Each cluster of descriptor vectors corresponds to a site type. Each site is assigned to the type corresponding to the descriptor cluster to which the majority of its descriptors were assigned. Small majorities (less than 70-80% agreement) typically indicate insufficient data, poorly chosen SOAP parameters, or very similar environments.

II.6 Discussion of design choices

The main motivation for a landmark based approach is its ability to significantly reduce noise resulting from thermal vibrations in the system while reducing the dimensionality and discretizing the trajectory of the mobile ions. The design described in this section is driven by physical intuition and trial-and-error. While developing the present approach, we attempted and discarded a number of approaches due to poor performance in trial systems. (1) Directly clustering the Cartesian coordinates of the mobile ions (the density-based approach discussed in Sec. I) was found to work poorly in some systems. We show this in more detail in Sec. III.1. (2) An analysis based on the $N$ nearest neighbors of the mobile atom was tried but discarded, since we could not determine $N$ without relying on the knowledge of the structure under investigation, in particular the expected size of the coordination shell of Li. (3) We tried various landmark representations, the most simple being the distance to each host-lattice atom, therefore taking the instantaneous positions of host-lattice atoms as landmarks. The results for different systems were not satisfactory.

We conclude this section by pointing out that passing from Cartesian coordinates to landmark vectors can significantly reduce the noise that comes mostly from thermal vibrations in the system. Different formulations of landmarks can be envisioned, and while the present framework performs well for Li-ionic diffusion, a different landmark framework might be needed to describe, for example, Grotthus-like proton diffusion in superprotonic \ceCsHSO4 Wood and Marzari (2007).

III Results and discussion

We apply the algorithms above to ten representative materials, \ceLi7La3Z2O12, \ceLiAlSiO4, \ceLi10GeP2S12, \ceLi32Al16B16O64, \ceLi24Sc8B16O48, \ceLi24Ba16Ta8N32, \ceLi20Re4N16, \ceLi12Rb8B4P16O56, \ceLi6Zn6As6O24 and \ceLi24Zn4O16. For the subsequent analysis, we also calculate radial distribution functions, mean-square displacements, and ionic densities. The mobile-ion densities $n_{\text{M}}({\bm{r}})$ are calculated from molecular dynamics trajectories as:

[TABLE]

where the index $m$ runs over all mobile ions $M$ and the angular brackets $\langle\cdot\rangle_{t}$ indicate a time average over the trajectory, which is equal to an ensemble average under the assumption of ergodicity. When applying Eq. (6) we replace the delta function by a Gaussian with a standard deviation of 0.3Å, and the summation is performed on a grid of ten points per Å in every direction. The tracer diffusion coefficient of the mobile species $D_{\mathrm{tr}}^{M}$ is computed from the mean-square displacement of the mobile ions as a function of time:

[TABLE]

where $\bm{r}_{m}(t)$ is the position of the mobile ion at time $t$. In practice, we fit a line to the mean-square displacement in the diffusive regime. The error of the tracer diffusion coefficients is estimated with a block analysis Allen and Tildesley (1987). The radial distribution function $g(r)_{M-S}$ of the mobile ions $M$ with species $S$ is calculated as:

[TABLE]

where $f(r)$ is the ideal-gas average number density at the same overall density. In addition, we integrate the average number density $\rho(r)$ to give the average coordination number as a function of distance Frenkel and Smit (1996).

III.1 Analysis of Li7La3Zr2O12

Garnet-type structures were proposed as lithium-ionic conductors by Thangadurai et al. Thangadurai et al. (2003). The general formula of garnets is \ceLi5La3M2O12 $\left(\mathrm{M=Ta,Nb}\right)$ Knauth (2009), but aliovalent substitutions of $M$ can change the lithium content. Xie et al. Xie et al. (2011) studied in more detail the distribution of Li+ in garnets. Their results indicate that increasing the lithium concentration in garnets leads to an increase in occupation of octahedral sites, which is confirmed by simulations Kozinsky et al. (2016) and also in experiments O’Callaghan et al. (2008). It has been established O’Callaghan and Cussen (2007); Adams and Rao (2012b); Thangadurai et al. (2014) for the garnet structure that Li ions can occupy tetragonal 24d sites, octahedral 48g sites and 96h distorted octahedral sites. The latter stem from a site splitting of the 48g sites to increase the Li-Li distances and occur at higher lithium concentrations. In this work, we study the Li-ion distribution of Zr-based cubic garnets with the stochiometric formula \ceLi7La3Zr2O12, referred to as LLZO in the remainder.

We sample the dynamics in the cell of 192 atoms in the canonical ensemble via a GLE thermostat Ceriotti et al. (2009) at a temperature of 500 K, using a lattice constant of 12.9872 Å, and a polarizable force field. We use LAAMPS Plimpton (1995) to perform the simulation for 10 ns, with the parameters of the force field taken from the work by Mottet et al. Mottet et al. (2019), which accurately reproduces the kinetics of the diffusing process in LLZO.

The estimate of the diffusion coefficient via Eq. (7) reveal that Li-ions indeed are diffusive in LLZO, with a tracer diffusion coefficient of $D_{tr}^{Li}=2.4\times\mathrm{10^{-6}\,cm^{2}\,s^{-1}}$. Application of the Nernst-Einstein equation gives the ionic conductivity $\sigma$:

[TABLE]

where $\left(Ze\right)$ the carrier’s charge, $N$ the carrier density, $k_{B}$ the Boltzmann constant, $T$ the absolute temperature and $H$ is the ratio between the tracer and charge diffusion coefficient, commonly referred to as Haven ratio: $H=\frac{D_{tr}}{D_{\sigma}}$. To account for the strong evidence for correlated motion in this material Jalem et al. (2013); Meier et al. (2014), we set the Haven ratio to $H=0.4$, reported in a study Morgan Benjamin J. (2017) for this Li-ion concentration in LLZO. We find $\mathrm{\sigma=0.58\,Scm^{-1}}$, which is one order of magnitude larger than the values reported by Murugan et al. Murugan et al. (2007) This is within the acceptable range, especially for a classical force-field, and not of concern since the focus of this work is the analysis method. The diffusive pathways can be illustrated by the Li-ion density, shown for three isosurfaces in Fig. 5. By visual inspection, the densities look similar to those presented by Adams and Rao in their computational study Adams and Rao (2012b). The splitting of 48g sites into 96h sites Thangadurai et al. (2014) cannot be seen from the isosurfaces at any isovalues, which is consistent with the conclusions drawn by Chen et al. Chen et al. (2017) that density-based clustering of real-space positions cannot resolve the two distinct 96h sites in LLZO from the 48g site. The radial distribution function $g(r)$, shown in the upper panel of Fig. 6, shows how the Li ions in our simulation are, as expected, coordinated closest by oxygen and then by other Li ions.

We use the site analysis presented to discretize the trajectory of lithium ions into meaningful states, as illustrated for one lithium ion in Fig. 7. The subsequent SOAP analysis produces two clearly resolvable clusters, which are detected by the clustering algorithm. We show the first two principal components in Fig. 8, with a color encoding representing the cluster detected. It is evident that the SOAP descriptor produces data that clusters well in this projection and that the clustering algorithm correctly assigns the clusters. The algorithm detects 24 sites of one kind (type 1) and 83 of another (type 0). We attribute the tetrahedral environment to the former, and the octahedral environment to the latter, since the expected values are 24 sites for the tetrahedral environment and 96 for the octahedral one. The last number is due to the site splitting inside each of the 48 octahedral cavities, leading to two sites inside each octahedral cavity. The under-prediction of the number of octahedral sites is due to the merging, in some cases, of octahedral sites into a single site. We stress that the numbers of sites presented as final results are after the site-merging step, presented in Sec. II.4. The proximity and fast ion exchange between the octahedral sites in the same cavity explains why our algorithm does not give the correct answer, but it is remarkably close to the correct result, without any encoding of prior information about possible site splitting. Comparing to the study by Chen et al. Chen et al. (2017), we can conclude that the landmark analysis is able to better distinguish minima in close proximity. We speculate that the main reason is a higher tolerance for thermal vibrations of the host lattice, that can lead to energetic minima being spread in real space.

To ensure that the analysis of LLZO provides reasonable and expected results, we calculate the Li-O radial distribution function $g(r)$ separately for each site type; these are shown in the bottom panel of Fig. 6. Li ions attributed to sites of type 1 have an environment characterized by a distinct nearest-neighbor peak stemming from four-fold coordination of lithium with oxygen, as evidenced by the integral plateauing at a value of four. The first peak for type 0, shown in red, has a shape compatible with a distorted octahedron, due to the appearance of a shoulder, and the weak, but distinguishable, plateau of the integral at a value of approximately 6. This is further evidence that the site types have been correctly attributed to the tetrahedral and octahedral site environments of LLZO, and that the SOAP descriptor can be used to cluster site types correctly. Additionally, we resolve the Li-ion density by site type in Fig. 9. The isosurfaces are compatible, by visual inspection, with reported work Adams and Rao (2012b). We also calculate the jump lag, which is the average residence time at a site A before jumping to a site B. If we average over all sites belonging to the same type, as shown in Fig. 10, we see that jumps between the octahedral sites are fastest. From the site splitting of the 48g into 96h sites follows that two sites are present inside each octahedron, and that there is a free energy barrier between the split sites. Our results are therefore in agreement with ab initio calculations Xu et al. (2012a); Jalem et al. (2013); Miara et al. (2013); Meier et al. (2014), that show that the minima in the Li-ion potential energy surface are displaced from the original central site in the octahedral site. The presence of two distinct but very close sites manifests in very high exchange rates between these two. LLZO also displays fast jumps from the tetrahedral into the octahedral environment, whereas the reverse jump takes three to four times longer. No jumps between tetrahedral environments are observed, as expected, since an ion needs to traverse octahedral sites to reach a different tetrahedral site. While the jump probabilities, or lag times, are non-symmetric, the fluxes are symmetric, which is necessary to observe local detailed balance.

The diffusive pathways estimated from the algorithms are shown in Fig. 11. The connectivity analysis reveals the existence of a single dominant pathway that allows mobile ions to diffuse through the entire simulation cell. The edge widths in the figure are proportional to the observed flux of particles, and we see that, where the octahedral site splitting is correctly determined, there is a large flux of ions between split octahedral sites compared to the smaller flux between the tetrahedral and octahedral environments.

Summarizing our results for this material, our site analysis finds the splitting of the 48g to 96h sites, which sets it apart from any density-based analysis. An analysis based on distance criteria to crystallographic sites would have worked as well or better, but obviously requires prior knowledge.

III.2 Analysis of LiAlSiO4

The structure of the $\beta$-eucryptite \ceLiAlSiO4 William W. Pillars and Donald R. Peacor (1973), referred to as LASO hereafter, is taken from COD Gražulis et al. (2012) entry 9000368. It has been studied for its anisotropic expansion coefficient Lichtenstein et al. (1998, 2000) and its ionic conductivity Schink and Löhneysen (1983); Donduft et al. (1988); Sartbaeva et al. (2005); Chen et al. (2018). The structure can be described as an ordered $\beta$-quartz solid solution, with alternating aluminum and silicon planes. Location and occupation of the sites for lithium have been contested. In the original reference William W. Pillars and Donald R. Peacor (1973), the difficulties in determining the lithium sites in previous and in the same work are explained very well. For example, earlier work Winkler (1948) concluded that the Li sites are coplanar with the Al sites, while Pillars and Peacor William W. Pillars and Donald R. Peacor (1973) show that the lithium sites are also present in the Si plane. Later work Donduft et al. (1988) shows that both sites are available to lithium and establishes the unidimensional chain of these sites as the mechanism for ionic diffusion in this material. There is now a better understanding of this structure and the sites available to lithium, but the original CIF-file in the COD originating from the experiments by Pillars and Peacor William W. Pillars and Donald R. Peacor (1973) does not list all sites. Any analysis that relies on this knowledge would therefore have failed. Our molecular dynamics simulations and subsequent site analysis yield results that are compatible with the latest literature regarding the ionic transport in this material.

We simulate with first principles \ceLi12Al12Si12O48, starting from the reported CIF-file Gražulis et al. (2012). A full atomic and cell optimization results in a volume increase of 3.6% without changing the cell angles in a significant way. We perform the subsequent dynamical simulations using Born-Oppenheimer molecular dynamics in the canonical ensemble at a temperature of 750 K for 291 ps, with further details given in Appendix D.

We show in Fig. 12 the Li-ion densities sampled during the dynamics. The unidimensional channels of ionic diffusion are compatible with published results Donduft et al. (1988); Chen et al. (2018). The diffusion coefficient is hard to converge for the short dynamics we obtained for this system, and so quantifying the diffusion coefficient and its error cannot be done rigorously. We plot the coordinates of Li-ions as a function of time in Fig. (2) in Ref. kah to show that motion along the z-coordinate is observed during the simulation, compatible with long-range diffusion. The RDFs of Li with all present species, shown in the upper panel of Fig. 13, display a first coordination shell composed by four oxygen ions, compatible with literature findings. A second and third shell are composed of silicon and aluminum, with the amplitude of Si being stronger in the second shell, and Al in the third shell. This hints that Li ions prefer sites in the Si plane to those in the Al plane.

When running the site analysis, we find 24 sites of two different types, twelve of type 0 and twelve of type 1, compatible with the latest literature results Donduft et al. (1988). The parameters are as given in Appendix E, except for a cutoff midpoint of 1.3 instead of 1.5, which is more robust with respect to the total number of sites obtained. The clustering analysis in Fig. (5) in Ref. kah shows that types can be distinguished easily. An analysis of the RDF for the individual site types, shown in the bottom panel of Fig. 13, reveals that the discriminant is the different coordination of Al and Si, which is detected by the SOAP descriptor. For type 0, the second shell is composed of two aluminum atoms; four silicon atoms are in the third shell. For type 1, the numbers are the same, but silicon is replaced by aluminum and vice versa. The RDF in the upper image of Fig. 13 hints at the fact that the Li ions prefer to occupy type 1 sites where the Si ions are closer than the Al ions. From the site analysis of our simulation, we calculate the mean occupation ratio to be 77% for site type 1 and 23% for site type 0. Literature reports give occupancies of 68% and 22% Press et al. (1980), respectively, or a 3:1 ratio Lichtenstein et al. (2000). In Fig. 14, we show the jump lag between the type sites. Jumps from type 1 to type 0 are about 3.5 times faster, which is necessary to preserve detailed balance. We observe no jumps within the site types, which is expected since the sites’ types are alternating along the diffusion channels (see Fig. 15).

We can thus show with first-principles molecular dynamics and an unsupervised analysis that the lithium ions occupy two different site types in LASO. This is done without any knowledge of the possible sites, since in the original CIF file only twelve sites (for twelve lithium ions) are given. This example highlights the challenges for algorithms that rely on prior knowledge of crystallographic sites. Such information might be missing or wrong, for example, because of the difficulties of resolving low occupancy sites for light elements when using XRD or neutron diffraction, or because of simulation conditions (e.g., temperature) differing from the experimental setup. Relying on all the sites being known can obviously be problematic in some cases. An unsupervised approach requiring minimal knowledge of the structures should be preferred in such cases. Unlike the case of LLZO, a density-based clustering on the lithium-ion positions would very likely also have given the same results, as can be conjectured from the lithium-ion densities in Fig. 12, where the highest isovalue clearly shows disconnected regions of high ionic density.

III.3 Analysis of tetragonal Li10GeP2S12

The superionic conductor \ceLi10GeP2S12 (LGPS) in its tetragonal phase was first reported by Kamaya et al. Kamaya et al. (2011) Its unprecedented ionic conductivity at room temperature motivated studies of its diffusion mechanisms using atomistic simulation techniques Adams and Rao (2012a); Mo et al. (2012); Xu et al. (2012b); Marcolongo and Marzari (2017). The original paper Kamaya et al. (2011) reports three site types in the unit cell: tetrahedrally coordinated 16h, tetrahedrally coordinated 8f, and octahedrally coordinated 4d sites, all coordinated with sulfur, with only the latter possessing full occupancy. The 16h and 8f sites denote edge-sharing \ceLiS4 tetrahedra that form one-dimensional channels along the c-axis, the main diffusive pathways Kamaya et al. (2011); Mo et al. (2012). Adams and Rao Adams and Rao (2012a) found evidence for an additional four-fold coordinated site – termed 4c – using classical simulations, which was validated in subsequent experiments by Kuhn et al. Kuhn et al. (2013a)

We analyze the first-principles molecular dynamics trajectories for LGPS that were produced for a recent work discussing the failure of the Nernst-Einstein relation in this structure Marcolongo and Marzari (2017). We refer to the reference for computational details, and state only that the trajectories were run with the cp.x code of the Quantum ESPRESSO distribution Giannozzi et al. (2009) with a PBE-exchange correlation functional Perdew et al. (1996). Using a unit cell of 50 atoms, 428 picoseconds of dynamics were obtained in the microcanonical ensemble, after an equilibration run at a target temperature of 500 K. We find a Li-ion density, shown in Fig. 16, that is compatible with literature results on the unidimensional channels Adams and Rao (2012a) that dominate the diffusion in this material. The diffusion in this material, calculated from the mean-square displacement, shown in Fig. (3) in Ref. kah , is $D_{tr}^{Li}=\mathrm{3.25\times 10^{-6}\,cm^{2}\,s^{-1}}$, compatible with literature results. For example, Kuhn et al. Kuhn et al. (2013b) report a value of the Li-ion tracer diffusion coefficient of $D_{tr}^{Li}\approx 10^{-6}\,\mathrm{cm^{2}\,s^{-1}}$ at 500 K, which is close to our estimate and certainly within the likely error bounds of FPMD that stem from, among other factors, short simulations in small unit cells.

The landmark analysis is applied to the equilibrated trajectory to determine statistics. We treat germanium and phosphorus atoms as one species since the 4d tetrahedral site is occupied by either species to avoid identifying extraneous site types due to the arbitrary choice of occupation of these sites. We will refer to both phosphorus and germanium as phosphorus hereafter. After site detection and SOAP clustering, shown in Fig. (4) in Ref. kah , we find 30 sites of type 0, 24 sites of type 1, and four sites of type 2. We see that type 2 corresponds to the octahedral environment of the 4d site. To understand the difference between the different site types, we calculate the RDF for every site type, shown in the middle and bottom panels of Fig. 17 for sulfur and phosphorus. A visual depiction of where the sites are located is shown in Fig. 18.

In the RDFs between lithium and phosphorus, key differences appear between the different site types. While site type 1 is compatible with four-fold coordination with sulfur, site types 0 and 2 tend to plateau towards a coordination with six sulfur atoms, which is expected only for the latter site type. There is no evidence for a six-fold coordinated site type inside the ion-conducting channel of LGPS. We should note that — to our knowledge — no analysis has yet been done on dynamically short-lived features of the coordination of lithium with sulfur in LGPS, so it is possible that the features we perceive in our analysis are not detected when studying averages. However, we also observe that the algorithm is less robust than for the other studied examples. The number of sites as well as the clustering to types depend in this case more strongly on the parameters chosen for the site analysis. The very similar atomic environment of different site types leads to large overlap of clusters of the SOAP vectors, shown in Fig. (4) in Ref. kah . The superionic behavior of Li ions in LGPS impedes the precise definition of a site for any single mobile ion in the dynamic potential energy landscape. LGPS, representative of superionic systems, can be seen as a worst-case scenario for the present site analysis.

When analyzing LGPS it also becomes evident that the classification via SOAP vectors can yield different results than the Wyckoff symbols resulting from symmetry analysis. Different Wyckoff positions can be classified as the same site type if their chemical and geometric environments are too similar to differentiate. Further, as a result of symmetry breaking during molecular dynamics, two sites with the same Wyckoff position can be classified as different types, especially in non-ergodic simulations. This is not necessarily a weakness of the analysis, but something to be aware of. We note that despite the obvious difficulties in detecting sites and site types reported in the literature, our analysis found four off-channel four-fold coordinated sites, which were termed 4c sites by Adams and Rao Adams and Rao (2012a). These sites had been missed in earlier FPMD simulations Mo et al. (2012), since they were not reported by preceding experiments. Thus an unsupervised and unbiased analysis can help when experimental data is lacking or incomplete.

III.4 Non-diffusive structures

To further validate the method, we additionally study seven non-conductive structures. The structures were selected from an ongoing screening effort intended to find new solid-state electrolytes, and chosen from the least conductive systems that had two site types in their CIF file. For every structure, a molecular dynamics simulation is run at a temperature of 1000 K, with simulation lengths long enough to estimate the diffusivity of the material. All other simulation parameters are the same as those presented for LASO in Sec. III.2; details can be found in Appendix D. The landmark analysis is run on every second frame of the trajectory (about every 60 fs). We use the same default landmark analysis parameters for all of the materials, with further details given in Appendix E.

The results can be seen for each of the seven materials in Table 1. For all but two materials, the landmark analysis produces the same number of sites and the same division of those sites into types as given in the corresponding CIF files. For these materials, unlike LGPS, the Wyckoff analysis and the SOAP analysis coincide. In \ceLi20Re4N16 and \ceLi24Zn4O16, however – like in \ceLiAlSiO4 – unsupervised landmark analysis identifies sites that are not present in the CIF files from ICSD (see Fig. 19). In \ceLi20Re4N16, these four sites complete the planar connected components in the material; they are transitional sites with low occupancy and residence time. Their existence is confirmed by an analysis of the Li-ion densities observed in the trajectories. In \ceLi24Zn4O16, Li-ions from neighboring sites occasionally and briefly jump to the additional sites and then back. The additional sites again have low occupancy and residence time and are confirmed by a density analysis of the real-space coordinates.

IV Implementation

The landmark analysis presented here is implemented as a component of SITATOR Musaelian and Kahle (2018), a modular, extensible, open-source Python framework for analyzing networks of sites in molecular dynamics simulations of solid-state materials. SITATOR provides two fundamental data structures: SITENETWORK, which represents possible sites for some mobile atoms in a host lattice and SITETRAJECTORY, which stores discretized trajectories for those mobile atoms. A SITENETWORK can also store arbitrary site and edge attributes. SITATOR includes an optimized implementation of landmark analysis as well as pre-processing utilities for trajectories and tools for analyzing and visualizing the results of site analyses.

V Conclusions

We presented a novel method to perform a site analysis of molecular dynamics trajectories to analyze ionic diffusion in solid-state structures. The method is robust and can run over a large range of materials with a minimal set of parameters and little human intervention. As we have shown, our landmark analysis performs well where other methods fail, whether because of very high exchange rates and/or close proximity between sites (as in LLZO), or because needed prior information is missing (as in LASO, \ceLi20Re4N16, and \ceLi24Zn4O16 where several sites that were occupied during our simulations are not given in the experimental CIF file). As became evident for LGPS, superionic conductors with a liquid-like, highly disordered lithium sublattice are hard to analyze with the tool, and the signals from the analysis need to be studied in further detail in subsequent work. A suggestion for subsequent work is automatically computing the configurational entropy descriptor $\tilde{S}$ described by Kweon et al. Kweon et al. (2017) While the method presented here will not necessarily outperform carefully chosen analysis tools with parameters specific to the system under investigation, it has advantages when comparing different systems and in high-throughput applications, such as the search for microscopic descriptors for ionic diffusion in the solid state. Another suggestion is to study the collective motion in common superionic conductors from occupation statistics given by the landmark analysis. Concerted motion is important for ionic diffusion in a wide class of systems He et al. (2017), and an analysis that can be used with the same set of parameters on a wide range of materials can be used to quantify collective effects rigorously.

Acknowledgments

We thank Aris Marcolongo for supplying the LGPS trajectories. We also express our gratitude to Matthieu Mottet who provided force-field parameters for LLZO and advice on running the simulations. We would like to thank Felix Musil, Piero Gasparotto, and Michele Ceriotti for their support with the SOAP descriptor and the QUIPPY interface as well as fruitful discussions. We gratefully acknowledge financial support from the Swiss National Science Foundation (SNSF) Project No. 200021-159198. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID mr0.

Appendix A Landmark vector clustering algorithm

When determining the landmarks in a system, we use the efficient implementation of the Voronoi decomposition from ZEO++ Willems et al. (2012); Martin et al. (2012), which accounts for periodic boundary conditions. A custom hierarchical agglomerative algorithm is used to cluster the landmark vectors. The algorithm is designed for streaming: no pairwise distance matrix is ever computed or stored, and the landmark vectors can be streamed from disk in the order they were written, avoiding random access. Clusters are represented by their average landmark vectors, called centers. At each iteration clusters whose centers are sufficiently similar are merged. After a small number of iterations, a steady state is reached when no clusters can be merged; this is taken as the final clustering. The original landmark vectors are then each assigned to the cluster whose center they are most similar to.

Two parameters control the characteristics of the landmark clustering: the clustering threshold, which determines how aggressively new clusters (sites) should be added, and the minimum cluster size, which filters out sites whose occupancy is extremely low (such clusters likely represent thermal noise or transitional states). These parameters allow the user to tune spatial and temporal resolution. Specifically:

Set the initial cluster centers $\mathbf{c}_{i}$ to the landmark vectors. The order of the landmark vectors does affect the clustering, but in practice we have found the effect to be minimal. We process the landmark vectors in the order they were generated: chronologically and in whatever order the mobile ions were numbered. 2. 2.

Take the first existing cluster center $\mathbf{c}_{0}$ as the first new cluster center $\mathbf{c}^{\prime}_{0}$. Then, for each remaining cluster center $\mathbf{c}_{i}$, $i\in[1,N)$:

1. (a)

Find the new cluster center $\mathbf{c}^{\prime}_{j}$ to which the old cluster center $\mathbf{c}_{i}$ is most similar:

[TABLE]

where $N^{\prime}$ is the current number of new cluster centers and $S(\cdot,\cdot)$ is the normalized cosine metric:

[TABLE] 2. (b)

If

[TABLE]

then merge the old cluster $\mathbf{c}_{i}$ into the new cluster $\mathbf{c}^{\prime}_{j}$:

[TABLE]

where $n$ is the total number of old clusters that have been merged to form $\mathbf{c}^{\prime}_{j}$ so far.

Otherwise, keep $\mathbf{c}_{i}$ as the center of its own cluster:

[TABLE] 3. 3.

Repeat the previous step until no further clusters can be merged; the $\mathbf{c}_{i}$, $i\in[0,N)$ are the final clusters. 4. 4.

Assign the landmark vectors to clusters. The assignment threshold controls how dissimilar a landmark vector can be to its cluster’s center before it is marked as unassigned. This parameter controls the trade-off between spatial accuracy and the proportion of unassigned mobile atom positions: high values will give greater spatial precision, while lower values will ensure that almost all mobile atoms are assigned to sites at all times.

For each landmark vector $\mathbf{l}$:

1. (a)

Find the most similar cluster center:

[TABLE] 2. (b)

If $s>\text{assignment threshold}$, then mark $\mathbf{l}$ as a member of the corresponding cluster with confidence $s$.

Otherwise, mark $\mathbf{l}$ as unassigned. 5. 5.

Remove clusters smaller than the minimum cluster size. 6. 6.

Repeat step 4 with the remaining clusters, yielding the final cluster assignments.

Appendix B Markov clustering

We apply Markov Clustering Van Dongen (2008) to the matrix $\bm{M}$ to simulate biased random walks through a graph, giving preference to high-probability routes. Once the process converges, a set of internally highly connected subgraphs remains. The sites in each resulting subgraph, if there are more than one, are merged into a single site. Their real-space positions are averaged, and the mobile ions that occupied any of the merged sites now occupy the new site.

We use typical Markov Clustering parameters of $2.0$ for both expansion and inflation. We do not add artificial self loops to the graph since $\bm{M}$ already contains appropriate nonzero values on the diagonal.

Appendix C Parameter estimation for Density-Peak clustering

Density-peak clustering Rodriguez and Laio (2014) defines the number of clusters as the number of data points with extreme outlier values of $\rho$ (density) and $\delta$ (distance to nearest neighbor with larger $\rho$), as determined by a user-specified threshold. Rodriguez and Laio Rodriguez and Laio (2014) suggest a simple heuristic for determining this threshold that we adopt and automate. First, we compute the values $\gamma_{i}=\rho_{i}\delta_{i}$ and sort them into decreasing order. In a well behaved clustering problem, a plot of $\gamma$ then has a recognizable “elbow,” and the points before the elbow – before the curve rapidly flattens out – are the outliers. Thus the problem of determining the thresholds is equivalent to finding the elbow of this curve.

We use a simplified version of the knee-finding algorithm presented by Satopaa et al. Satopaa et al. (2011) A straight line is taken between $(0,\gamma_{0})$ and $(n,\gamma_{n})$, and the point $(i,\gamma_{i})$ with the maximum distance to that line is taken as the elbow. The $\rho$ and $\delta$ values corresponding to that point are then used as the thresholds for the density-peak clustering.

Appendix D Molecular dynamics parameters

The simulations for LASO and the seven non-diffusive structures are performed with the pw.x module in the Quantum ESPRESSO distribution Giannozzi et al. (2009), using pseudopotentials and cutoffs from the SSSP Efficiency library 1.0 Prandini et al. (2018). The exchange-correlation used in the DFT is PBE Perdew et al. (1996). The materials informatics platform AiiDA Pizzi et al. (2016) is used to ensure full reproducibility of the results and achieve a high degree of automation.

We always perform a variable-cell relaxation prior to the molecular dynamics, with a uniform k-point grid of 0.2 Å*-1* and no electronic smearing since we consider only electronic insulators. The energy and force convergence thresholds are $0.5\times 10^{-4}$ and $0.25\times 10^{-5}$ in atomic units, respectively. We set the pressure threshold to 0.5 kbar. A meta-convergence threshold on the volume, which specifies the relative volume change between subsequent relaxations, is set to 0.01.

We create supercells with the criterion that the minimal distance between opposite faces is always larger than 6.5 Å. We run the molecular dynamics simulations with a stochastic velocity rescaling thermostat Bussi et al. (2007) which we implemented into Quantum ESPRESSO, with a characteristic time of the thermostat set to 0.2 ps at constant volume and number of particles (NVT ensemble). The timestep is set to 1.45 fs, and snapshots of the trajectory are taken every 20 time steps.

The origin of the non-diffusive structures is given in Table 2, together with the simulation time. The structures are taken from the Inorganic Crystallography Open Database (ICSD) Belsky et al. (2002) and the Open Crystallography database (COD) Gražulis et al. (2012).

Appendix E Site analysis parameters

Unless otherwise indicated, the landmark analysis uses a cutoff midpoint of $d_{0}=1.5$ and steepness of $k=30$, a minimum site occupancy of 1%, and landmark clustering and assignment thresholds of $0.9$. For computing SOAP descriptors, unless otherwise specified, we use a Gaussian width of 0.5 Å on the atomic positions, a cutoff transition width of 0.5 Å, and spherical harmonics up to $n_{max}=l_{max}=6$. The radial cutoff is set to always include the nearest neighbor shell of all other species (excluding the mobile species). We calculate SOAP vectors for mobile ions every tenth frame, and average every 10 SOAP vectors to reduce noise. The principal components of the averaged SOAP vectors are extracted using Principal Component analysis (PCA) to retain at least 95% of the variance, and the clustering is performed in this reduced space.