# Atomistic modeling of bulk and grain boundary diffusion in solid electrolyte
Li6PS5Cl using machine-learning interatomic potentials

Yongliang Ou<EMAIL_ADDRESS>Institute for Materials
Science, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany
Yuji Ikeda<EMAIL_ADDRESS>Institute for Materials Science,
University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany  Lena
Scholz  Institute of Applied Mechanics, University of Stuttgart,
Universitätsstr. 32, 70569 Stuttgart, Germany  Sergiy Divinski  Institute of
Materials Physics, University of Münster, Wilhelm-Klemm-Str. 10, 48149
Münster, Germany  Felix Fritzen  Institute of Applied Mechanics, University of
Stuttgart, Universitätsstr. 32, 70569 Stuttgart, Germany  Blazej Grabowski
Institute for Materials Science, University of Stuttgart, Pfaffenwaldring 55,
70569 Stuttgart, Germany

###### Abstract

Li6PS5Cl is a promising candidate for the solid electrolyte in all-solid-state
Li-ion batteries due to its high ionic conductivity. In applications, this
material is in a polycrystalline state with grain boundaries (GBs) that can
affect ionic conductivity. While atomistic modeling provides valuable
information on the impact of GBs on Li diffusion, such studies face either
high computational cost (when using ab initio methods) or accuracy limitations
(when using classical potentials) as challenges. Here, we develop a quality-
level-based active learning scheme for efficient and systematic development of
ab initio-based machine-learning interatomic potentials, specifically moment
tensor potentials (MTPs), for large-scale, long-time, and high-accuracy
simulations of complex atomic structures and diffusion mechanisms as
encountered in solid electrolytes. Based on this scheme, we obtain MTPs for
Li6PS5Cl and investigate two tilt GBs, $\Sigma 3(1\bar{1}2)[110]$, $\Sigma
3(\bar{1}11)[110]$, and one twist GB, $\Sigma 5(001)[001]$. All three GBs
exhibit low formation energies of less than
$20\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{2}$,
indicating their high stability in polycrystalline Li6PS5Cl. Using the MTPs,
diffusion coefficients of the anion-ordered and anion-disordered bulk, as well
as the three GBs, are obtained from molecular dynamics simulations of
atomistic models with more than
$16\,000\text{\,}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}\mathrm{s}$ for
$5\text{\,}\mathrm{ns}$. At $300\text{\,}\mathrm{K}$, the GB diffusion
coefficients fall between the ones of the anion-ordered bulk structure
($1.2\text{\times}{10}^{-9}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$,
corresponding ionic conductivity about
$0.2\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}$) and the
anion-disordered bulk structure ($50\text{\,}\mathrm{\char 37\relax}$
Cl/S-anion disorder;
$2.2\text{\times}{10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$,
about $29.8\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}$) of
Li6PS5Cl. Experimental data fall between the Arrhenius-extrapolated diffusion
coefficients of the investigated atomic structures, supporting our
quantitative in silico predictions.

##  I Introduction

All-solid-state Li-ion batteries have attracted attention for their improved
safety compared to conventional Li-ion batteries by virtue of their solid
instead of a flammable liquid electrolyte [1]. However, finding a suitable
material for the solid electrolyte is not trivial [2]. In 2008, the
argyrodite-type Li6PS5Cl was first reported, featuring an unusually high Li-
ion mobility [3]. Later, intensive experimental investigations [4, 5, 6]
affirmed the superior ionic conductivity, and ever since Li6PS5Cl is
considered a candidate for the solid electrolyte.

Diffusion of Li ions in Li6PS5Cl was measured using nuclear magnetic resonance
[7, 8]. The room-temperature diffusion coefficients reported in the two
studies are of the same order of magnitude
($3.87\text{\times}{10}^{-8}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
[7] vs.
$2.5\text{\times}{10}^{-8}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
[8]), and the derived activation energies differ by about 20%, $0.35\pm
0.01\text{\,}\mathrm{e}\mathrm{V}$ [7] vs. $0.28\pm
0.01\text{\,}\mathrm{e}\mathrm{V}$ [8]. Differences in the activation energies
can come from the relatively narrow temperature intervals of the diffusion
measurements and difficulties in reproducibly synthesizing such materials.
Depending on the synthesis conditions, polycrystalline materials with
Cl/S-anion disorder [9] and potentially different grain boundary (GB)
distributions [10, 11] are obtained. Further, the contributions of crystalline
bulk and GBs to the Li diffusion can hardly be separated in experiments [10],
which might affect the measured coefficients.

Atomistic simulations provide an important complementary tool for an improved
understanding of diffusion mechanisms and the impact of GBs. Table 1
summarizes molecular dynamics (MD) simulations for Li self-diffusion in
Li6PS5Cl. The diffusion coefficient for the anion-ordered bulk structure was
computed in several previous studies [12, 13, 14, 15], all of which gave
values that are orders of magnitude smaller than those in experiments. Further
investigations (both simulations [16, 17] and experiments [9]) showed that the
origin of superionic Li diffusion is most likely related to the Cl/S-anion
disorder in the bulk structure of Li6PS5Cl. This anion disorder triggers Li
inter-cage diffusion. Recently, quantitative investigations of Li diffusivity
in the anion-disordered structure were also performed by ab initio MD
simulations [18, 19, 15]. With $50\text{\,}\mathrm{\char 37\relax}$ Cl/S-anion
disorder in the bulk structure, the Arrhenius-extrapolated Li diffusion
coefficients at room temperature are consistent with experiments [18, 15]. The
simulated activation energies of Li diffusion are lower
($0.200.26\text{\,}\mathrm{e}\mathrm{V}$ [18] and
$0.25\text{\,}\mathrm{e}\mathrm{V}$ [15]) than the above-listed experimental
data. Due to the typical high computational requirements of ab initio MD
simulations [20], the sampling time was relatively short in these previous
studies, and defects (e.g., GBs) were often neglected.

Table 1: Theoretical studies which provided calculated diffusion coefficients for Li self-diffusion in solid electrolyte Li6PS5Cl. A closely related Li6PS5Br study [19] including grain boundaries (GBs) is also listed. Previous studies used ab initio molecular dynamics (MD) simulations, while in the present work, a machine-learning interatomic potential (MLIP), specifically moment tensor potential (MTP), fitted to ab initio data is used to accelerate MD. The number of atoms in the simulation cell is denoted by $N_{\textrm{at}}$. “Yes” or “no” means that the corresponding structures have or have not been considered in the given study.  Year | Sampling | Size ($N_{\mathrm{at}}$) | Time (ns) | Anion-ordered bulk | Anion-disordered bulk | GBs | Reference
---|---|---|---|---|---|---|---
2017 |  ab initio MD | $52$ | $0.1$ | yes | no | no | [12]
2019 |  ab initio MD | $52$ | $0.3$ | yes | partially111Only one configuration with a Cl/S antisite defect was considered in Ref. [13]. | no | [13]
2022 |  ab initio MD | $52$ | $0.12$ | yes | no | no | [14]
2022 |  ab initio MD | $52$ | $0.15$ | yes | yes | no | [18]
2022 | classical force field | $416$ | $20$ | no | yes | no | [21]
2023222A structurally similar material, Li6PS5Br, was investigated in Ref. [19]. |  ab initio MD | $312$ | $0.04$ | yes | yes | yes | [19]
2024 |  ab initio MD | $52$ | $0.3$ | yes | yes | no | [15]
2024 |  ab initio $\rightarrow$ MLIP MD | $>16\,000\text{\,}$ | $5$ | yes | yes | yes | this work

Interatomic potentials can be utilized instead of ab initio simulations to
reduce the computational cost. This enables large-scale and long-time MD
simulations, equipping us with tools for a statistically reliable prediction
of diffusion properties. For example, classical force fields parameterized by
ab initio data were utilized for simulating solid electrolytes [22, 23, 21].
Due to the complex structure of these materials and the limited number of
fitting parameters in the force fields, simulation results with lower accuracy
have to be expected. A good alternative is given by machine-learning
interatomic potentials calibrated to ab initio data, which have recently
emerged as a powerful tool to accelerate MD simulations while preserving near
ab initio accuracy [24, 25]. Machine-learning interatomic potentials were
shown to accurately describe diffusion [26, 27] and work well even for complex
electrolytes [28]. Li diffusion in Li6PS5Cl has not been investigated
systematically with machine-learning interatomic potentials to our knowledge.

The present study investigates the impact of three GBs with different
structural characteristics ($\Sigma 3(1\bar{1}2)[110]$, $\Sigma
3(\bar{1}11)[110]$, and $\Sigma 5(001)[001]$) on Li diffusion in the solid
electrolyte Li6PS5Cl, aiming to provide accurate GB diffusion coefficients and
thereby enhance the understanding of experimentally measured data. To this
end, we propose and apply an active learning scheme that systematically
exploits the quality levels of machine-learning interatomic potentials,
specifically of moment tensor potentials (MTPs) [29]. Based on the scheme, we
fit machine-learning potentials to ab initio data and systematically analyze
their performance using the target quantity, i.e., the diffusion coefficient,
as a measure. Accelerated by the thus obtained and validated machine-learning
interatomic potentials, large-scale and long-time MD simulations are performed
to optimize the GB structures and analyze the difference in GB diffusion
mechanisms compared to the bulk. From these simulations, we extract the
diffusion coefficients for the GBs as well as for the anion-ordered and anion-
disordered bulk structures. The results provide a theoretically admissible
range of Li diffusivity in Li6PS5Cl.

##  II Methodology

###  II.1 Grain boundary construction

Li6PS5Cl exhibits an argyrodite-type structure and belongs to the cubic
crystal system with the space group $F\bar{4}3m$ (No. 216) [3, 30]. Figure
1(a) shows the conventional unit cell of the bulk structure, which contains
four formula units (52 atoms). The shown bulk structure is anion-ordered,
i.e., the Wyckoff sites $4a$ [the sites symmetrically equivalent to (0, 0, 0)]
and $4c$ [equivalent to (0.25, 0.25, 0.25)] are fully occupied by Cl and S
atoms, respectively. Further, PS4 tetrahedral units are present with the P and
S atoms located at $4b$ [equivalent to (0.5, 0.5, 0.5)] and $16e$ sites,
respectively. The Cl and the P atoms each form a face-centered cubic
sublattice. These relatively heavy atoms (P, S, and Cl) comprise a three-
dimensional backbone that provides multiple interstitial sites to accommodate
the lightweight Li atoms. The Li atoms are positioned at $24g$ sites in the
ideal model with high symmetry. As highlighted in Fig. 1(a), six Li atoms
octahedrally coordinate an S atom at a $4c$ Wyckoff site, forming a so-called
Li cage.

It was reported both experimentally [30, 31] and theoretically [32, 33, 34]
that the high-symmetry $F\bar{4}3m$ structure is not stable at low
temperatures and that neighboring interstitial sites, e.g., $48h$ Wyckoff
sites, offer lower-energy states for the Li atoms. This results in symmetry-
breaking displacements of the Li atoms, referred to as tilting of the Li
cages. The dynamic stabilization of the Li atoms at the high-symmetry
positions (phase transition) with increasing temperature was investigated in
Ref. [31]. Such a temperature-induced phase transition is also observed in
many other materials. For example, the high-symmetry structure of BaFeO3
perovskite is dynamically unstable at lower temperatures and becomes
dynamically stabilized at elevated temperatures [35].

Figure 1: (a) Conventional unit cell of anion-ordered bulk Li6PS5Cl with 52
atoms (four formula units). One octahedral cage consisting of six Li atoms is
highlighted in blue. (b) Coincidence-site lattice (CSL) unit cell with 78
atoms (six formula units). The planes for constructing the two $\Sigma 3$ GBs
are highlighted in blue. (c) CSL unit cell with 130 atoms (ten formula units).
The plane for constructing the $\Sigma 5$ GB is highlighted in blue. The CSL
unit cells are used to construct the GB simulation models (Table 2).
Visualization performed using vesta [36].

Based on the cubic symmetry of Li6PS5Cl and the coincidence-site lattice
theory, periodic and commensurate GB structures can be generated [37].
Considering the complexity of the bulk structure, three relatively simple GBs,
yet featuring different geometrical arrangements, namely $\Sigma
3(1\bar{1}2)[110]$, $\Sigma 3(\bar{1}11)[110]$, and $\Sigma 5(001)[001]$, were
chosen for the present study. Table 2 lists information about the simulation
models for the GBs and also the bulk. $\Sigma 3(1\bar{1}2)[110]$ and $\Sigma
3(\bar{1}11)[110]$ are tilt GBs symmetric for the Cl and P atoms with
misorientation angles of $70.53\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$109.47\text{\,}\mathrm{\SIUnitSymbolDegree}$, respectively. The tilt axis is
$[110]$, and the GB planes are $(1\bar{1}2)$ and $(\bar{1}11)$. $\Sigma
5(001)[001]$ is a twist GB with a rotational axis of $[001]$ and a
misorientation angle of $36.87\text{\,}\mathrm{\SIUnitSymbolDegree}$. The GB
plane is $(001)$ and lies perpendicular to the rotational axis.

Table 2: Information on the anion-ordered bulk and the GB models. The $x$, $y$, and $z$ entrees give the relation of the Cartesian axes of each model with the crystallographic directions of the conventional unit cell of the bulk structure [Fig. 1(a)]. For all GB simulation cells, $x$ and $y$ are parallel ($\parallel$) and $z$ is perpendicular ($\perp$) to the GB plane. The tilt axis of the $\Sigma 3$ tilt GBs is along the $y$ direction, and the rotational axis of the $\Sigma 5$ twist GB is along the $z$ direction. The total number of atoms is denoted by $N_{\textrm{at}}$. The dimensions of the optimized structures, the GB formation energies [$\gamma$, Eq. (8)] calculated at $0\text{\,}\mathrm{K}$ (see Sec. III.2), and the GB widths ($\delta$) evaluated at $600\text{\,}\mathrm{K}$ (see Secs. III.2 and IV) are shown. The excess volumes per unit GB area of all three GBs are less than $1\text{\,}\mathrm{\text{Å}}$.  Model |  $x$ $(\parallel)$ |  $y$ $(\parallel)$ |  $z$ $(\perp)$ | $N_{\mathrm{at}}$ | Dimension ($\text{\,}\mathrm{\text{Å}}$) |  $\gamma$ (meV/Å2) |  $\delta$ (Å)
---|---|---|---|---|---|---|---
anion-ordered bulk | $[100]$ | $[010]$ | $[001]$ | $17\,836\text{\,}$ | $70\times 70\times 70$ | n/a | n/a
GB $\Sigma 3(1\bar{1}2)[110]$ tilt $70.53\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $[\bar{1}11]$ | $[110]$ | $[1\bar{1}2]$ | $19\,656\text{\,}$ | $53\times 51\times 148$ | $10.44$ | $28.8$
GB $\Sigma 3(\bar{1}11)[110]$ tilt $109.47\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $[1\bar{1}2]$ | $[110]$ | $[\bar{1}11]$ | $17\,472\text{\,}$ | $50\times 51\times 133$ | $7.12$ | $16.7$
GB $\Sigma 5(001)[001]$ twist $36.87\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $[130]$ | $[3\bar{1}0]$ | $[001]$ | $16\,380\text{\,}$ | $49\times 49\times 139$ | $18.78$ | $26.5$

To obtain the three GB supercells, two coincidence-site lattice unit cells
[Fig. 1(b) for the two $\Sigma 3$ GBs, and Fig. 1(c) for the $\Sigma 5$ GB]
were constructed by coordinate transformations of the conventional unit cell
[Fig. 1(a)]. The GB planes used for constructing the GBs are highlighted in
blue. For $\Sigma 3(1\bar{1}2)[110]$, the GB plane cuts through some PS4 units
and some Li cages. For $\Sigma 3(\bar{1}11)[110]$, in contrast, the integrity
of all Li cages and PS4 units is preserved by setting the GB plane to the
position of the Cl atomic layers. The integrity of all Li cages is also
preserved for $\Sigma 5(001)[001]$, but some PS4 units are cut through by the
GB plane. With these distinct structural features of the GBs, different
diffusion behaviors of Li is expected.

The GB simulation models were constructed as periodic bicrystals, where two
grains of different crystallographic orientations are stacked in the $z$
direction. The two grains were slightly separated to ensure that the smallest
atomic distance in the GB simulation cells is larger than
$1.5\text{\,}\mathrm{\text{Å}}$. In each GB simulation cell, $x$ and $y$ lie
within the GB plane ($\parallel$), and $z$ is perpendicular ($\perp$) to the
GB plane. The tilt axis of the $\Sigma 3$ GBs is in the $y$ direction, and the
rotational axis of the $\Sigma 5$ GB is in the $z$ direction. There are two
GBs of the same type per simulation cell but with inverted orientations due to
the coincidence-site lattice construction. To avoid the interaction between
these GBs, supercells were made so that the distances between the GBs in the
direction perpendicular to the GB planes ($z$) are larger than
$65\text{\,}\mathrm{\text{Å}}$. Table 2 lists the total number of atoms for
each GB simulation cell. All simulation cells maintain the stoichiometry of
Li6PS5Cl, ensuring formal charge neutrality.

###  II.2 Anion disorder in bulk structure

Anion disorder in Li6PS5Cl refers to the mixing of Cl and S atoms at the $4a$
and the $4c$ Wyckoff sites of the bulk structure, as observed in experiments
[3, 30]. In the present study, $50\text{\,}\mathrm{\char 37\relax}$ mixing of
Cl and S was considered. The atomistic model with the anion-disordered bulk
structure was constructed based on that with the anion-ordered bulk structure
(listed in Table 2). For each anion-disordered bulk model, half of the Cl
atoms at the $4a$ sites and half of the S atoms at the $4c$ sites were
randomly selected and exchanged.

###  II.3 Self-diffusion coefficient calculation

The self-diffusion coefficients were computed from the mean square
displacements (MSDs) of Li atoms in MD simulations. The MSD at simulation time
$t$ and temperature $T$ reads,

$\textrm{MSD}(t,T)=\langle|\boldsymbol{r}(t)-\boldsymbol{r}_{0}|^{2}\rangle_{T},$
(1)

where $\boldsymbol{r}(t)$ is the position vector of a Li atom at time $t$,
$\boldsymbol{r}_{0}$ is its position vector at the reference time $t=0$, and
$\langle\cdots\rangle_{T}$ denotes the average at $T$ over all Li atoms in the
simulation cell. Generally, the MSD is related to the simulation time $t$ by

$\textrm{MSD}(t,T)\propto t^{\alpha},$ (2)

where the parameter $\alpha$ captures the type of diffusion regime [38]. For
$\alpha\neq 1$, diffusion is considered to be anomalous. Specifically, the
$\alpha<1$ regime is referred to as subdiffusion and the $\alpha>1$ regime as
superdiffusion [38]. For $\alpha=1$, i.e., when the MSD depends linearly on
$t$, a normal diffusion regime is identified, and the diffusion coefficient
$D$ can be calculated by the Einstein equation [39],

$D(T)=\frac{1}{2n}\lim_{t\to\infty}\left(\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{MSD}(t,T)\right),$
(3)

where $n$ is the number of dimensions. In simulations, the $t\to\infty$ limit
in Eq. (3) cannot be reached and has to be replaced by a finite period. In the
present work, the MSD from the simulation time interval
$35\text{\,}\mathrm{ns}$ was used to fit based on Eq. (2) and to determine
$\alpha$. Considering statistics, normal diffusion was assumed when
$\alpha\in[0.9,1.1]$.

Further, within the normal diffusion regime, the diffusion coefficient
typically follows an Arrhenius relation,

$D(T)=D_{0}\exp\left(-\frac{E_{\textrm{a}}}{k_{\textrm{B}}T}\right),$ (4)

where $D_{0}$ is a constant, $E_{\mathrm{a}}$ is the activation energy for
diffusion, and $k_{\textrm{B}}$ is the Boltzmann constant. In the present
study, the bulk and GB diffusion coefficients of Li6PS5Cl in the normal
diffusion regime were assumed to obey the Arrhenius relation (see Sec. III.4
for the confirmation of this assumption for the bulk).

Previous studies extracted GB self-diffusion coefficients by considering the
MSD only within the rather small GB region (typical widths of about
$10\text{\,}\mathrm{\text{Å}}$) [40, 41, 42]. However, the limited number of
jumping events captured in the corresponding MD simulations can lead to
significant uncertainties in the calculated GB diffusion coefficients.
Therefore, we applied a different approach to improve the accuracy of the GB
diffusion coefficients. Specifically, we first computed the _effective_
diffusion coefficient $D^{\textrm{eff}}$ for a bicrystal supercell made of two
symmetrically equivalent GBs and the bulk separating them (see Sec. II.1). To
extract from $D^{\textrm{eff}}$ the GB diffusion coefficient,
$D^{\textrm{GB}}$, we assumed that $D^{\textrm{eff}}$ is linearly composed
from $D^{\textrm{GB}}$ and a bulk diffusion contribution, $D^{\textrm{bulk}}$.
The procedure was applied separately to the GB diffusion component
perpendicular to the GB plane, $D^{\textrm{GB}}_{\perp}$, and the component
within the GB plane, $D^{\textrm{GB}}_{\parallel}$. To separate
$D^{\textrm{GB}}$ into $D^{\textrm{GB}}_{\perp}$ and
$D^{\textrm{GB}}_{\parallel}$, the vectors $\boldsymbol{r}$ entering Eq. (1)
were projected onto the $z$ direction or the $x$-$y$ plane, yielding the
respective projected diffusion coefficients (subscripts $\perp$ and
$\parallel$, respectively). For the diffusion direction within the GB plane,
the so-called Hart equation [43, 44],

$D^{\textrm{eff}}_{\parallel}(T)=(1-\tau)D^{\textrm{bulk}}_{\parallel}(T)+\tau
D^{\textrm{GB}}_{\parallel}(T),\quad n=2,$ (5)

where $\tau$ is the fraction of time spent by the diffusing atoms in the GB
region [45] and $n$ indicates the dimension as used in Eq. (3), manifests the
linear decomposition of the effective bicrystal diffusion coefficient. For the
diffusion direction perpendicular to the GB plane, long-range diffusion of Li
atoms crossing the bulk and the GB area may become relevant. However, the
present MD simulations indicate that most Li atoms diffuse around a local
area, i.e., there is no significant long-range mass transport of Li (see Sec.
III.3). Thus, approximately, the Hart equation can also be used to extract the
perpendicular GB diffusion coefficient,

$D^{\textrm{eff}}_{\perp}(T)\approx(1-\tau)D^{\textrm{bulk}}_{\perp}(T)+\tau
D^{\textrm{GB}}_{\perp}(T),\quad n=1.$ (6)

For consistency, the projected diffusion coefficients for the bulk,
$D^{\textrm{bulk}}_{\parallel}(T)$ and $D^{\textrm{bulk}}_{\perp}(T)$, were
obtained from the anion-ordered bulk models with supercell geometries
corresponding to the respective GB models (Table 2).

The linear mixing parameter $\tau$ can be approximated by the volume fraction
of the GB region in the GB model, i.e.,

$\tau\approx\frac{\delta}{l},$ (7)

where $\delta$ is the GB width, and $l$ is the distance between the GBs.
Equation (7) applies to a steady state with negligible segregation (see Sec.
III.2) and when GBs are parallel with each other in the periodic simulation
model, as in the present case (see Sec. II.1). The GB width $\delta$ was
determined based on the GB profiles (see Sec. III.4).

###  II.4 Computational details

Ab initio simulations were carried out under the density-functional theory
(DFT) framework using the projector augmented wave method [46] and the
generalized gradient approximation in the PBE parametrization [47], as
implemented in vasp [48, 49, 50]. The $1s^{2}2s^{1}$, $3s^{2}3p^{3}$,
$3s^{2}3p^{4}$, and $3s^{2}3p^{5}$ electrons were treated as valence electrons
for Li, P, S, and Cl, respectively. The plane-wave cutoff was set to
$500\text{\,}\mathrm{eV}$. For calculations of the conventional unit cell
[Fig. 1(a)], the reciprocal space was sampled using the tetrahedron method
with Blöchl corrections [51] and a $\Gamma$-centered $2\times 2\times 2$
$\boldsymbol{k}$-point mesh
($416\text{\,}\mathrm{k}\mathrm{p}\cdot\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$).
For calculations of the GB structures, the reciprocal space was sampled at the
$\Gamma$-point using the Gaussian smearing with a width of
$0.03\text{\,}\mathrm{eV}$. For structural optimization, the energy and the
maximum residual force were converged to better than
${10}^{-5}\text{\,}\mathrm{eV}$ per simulation cell and
${10}^{-2}\text{\,}\mathrm{e}\mathrm{V}\,\mathrm{\text{Å}}\textsuperscript{\textminus
1}$, respectively.

Large-scale MD simulations were performed using the validated MTP (see Sec.
III.1) within lammps [52] for the models listed in Table 2. The MD time step
was set to $2\text{\,}\mathrm{fs}$. To fully relax the structures, the so-
called annealing-and-quenching approach (a$+$q), which was used previously for
polycrystal relaxation in Ref. [53], was employed in the present study, i.e.,
the structures were first equilibrated at $600\text{\,}\mathrm{K}$ for
$0.2\text{\,}\mathrm{ns}$ and then cooled to $1\text{\,}\mathrm{K}$ within
$1\text{\,}\mathrm{ns}$. Thermal vibrations at the annealing stage enable Li
diffusion, and the slow quenching stage ensures that the system is maintained
close to a thermodynamically favorable state. The a$+$q approach was performed
in the NPT ensemble. Lastly, geometry optimization was performed using the
conjugate gradient algorithm implemented in lammps. Only the lattice variation
along the $z$ direction was allowed during relaxation. The GB formation energy
$\gamma$ at $0\text{\,}\mathrm{K}$ was calculated by,

$\gamma=\frac{E_{\textrm{GB}}-E_{\textrm{bulk}}}{A_{\textrm{tot}}},$ (8)

where $E_{\textrm{GB}}$ is the total energy of the respective GB models (Table
2) after relaxation, and $E_{\textrm{bulk}}$ is the total energy of the
reference anion-ordered bulk model after relaxation. Further,
$A_{\textrm{tot}}$ is the total GB area in the simulation cell, which is twice
as large as the $x$–$y$ cross-sectional area because of two GB planes in the
present simulation cells. For consistency, anion-ordered bulk models with a
supercell geometry corresponding to the GB models were used in Eq. (8). For
the diffusion investigations, the structures were first relaxed and
equilibrated using the Nosé–Hoover thermostat for $0.2\text{\,}\mathrm{ns}$ at
the target temperature. The NPT ensemble with a temperature rescaling every
100 time-steps at zero pressure was used. Next, to avoid any spurious effects
of the thermostat on diffusion, a microcanonical (NVE) ensemble was used, and
the system was sampled for $5\text{\,}\mathrm{ns}$. The GB profiles were
created by binning Li atoms into slices of a width of about
$5\text{\,}\mathrm{\text{Å}}$ in the $z$ direction of the simulation cell and
averaging over the entire NVE MD sampling period.

##  III Results

###  III.1 Training scheme and validation of MTP

A special training and active learning scheme for the MTP, as summarized in
Fig. 2 and outlined in the following, is proposed in the present study. The
initialization involves ab initio MD simulations for the anion-ordered bulk
structure. The trajectories obtained from these initial runs are used as the
basis for the training set. Next, the first pre-training and active learning
cycle is performed for the target training structure (e.g., a GB structure)
for an MTP of level 4. The “level” indicates the quality, i.e., the number of
fitting parameters, in the MTP [54]. The new configurations obtained during
standard active learning are labeled by DFT calculations (i.e., energies,
forces, and stresses are computed) [55] and then added to and accumulated in
the training set. In the following steps, MTPs of higher levels (6, 8, 10, and
so on) are subsequently used for pre-training and standard active learning.
Configurations obtained after each pre-training or standard active learning
cycle are labeled and added to the training set. The accuracy of the trained
MTPs at different levels is estimated by the fitting errors with respect to
DFT. Finally, when a high-accuracy MTP is obtained (as quantified below), the
active-learning scheme finishes. Further details on the proposed scheme are
given in Appendix A.

Figure 2: Proposed quality-level-based active learning scheme for the MTP.
After initialization, pre-training and standard active learning are performed
iteratively for MTPs with increasing levels, and the generated configurations
are accumulated in the training set. When the targeted level $M$ is reached,
the final MTP with high accuracy is outputted for subsequent simulations.

In the first few cycles of the proposed active learning scheme, a relatively
large number of new configurations are being sampled and added to the training
set. The corresponding information (energies, forces, stresses) is required to
provide a larger-scale estimate of the potential energy surface of the target
structure. Consequently, multiple iterations within each standard active
learning cycle of the initial MTPs are required. It is, therefore, beneficial
to use low-level MTPs at the beginning of the active learning scheme since
they can be generally constructed within a smaller number of iterations than
high-level MTPs due to their fewer fitting parameters. The low-level MTPs are
also more robust in phase-space extrapolation, requiring less computational
cost for refitting. Starting the whole active learning scheme with low-level
MTPs thus generates, in a computationally efficient manner, an ample training
set for the high-level MTPs at the later stages of the active learning scheme.
In this way, the number of standard active learning iterations is kept small
for the computationally expensive high-level MTPs. Overall, the proposed
training and active learning scheme systematizes and accelerates the
construction of precise MTPs for chemically and structurally complex
materials.

In the present study, various MTPs were fitted from level 4 up to level 20 for
the anion-ordered and the anion-disordered bulk structures. Based on the bulk
results, single MTPs were fitted for each GB from level 4 up to level 18.
Figures 3(a) and (b) show the validation root-mean-square errors (RMSEs) in
energies and forces for the anion-ordered and anion-disordered bulk,
respectively. Additionally, Fig. 3(c) shows the dependence of the target
quantity, i.e., the diffusion coefficient, on the MTP level. For the anion-
ordered bulk structure, the RMSE in energy is less than
$10\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$
for MTPs from level 14 on. The RMSE in force slowly decreases with the MTP
level and reaches a small final value of about
$0.05\text{\,}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{Å}}$ for level 20.
Larger RMSEs in both energy and force are observed for the anion-disordered
bulk structure, likely due to its complex anionic arrangement. With MTPs from
level 10 on, the calculated diffusion coefficients for the ordered and the
disordered bulk structures converge to final values of about
$20\text{\times}{10}^{-7}180\text{\times}{10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$,
respectively. The ordered bulk structure also shows less variation among
multiple MTPs for both the RMSEs and the diffusion coefficients.

Figure 3: Validation root-mean-square error (RMSE) in (a) energy and (b) force
of the MTPs obtained for the anion-ordered and the anion-disordered bulk
structures from level 4 to 20 with the proposed scheme (Fig. 2). (c) Diffusion
coefficients $D$ at $600\text{\,}\mathrm{K}$ calculated with the obtained
MTPs. Each line represents the result of an independently fitted MTP. Multiple
MTPs were fitted for each structure to show the statistical variation.
Simulations were performed with supercells with $3\,328$ atoms [equivalent to
a $4\times 4\times 4$ expansion of the unit cell shown in Fig. 1(a)].

Variation in the calculated diffusion coefficients may originate from the
following sources: 1. Limited diffusion sampling time; 2. Different local
minima within the MTP parameter space; 3. Different accumulated training sets.
To systematically investigate these sources, the diffusion coefficients of the
anion-ordered bulk structure were calculated at $600\text{\,}\mathrm{K}$ for
three different sets of simulations:

  1. A.

The complete active learning scheme was run once, thereby generating one MTP
at each level (from 4 to 20). For each MTP, five independent runs for
$5\text{\,}\mathrm{ns}$ were performed, and the corresponding standard
deviation for the diffusion coefficient was determined. The standard deviation
for set A represents the contribution to the variation due to the finite
sampling time of the diffusion coefficient.

  2. B.

The final training set obtained after the active learning in A was used to fit
ten MTPs for each level. Different MTPs (in terms of the basis set
coefficients) can be obtained even for the same training set and for the same
MTP level because of multiple local minima in the MTP parameter space visited
due to randomized initial starting conditions. At each level, the 10 MTPs were
used to run diffusion coefficient simulations for $5\text{\,}\mathrm{ns}$. The
corresponding standard deviation in the diffusion coefficient was determined
for each level. The standard deviation for set B contains the contribution
from the finite sampling time (set A) and additionally the variation from the
different MTP parameters.

  3. C.

The whole active learning scheme was run ten times, thereby generating 10
different MTPs at each level. These MTPs vary not only in the MTP parameters
but, importantly, also in the training set used for the fitting. Thus, set C
includes, in addition to the previous variation contributions (as for set B),
the variation due to different training sets.

The standard deviations of the diffusion coefficients obtained for these three
sets of simulations are shown in Fig. 4(a). The results allow us to draw a
very important conclusion. All curves decrease with an increasing MTP level
and reach very small values for the highest investigated levels. This means
that the final diffusion-coefficient values for the highest MTP levels do not
depend (beyond the remaining statistical variation) on the initial starting
conditions of the active learning scheme. Thus, it is possible to
systematically increase the precision in the predicted diffusion coefficient
values of the MTPs. However, one should be aware of the accompanying strongly
increasing computational cost with MTP level as shown in Fig. 4(b). To balance
the accuracy and simulation cost, MTPs of level 18 with 807 fitting parameters
were used in the present study.

Figure 4: (a) Standard deviation of the diffusion coefficients $D$ of the
anion-ordered bulk structure at $600\text{\,}\mathrm{K}$ for the three
different sets of simulations A, B, and C as described in the main text. As a
reference, the diffusion coefficient of the anion-ordered bulk structure
calculated with a level-20 MTP at $600\text{\,}\mathrm{K}$ is about
$20\text{\times}{10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
[Fig. 3(c)]. (b) CPU time required to calculate $D$ for the ordered bulk
structure. An exponential function was fitted to show the CPU time and the MTP
level relation. To balance the accuracy and simulation cost, MTPs of level 18
(Table 3) were used in the subsequent simulations. Simulations were performed
with supercells with $3\,328$ atoms.

Table 3 lists the information on the final level-18 MTPs trained for different
target structures. Larger training RMSEs in both energies and forces are
observed for the GBs compared to the bulk structures (more than
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$
in energy). Despite applying a different anion disorder in constructing the
validation set for the MTP of the disordered bulk structure, the validation
RMSE remains consistent with the training RMSE. Details of the construction of
the validation set are given in Appendix A.

Table 3: Information of the final MTPs. The number of atoms in the simulation cell of the training structures is denoted by $N_{\textrm{at}}$, and the number of configurations in the final accumulated training set is denoted by $N_{\textrm{conf}}$. Different MTPs were independently trained to different structures according to the proposed scheme up to level 18 (i.e., $M=18$ in Fig. 2; cf. Fig. 4) with 807 fitting parameters. Training and validation root-mean-square error (RMSE) is shown for the bulk structures, and training RMSE is shown for the GB structures.  Training structure | $N_{\textrm{at}}$ | $N_{\textrm{conf}}$ | Training RMSE | Validation RMSE
---|---|---|---|---
|  |  | Energy (meV/atom) | Force (eV/Å) | Energy (meV/atom) | Force (eV/Å)
anion-ordered bulk | 104 | 2219 | 4.7 | 0.132 | 3.1 | 0.056
anion-disordered bulk | 208 | 2210 | 5.1 | 0.153 | 6.0 | 0.095
GB $\Sigma 3(1\bar{1}2)[110]$ | 312 | 2350 | 6.6 | 0.164 | n/a | n/a
GB $\Sigma 3(\bar{1}11)[110]$ | 156 | 2293 | 6.9 | 0.157 | n/a | n/a
GB $\Sigma 5(001)[001]$ | 260 | 2241 | 6.1 | 0.157 | n/a | n/a

Figure 5 shows the radial distribution functions of the anion-ordered bulk
structure in the conventional unit cell [Fig. 1(a)] and the $\Sigma
5(001)[001]$ GB obtained from MD simulations at $600\text{\,}\mathrm{K}$. The
accuracy of the MTPs is validated by the small difference in the radial
distribution functions obtained by MTP and DFT.

Figure 5: Comparison of radial distribution functions obtained by DFT and MTP
MD at $600\text{\,}\mathrm{K}$ for the anion-ordered bulk and the $\Sigma
5(001)[001]$ GB averaged over 4000 and 2000 MD steps, respectively. Analysis
was performed using ovito [56].

###  III.2 Structures and energies

Table 4 shows structural and energetic information on the anion-ordered bulk
Li6PS5Cl at $0\text{\,}\mathrm{K}$ as obtained with various optimization
approaches. For the optimization with symmetry constraint (i.e., preserving
the space group $F\bar{4}3m$), DFT and MTP show similar lattice constants of
about $10.25\text{\,}\mathrm{\text{Å}}$. The symmetry-constrained structure is
dynamically unstable at $0\text{\,}\mathrm{K}$ as confirmed by additional
phonon calculations. The related double-well potentials of two imaginary modes
were analyzed previously in Ref. [34]. The dynamical instability is also
consistent with previous experimental [30, 31] and computational [33, 32]
studies reporting that the high symmetry Li6PS5Cl structure [Fig. 1(a)] is not
the most stable phase at low temperatures and $0\text{\,}\mathrm{K}$,
respectively.

Table 4: Structural and energetic data for the anion-ordered Li6PS5Cl bulk at $0\text{\,}\mathrm{K}$ based on different methods.444A lattice constant of $9.818\text{\,}\mathrm{\text{Å}}$ was obtained in experiments at $150\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{a}\mathrm{t}\mathrm{m}$ for a sample with about $56.2\text{\,}\mathrm{\char 37\relax}$ of Cl/S-anion disorder showing the space group $F\bar{4}3m$ [9]. Note that this experimental value is not directly comparable with the simulation results given in the table because of the differences in the space group, temperature, and anion ordering. The total number of atoms in the simulation cell is denoted by $N_{\textrm{at}}$. The structure before optimization [Fig. 1(a)] has a cubic symmetry with the space group $F\bar{4}3m$. The cubic shape of the simulation cell is constrained during the a$+$q optimization, and $a$ refers to the cubic lattice constant. The energy $\Delta E$ is the energy of the optimized structure referenced with respect to the symmetry-constrained structure. Source | $N_{\mathrm{at}}$ | Method |  $a$ (Å) | Space group | Stability555The Stability column indicates the dynamical stability as determined by the presence or absence of imaginary phonons at $0\text{\,}\mathrm{K}$. |  $\Delta E$ (meV/atom) | Reference
---|---|---|---|---|---|---|---
DFT | $52$ | CG666The conjugate gradient (CG) algorithm implemented in vasp was used.  with sym. constraint | $10.246$ | $F\bar{4}3m$ | unstable | $0.0$ | this work
MTP | $52$ | BFGS777The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm implemented in ase [57] was used.  with sym. constraint | $10.264$ | $F\bar{4}3m$ | unstable | $0.0$ | this work
MTP | $52$ | a$+$q888The annealing-and-quenching (a$+$q) approach described in Sec. II.4 was used. | $9.950$ | $P1$ | stable | $-55.3$ | this work
DFT | $52$ | MTP a$+$q $\rightarrow$ DFT CG999The optimized structure with MTP using the a$+$q44footnotemark: 4 method was further optimized with DFT using the CG22footnotemark: 2 method. |  $9.940$101010The cubic cell shape was changed after this optimization, so the effective cubic lattice constant is shown. | $P1$ | n/a | $-58.9$ | this work
DFT | $13$ | USPEX111111The universal structure predictor: evolutionary xtallography (USPEX) described in Refs. [58, 32] was used. |  $10.209$77footnotemark: 7 | $P1$ | stable | $-50.0$ | [32]
MTP | $17\,836$ | a$+$q | $9.954$ | $P1$ | n/a | $-57.7$ | this work

Due to the complicated potential energy surface of Li6PS5Cl with multiple
local minima, it is difficult to find the fully optimized structure based on
direct optimization methods (e.g., the conjugate gradient algorithm). To
enhance the search for the global energy minimum, an MD-based
annealing$+$quenching (a$+$q) optimization approach with volume relaxation was
carried out with the MTP (see Sec. II.4). The MTP-optimized structure was
further relaxed with DFT utilizing the conjugate gradient algorithm to
validate the MTP result. The resulting structure has an energy of
$-58.9\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$
with respect to the $F\bar{4}3m$ high-symmetry state, which is energetically
more favorable than the structure obtained previously based on optimization
with DFT and the USPEX method
($-50.0\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$)
[32]. The corresponding MTP energy of our optimized structure is
$-55.3\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$
and thus in agreement with the corresponding DFT value with an error of
$+3.6\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$,
consistent with the validation RMSE in energy [Fig. 3(a) and Table 3]. For a
larger simulation cell with $17\,836$ atoms (the anion-ordered bulk model
listed in Table 2), the optimization with the MTP yields an energy that is
further reduced by
$2.4\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$.
This result indicates that there are favorable arrangements of the Li atoms
that cannot be represented with the small simulation cell. Tests show that
repeating relaxations using the a$+$q approach may lead to different
arrangements of Li atoms but have a small impact on the total energy of the
relaxed structure (less than
$0.1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$).

For the GB models, an additional step, the so-called $\gamma$-surface search
[59], should be considered to optimize the interface between the two grains
before using the a$+$q approach. Specifically, in the $\gamma$-surface search,
different rigid shifts of the two grains in directions parallel to the GB
plane ($x$-$y$) are investigated in order to probe different initial GB
structures for further optimization. For each GB model, we investigated a mesh
of $20\times 20$ points on the respective $\gamma$-surface, by performing
direct optimization for each such point. The structure of the lowest-energy
point in the $\gamma$-surface was selected for further relaxation with the
a$+$q approach. Additionally, the original GB structures without a shift were
also used for relaxation with the a$+$q approach. Full relaxation of the
structures with or without the $\gamma$-surface search leads to similar final
optimized structures. This finding indicates that the shifting of two grains
is implicitly included in the a$+$q optimization process, at least for the
three GBs considered in the present study.

The atomic structures of the three GB models after optimization are shown in
Fig. 6. To better emphasize the characteristic features, Li atoms are not
visualized, and only the atoms close to the GB plane are shown. The dimensions
of the optimized structures, and GB energies [$\gamma$, Eq. (8)] calculated at
$0\text{\,}\mathrm{K}$ are shown in Table 2. For $\Sigma 3(1\bar{1}2)[110]$
[Figs. 6(a) and (b)], voids (as large as $6\text{\,}\mathrm{\text{Å}}$ in
diameter) form at the GB plane. The formation of these voids requires the
reordering of a range of atoms, which in turn requires enough thermal energy
that is available through the a$+$q process. In contrast, direct optimization
could not reproduce the voids, emphasizing again the importance of the a$+$q
process. For the $\Sigma 3(\bar{1}11)[110]$ GB [Figs. 6(c) and (d)], for which
the cage structure is preserved during the GB construction, only a small
change of the cage arrangement is visible at the GB plane after optimization.
In particular, no void formation is observed. Correspondingly, $\Sigma
3(\bar{1}11)[110]$ shows a GB energy of only
$7.12\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{2}$,
which is the smallest among the three investigated GBs. Such a small value
indicates that the GB structure does not differ much from the bulk structure.
For the $\Sigma 5(001)[001]$ GB [Figs. 6(e) and (f)], for which the cage
structure is preserved and the GB plane is densely filled with atoms, yet
other results are observed than for the two $\Sigma 3$ GBs. For the $\Sigma 5$
GB, optimization leads to an amorphous-like area at and near the GB plane. The
amorphous area extends from the Cl atoms at the GB plane to the next layer of
the Cl atoms away from the GB plane. The atoms inside the affected area appear
to be structurally disordered. Since the bonds at the GB plane are severely
modified compared to the bulk, the largest GB energy
($18.78\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{2}$)
is found for the $\Sigma 5(001)[001]$ GB. Results show that the excess volumes
per unit GB area of all GBs are small, with values less than
$1\text{\,}\mathrm{\text{Å}}$.

Figure 6: Atomic structures of the three GB models (Table 2) after
optimization. Each GB model is visualized in two views (top and bottom row),
with the GB normal lying in the paper plane for both views. The GB plane is
indicated by the vertical black line in the middle of each structure. For
visibility, Li atoms and cages are not visualized, and only the atoms and PS4
units close to the GB plane (within about $\pm$10\text{\,}\mathrm{\text{Å}}$$)
are shown. The smallest repeating GB structural units are highlighted for the
$\Sigma 3$ GBs.

Figure 7 shows the atomic energy and the atomic concentration profiles of the
Li atoms for the three GB models. To emphasize the impact of the GBs, the
energy profiles are referenced with respect to their bulk counterparts, and
the concentration profiles are normalized with the average Li concentration of
the GB simulation cells. According to the energy profiles (upper panel), the
atomic energy of the Li atoms increases in the GB region for the $\Sigma
3(1\bar{1}2)[110]$ (up to
$0.2\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$)
and the $\Sigma 5(001)[001]$ (up to
$0.1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$)
GBs. For the $\Sigma 3(\bar{1}11)[110]$ GB, the change in atomic energy of Li
in the GB region is negative (about
$-0.1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$).
For the $\Sigma 3(1\bar{1}2)[110]$ and $\Sigma 5(001)[001]$ GBs, the impact of
the GB structure on the Li atomic energy is slightly reduced as temperature
increases. As regards the concentration profiles (lower panel), we observe
layers with increased and decreased Li concentrations around the GB plane.
Especially for $\Sigma 5(001)[001]$, Li atoms tend to segregate from the
center to the edge of the amorphous-like GB region. With higher temperatures,
smaller Li concentration fluctuations are found.

The atomic energy profiles in Fig. 7 were utilized to determine the GB widths
as needed for the GB diffusion coefficient calculations (see Sec. III.4).
Specifically, a Gaussian function was fitted to the peaks in the Li atomic
energy, and the GB width was defined to be six times the standard deviation of
the fitted Gaussian function. A similar approach to determine the GB width was
used in Ref. [60]. Since the temperature dependence of the standard deviation
of the fitted Gaussian function is very small, the GB widths obtained from the
$600\text{\,}\mathrm{K}$ profiles (shown in Table 2) were used throughout the
GB diffusion coefficient calculations. The estimated GB widths are similar
between $\Sigma 3(1\bar{1}2)[110]$ ($28.8\text{\,}\mathrm{\text{Å}}$) and
$\Sigma 5(001)[001]$ ($26.5\text{\,}\mathrm{\text{Å}}$) and smaller for
$\Sigma 3(\bar{1}11)[110]$ ($16.7\text{\,}\mathrm{\text{Å}}$), consistent with
visual inspection of the Li diffusion trajectories (see Sec. III.3). The
volume fractions [$\tau$ in Eq. (7)] of the three GBs in the GB models are
between $2539\text{\,}\mathrm{\char 37\relax}$.

The determined GB regions for the three GBs are indicated in the atomic
concentration profiles in Fig. 7. Analysis shows that the Li concentration
averaged within the GB region is comparable to the average Li concentration of
the reference bulk structure. It indicates that the segregation of Li in the
GB region is negligible, supporting the approximation in Eq. (7). Li
segregation in the present simulations may be hindered by the restriction of
all simulation cells to the stoichiometric composition of Li6PS5Cl (see Sec.
II.1).

Figure 7: Temperature-dependent GB profiles of the three GB models (Table 2)
averaged over $5\text{\,}\mathrm{ns}$ of MD simulation time. The upper panel
shows the atomic energy of Li atoms referenced with respect to the bulk. A
Gaussian fit (gray line) of the data from $600\text{\,}\mathrm{K}$ was used
for GB width analysis. The lower panel shows the atomic concentration of Li
atoms normalized with the average Li concentration of the GB simulation cells.
The kernel-density estimation was used for smoothing. The gray area indicates
the GB regions, which were defined to be six times the standard deviations
(stand. dev.) of the Gaussian fits in the upper panel. The Li concentration
averaged over the GB region was used to confirm that the segregation of Li in
the GB region is negligible.

###  III.3 Bulk and GB diffusion

Trajectories of all Li atoms obtained from $5\text{\,}\mathrm{n}\mathrm{s}$ MD
runs at $300600\text{\,}\mathrm{K}$ for the anion-ordered bulk structure are
shown in Figs. 8(a) and (b), respectively. The trajectories of a few selected
Li atoms are highlighted in blue. At $300\text{\,}\mathrm{K}$, the Li
trajectories are mostly localized within the cages, and only very few inter-
cage jumps occur during the simulation time. At $600\text{\,}\mathrm{K}$, the
Li trajectories are distributed much more homogeneously in the simulation cell
due to the enhanced inter-cage (long-range) diffusion. Clearly, temperature
strongly affects Li diffusion, especially inter-cage diffusion. It is worth
emphasizing that the macroscopic conductivity measured on the experimental
scale mainly results from inter-cage diffusion.

Figure 8: Trajectories of all Li atoms in the ordered bulk structure (Table 2)
at (a) $300\text{\,}\mathrm{K}$ and (b) $600\text{\,}\mathrm{K}$ after
$5\text{\,}\mathrm{ns}$ of MD simulations. The trajectories of a few selected
Li atoms are highlighted in blue, and the trajectories of the other Li atoms
are shown in gray. Magnified views of one of the highlighted Li trajectories
at (c) $300\text{\,}\mathrm{K}$ and (d) $600\text{\,}\mathrm{K}$ relative to
the initial position are also shown. The color mapping indicates the
probability density smoothed by kernel-density estimation. The elapsed time
between the snapshots is $1\text{\,}\mathrm{ps}$. (e) Inter-cage jumping rate
and (f) cage occupancy of Li at different temperatures. The color mapping
indicates the percentage of Li atoms that exhibit the corresponding jumping
rate or cage occupancy. The result of each temperature is averaged over
$5\text{\,}\mathrm{ns}$ of MD simulation. Trajectories in (a) and (b) are
visualized by ovito [56].

Figures 8(c) and (d) are zoom-ins of the trajectories of a single Li atom at
$300600\text{\,}\mathrm{K}$, respectively. At $300\text{\,}\mathrm{K}$, the
selected Li atom shows frequent _intra_ -cage diffusion visiting multiple
interstitial sites of the same cage within the
$5\text{\,}\mathrm{n}\mathrm{s}$ of simulation time but no _inter_ -cage
diffusion. At $600\text{\,}\mathrm{K}$, in contrast, significant intra- and
inter-cage diffusion is observed for the selected Li atom. After jumping to a
neighboring cage, the Li atom typically resides in the new cage for some time,
during which various interstitial sites are visited via intra-cage diffusion.
Occasionally, the residence time is shorter, and the Li atom diffuses quickly
to another cage. Short residence times in the cages are more frequently
observed at higher temperatures, and they further enhance long-range diffusion
and conductivity.

To quantify the frequency of inter-cage jumps, Fig. 8(e) shows the inter-cage
jump rates per Li atom and per nanosecond as a function of temperature. At the
lowest investigated temperature ($200\text{\,}\mathrm{K}$), almost no inter-
cage jumps are observed within the $5\text{\,}\mathrm{n}\mathrm{s}$ simulation
time. The average inter-cage jumping rate increases strongly with increasing
temperature. Consequently, the macroscopic conductivity likewise increases
strongly with temperature.

In the anion-ordered bulk structure, each cage consists of six Li atoms around
one S atom. When inter-cage jumps of Li occur, it can be anticipated that
there will be an imbalanced distribution of Li in the cages, i.e., single
cages may show a Li occupancy of more or less than six. Figure 8(f) shows the
cage occupancy averaged over the entire MD simulation time at different
temperatures. The color mapping in the plot indicates the percentage of the
cages in the simulation cell that exhibit the corresponding occupancy. As
expected, when very few inter-cage jumps are observed, e.g., at
$200\text{\,}\mathrm{K}$, almost all cages are occupied by six Li atoms, i.e.,
neutral occupation. With increasing temperature and more inter-cage jumps
[Fig. 8(e)], the cage occupancy distribution broadens, i.e., cages with more
or less than six Li atoms are observed. Taking the highest investigated
temperature of $800\text{\,}\mathrm{K}$, about 18% of cages have an imbalance
of $\pm 1$ Li atoms (i.e., cages with five or seven Li atoms). The majority
(82%) of the cages are still occupied with six Li atoms.

Figure 9 shows an analysis of Li diffusion for the three GB models, with a
focus on $300\text{\,}\mathrm{K}$. Similar to the bulk analysis [Figs. 8(a)
and (b)], trajectories of all Li atoms obtained from
$5\text{\,}\mathrm{n}\mathrm{s}$ MD simulations for the three GB models are
shown in Figs. 9(a)–(c). Trajectories far from the GB regions are similar to
the bulk results discussed above. In the GB regions, differences can be seen
depending on the GB type. For $\Sigma 3(1\bar{1}2)[110]$ [Fig. 9(a)], broken
cages are connected by Li trajectories across the GB plane, leading to
multiple complex diffusion paths. Additionally, small voids (regions not
covered by Li trajectories) in diameter of about $5\text{\,}\mathrm{\text{Å}}$
are found at the GB plane. These voids result from the intersection of rows of
cages from the two grains, oriented with a rotational angle of
$70.53\text{\,}\mathrm{\SIUnitSymbolDegree}$ (the tilt angle). The impact on
Li diffusion is substantially different when the GB planes do not break any
cages: For $\Sigma 3(\bar{1}11)[110]$ [Fig. 9(b)], the GB plane is located
between two rows of cages from each grain. Most Li atoms show only intra-cage
diffusion, even in the GB region. Only a few trajectories cross the GB planes.
For $\Sigma 5(001)[001]$ [Fig. 9(d)], the analysis of the Li trajectories is
complicated by the projection plane. Nevertheless, one can clearly distinguish
different trajectories within the GB region. The trajectories appear
homogeneously distributed, with no clear voids visible.

Figure 9: Analysis of Li diffusion for the three GB models (Table 2) at
$300\text{\,}\mathrm{K}$ and after $5\text{\,}\mathrm{ns}$ MD simulations.
(a), (b), and (c) Trajectories of all Li atoms in the simulation cells. The
trajectories of a few selected Li atoms in the GB region are highlighted in
blue, and the trajectories of the other Li atoms are shown in gray. For each
GB model, trajectories are visualized in two views. GB planes are located in
the middle and at the edge of the GB simulation cells in the $z$-direction.
(d) Distribution of the Li displacements for the GB models (gray-shaded) in
comparison with that for the anion-ordered bulk model (black-edged).

The trajectories of a few selected Li atoms in the GB regions are highlighted
in Figs. 9(a), (b), and (c). Relatively long trajectories are found, which
indicates that the GB structure increases the probability of long-range
diffusion compared to the anion-ordered bulk. For $\Sigma 3(1\bar{1}2)[110]$,
the long-range trajectories connect the broken cages at the GB plane.
Trajectories combined with both intra- and inter-cage jumps are observed for
the $\Sigma 3(\bar{1}11)[110]$ GB. Similar trajectories are also found for the
anion-ordered bulk but at a higher temperature, e.g., $600\text{\,}\mathrm{K}$
[Figs. 8(c) and (d)]. For $\Sigma 5(001)[001]$, trajectories reveal straighter
diffusion paths for Li crossing the GB plane, suggesting a weaker trapping
effect of the cages for Li. Trajectories within the GB plane [right column of
Figs. 9(a)–(c)] suggest negligible in-plane anisotropy in the GB diffusion for
the three GBs.

In Fig. 9(d), distributions of Li according to displacements from the GB
simulations are compared to those derived from simulations of the anion-
ordered bulk structure. For the bulk structure, most Li atoms show
displacements up to about $6\text{\,}\mathrm{\text{Å}}$. Considering that the
maximum atomic distance within one cage at $0\text{\,}\mathrm{K}$ is about
$4.6\text{\,}\mathrm{\text{Å}}$, it can be concluded that most Li atoms at
$300\text{\,}\mathrm{K}$ are restricted to intra-cage diffusion, consistent
with the previous analysis (Fig. 8). At the same temperature and simulation
time, longer Li displacements and a higher probability of long-distance
displacements are observed for the GB models compared to the bulk.
Importantly, a different impact on Li mobility is observed for the different
GB types. A strong enhancement in Li diffusion is shown for $\Sigma
3(1\bar{1}2)[110]$ and $\Sigma 5(001)[001]$ GBs. A relatively small
enhancement in Li diffusion is seen for $\Sigma 3(\bar{1}11)[110]$, which
mainly results from a higher probability of Li atoms that show displacements
between 7 and $11\text{\,}\mathrm{\text{Å}}$.

###  III.4 Diffusion coefficients

Diffusion coefficients were calculated from the MSDs according to Eq. (3) for
normal diffusion. Figure 10(a) shows the MSDs of the Li atoms in the anion-
ordered bulk as a function of time and for different temperatures at intervals
of $50\text{\,}\mathrm{K}$ as a representative example. Note that a double-
logarithmic scale is used to reveal the diffusion type. At low temperatures
(up to $300\text{\,}\mathrm{K}$), a plateau in the MSD is visible after a
certain simulation time. The plateau indicates that most Li atoms in the
simulation cell exhibit only intra-cage (short-range) diffusion within the
simulation time, i.e., Li atoms are trapped in the cages. This is fully
consistent with the above analysis of the Li trajectories [Figs. 8(a) and
(b)]. At elevated temperatures ($300600\text{\,}\mathrm{K}$), the MSD curves
show a characteristic increase–plateau–increase shape. This shape indicates a
combination of intra- and inter-cage jumps (short- and long-range diffusion)
of Li atoms in the simulation cell, as also implied by the increased inter-
cage jump rates at higher temperatures shown in Fig. 8(e). At temperatures
exceeding $600\text{\,}\mathrm{K}$, the plateau in the MSD disappears because
the long-range diffusion of Li sets in already in the initial stages of the
simulation. Similar MSD curves are found for the anion-disordered bulk and the
three GB models.

Figure 10: (a) Mean square displacements (MSD) of all Li atoms in MD
simulations of the anion-ordered bulk models (Table 2) as a function of
simulation time up to $5\text{\,}\mathrm{n}\mathrm{s}$. The color of the MSD
curves indicates the simulation temperature. The slope corresponding to
$\alpha=1$ [Eq. (2), normal diffusion] is shown by dashed black lines for
reference. (b) $\alpha$ values as a function of temperature obtained by
fitting the MSD within the simulation time interval
$35\text{\,}\mathrm{n}\mathrm{s}$ according to Eq. (2). The normal diffusion
regime, defined by the criterion $\alpha\in[0.9,1.1]$ (Sce. II.3), is
indicated. Note that the MSDs from the anion-disordered bulk and the GB models
are not shown here, and the GB models include both the bulk and the GB regions
(see Sec. II.1).

A linear dependence of the MSD with a slope of $\alpha=1$ [see Eq. (2)]
identifies normal diffusion, while the respective offset corresponds to the
diffusion–coefficient value. The dotted line in Fig. 10(a) indicates a slope
of one as a visual guide. According to the MSDs, subdiffusion ($\alpha<1$) is
found for simulations at low temperatures in the whole simulation time frame
up to $5\text{\,}\mathrm{n}\mathrm{s}$. At high temperatures, subdiffusion is
also observed for short observation times as evidenced by the initial slope
below one (e.g., at $800\text{\,}\mathrm{K}$ in the time interval up to
$10\text{\,}\mathrm{ps}$). During these short observation times, there are not
enough inter–cage jumps in the simulation to trigger long-range diffusion.
These results emphasize the importance of carefully checking the slope of the
MSD curve to ensure a correct application of the Einstein equation [Eq. (3)].
Figure 10(b) shows the $\alpha$ values as a function of temperature obtained
by fitting the MSD within the simulation time interval of
$35\text{\,}\mathrm{n}\mathrm{s}$ according to Eq. (2). We define normal
diffusion by the criterion $\alpha\in[0.9,1.1]$ (see Sec. II.3). Based on this
definition, normal diffusion occurs for the anion-ordered bulk at temperatures
higher than $500\text{\,}\mathrm{K}$. For the GBs that enhance Li diffusion,
normal diffusion is already observed at lower temperatures. The anion-
disordered bulk structure with $50\text{\,}\mathrm{\char 37\relax}$ anion
disorder shows normal diffusion at the lowest temperatures, i.e., from
$350\text{\,}\mathrm{K}$ on. Diffusion coefficients were calculated for
temperatures within the defined normal–diffusion regime.

The calculated diffusion coefficients and the resulting Arrhenius fits are
shown in Fig. 11 for the anion-ordered and anion-disordered bulk structures
(open symbols and black lines). For the three GBs, only the Arrhenius fits are
shown, separated into diffusion coefficients perpendicular to and within the
GB plane. The corresponding Arrhenius fitting parameters are listed in Table
5. Clear linear Arrhenius relations are observed for the diffusion
coefficients of both the anion-ordered and the anion-disordered bulk
structures with respect to the inverse temperature. The ordered bulk structure
shows a high diffusion activation energy of
$384\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, which is lowered to
$216\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ by the random
$50\text{\,}\mathrm{\char 37\relax}$ anion disorder of the Cl and S anions.
This is in good agreement with previous experiments [17] and simulations [9]
focusing on the effect of anion disorder. For comparison, Fig. 11 includes
results of anion-ordered [12, 13, 14] and anion-disordered [18, 15] bulk
structures obtained by ab initio MD simulations from previous studies, all of
which show significant deviations compared to the values obtained in the
present work. This is likely due to smaller simulation cells
{$52\text{\,}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}\mathrm{s}$ [12, 13, 14,
18, 15] vs. more than
$16\,000\text{\,}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}\mathrm{s}$ in the
present work (Table 1)} and shorter simulation times {up to
$300\text{\,}\mathrm{p}\mathrm{s}$ [12, 13, 14, 18, 19, 15] vs.
$5\text{\,}\mathrm{ns}$ in the present work (Table 1)} due to the expensive ab
initio MD simulations utilized in the previous studies.

Figure 11: Arrhenius diagram for diffusion coefficients $D^{\mathrm{bulk}}$, $D^{\mathrm{GB}}_{\parallel}$, and $D^{\mathrm{GB}}_{\perp}$ obtained by atomistic simulations. Anion-ordered and anion-disordered ($50\text{\,}\mathrm{\char 37\relax}$ Cl/S-anion disorder) bulk results are included. The fitted Arrhenius lines for the three GBs and directions parallel ($\parallel$, solid lines) and perpendicular ($\perp$, dashed lines) to the GB plane are shown. For comparison, bulk diffusion coefficients obtained by ab initio MD [shorter sampling, up to $300\text{\,}\mathrm{p}\mathrm{s}$ vs. $5\text{\,}\mathrm{ns}$ in the present work (Table 1)] from Refs. [12, 13, 14, 15] (ordered) and [18, 15] (disordered, $50\text{\,}\mathrm{\char 37\relax}$ Cl/S-anion disorder) are shown. The experimental data marked by $+$ (Adeli  _et al._ [7]), × (Schlenker  _et al._ [8]), and $\star$ (Hanghofer  _et al._ [61]) are based on nuclear magnetic resonance measurements, while $\ast$ (Deiseroth  _et al._ [31]) are based on polarization measurements. The samples corresponding to $+$ and × showed $61.5\text{\,}\mathrm{\char 37\relax}$ (Adeli  _et al._ [7]) and $53.8\text{\,}\mathrm{\char 37\relax}$ (Schlenker  _et al._ [8]) Cl/S-anion disorder, respectively.  Table 5: Fitted Arrhenius parameters [$D_{0}$, $E_{\textrm{a}}$, see Eq. (4)] and self-diffusion coefficients at $300\text{\,}\mathrm{K}$ of bulk and pure GBs, i.e., $D^{\mathrm{bulk}}$, $D^{\mathrm{GB}}_{\parallel}$, and $D^{\mathrm{GB}}_{\perp}$ obtained in the present work. Note that $D^{\mathrm{GB}}_{\parallel}$ and $D^{\mathrm{GB}}_{\perp}$ in Eqs. (5) and (6), respectively, are evaluated with the estimated GB widths listed in Table 2 (for the impact of GB widths, see Sec. IV). The experimental results obtained by nuclear magnetic resonance measurements in previous studies [7, 8] are shown for comparison. The ionic conductivities, $\sigma$, are converted from the calculated diffusion coefficients, $D$, based on a relation extracted from the experimental data in Ref. [8] (see Appendix B). Structure | Direction |  $D_{0}$ (${10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$) |  $E_{\mathrm{a}}$ (meV) |  $D$ at $300\text{\,}\mathrm{K}$ (${10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$) |  $\sigma$ at $300\text{\,}\mathrm{K}$ ($\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}$)
---|---|---|---|---|---
anion-ordered bulk | all | $3.3\text{\times}{10}^{4}\text{\,}$ | $384$ | $0.012$ | $0.2$
anion-disordered bulk121212The structure has $50\text{\,}\mathrm{\char 37\relax}$ of Cl/S-anion disorder. | all | $9.4\text{\times}{10}^{3}\text{\,}$ | $216$ | $2.203$ | $29.8$
GB $\Sigma 3(1\bar{1}2)[110]$ | $\parallel$ | $9.4\text{\times}{10}^{3}\text{\,}$ | $274$ | $0.236$ | $3.2$
| $\perp$ | $9.0\text{\times}{10}^{3}\text{\,}$ | $275$ | $0.218$ | $2.9$
GB $\Sigma 3(\bar{1}11)[110]$ | $\parallel$ | $3.9\text{\times}{10}^{4}\text{\,}$ | $360$ | $0.035$ | $0.5$
| $\perp$ | $4.0\text{\times}{10}^{4}\text{\,}$ | $369$ | $0.026$ | $0.3$
GB $\Sigma 5(001)[001]$ | $\parallel$ | $1.1\text{\times}{10}^{4}\text{\,}$ | $270$ | $0.307$ | $4.1$
| $\perp$ | $1.1\text{\times}{10}^{4}\text{\,}$ | $280$ | $0.222$ | $3.0$
polycrystal in experiment131313The sample showed $61.5\text{\,}\mathrm{\char 37\relax}$ of Cl/S-anion disorder [7]. | n/a | n/a |  $350$ p m 10 | $0.387$ | $2.5$
polycrystal in experiment141414The sample showed $53.8\text{\,}\mathrm{\char 37\relax}$ of Cl/S-anion disorder [8]. | n/a | n/a |  $280$ p m 10 | $0.25$ | $3.4$

Enhanced diffusion is observed for all three GBs as compared to the anion-
ordered bulk
($0.026\text{\times}{10}^{-7}0.307\text{\times}{10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
for GBs vs.
$0.012\text{\times}{10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
for the anion-ordered bulk at $300\text{\,}\mathrm{K}$). The enhancement is
sensitive to the GB width (details in Sec. IV). Based on the current width
estimates [Figs. 7 and Table 2], varying increases are shown at
$300\text{\,}\mathrm{K}$: from $\Sigma 3(\bar{1}11)[110]$ (about two to three
times) to $\Sigma 3(1\bar{1}2)[110]$ or $\Sigma 5(001)[001]$ (more than an
order of magnitude). For all GBs, higher diffusivity is shown for Li diffusing
along rather than across the GB plane at $300\text{\,}\mathrm{K}$, i.e.,
$D^{\mathrm{GB}}_{\parallel}>D^{\mathrm{GB}}_{\perp}$. The calculated GB
diffusion coefficients are within the interval set by the anion-ordered and
anion-disordered bulk diffusion coefficients for temperatures below
$1\,000\text{\,}\mathrm{K}$. Table 5 shows that the diffusion activation
energies of all three GBs are also within the activation energy range of the
anion-ordered and anion-disordered bulk structures
($216384\text{\,}\mathrm{meV}$).

Experimental data obtained by nuclear magnetic resonance measurements from
Refs. [7, 61, 8] and obtained by polarization measurements from Ref. [31] are
shown in Fig. 11 and Table 5 for comparison. For reference, ionic
conductivities, $\sigma$, which are converted from the calculated diffusion
coefficients from atomistic simulations (see Appendix B), are also shown in
Table 5. The available experimental data, including all the activation
energies, the diffusion coefficients, and the conductivities, fall between or
together with the computed and Arrhenius-extrapolated values of the
investigated atomic structures.

##  IV Discussion

We have proposed an active learning scheme based on MTPs of progressively
increasing quality. A systematic validation for the thus-obtained MTPs has
also been conducted. Our results demonstrate that the proposed active learning
scheme is able to automatically and efficiently sample configurations for
constructing the training set of complex atomic structures in Li6PS5Cl,
including GBs, starting from ab initio MD simulations of only the anion-
ordered bulk structure. Importantly, diffusion coefficients obtained from the
final high-level MTPs show small variations with respect to the training set
and MTP parameter space, guaranteeing the reproducibility of atomistic
simulation results using machine-learning interatomic potentials. We note that
active learning ideas akin to our current scheme were investigated recently
for structurally and chemically simpler systems, i.e., a unary metallic system
(Zr) with defects [62] and a defect-free ternary system (Ti0.5Al0.5N) [63].
With the verified applicability to a structurally and chemically complex
material system, Li6PS5Cl, the present work systematically brings active
learning to the next development stage.

The high-level MTPs trained with the proposed scheme show energy RMSEs from
$67\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$,
which are higher than that found for metallic systems, either with defects (up
to
$4\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$)
[64, 65, 66] or without defects
($13\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$)
[67, 68, 69, 70, 71]. This comparison indicates the inherent complexity of the
electronic interactions in Li6PS5Cl. Since the computational time grows
quickly for high MTP levels without substantial reduction of the RMSEs, a
further increase of the MTP level does not seem to be a viable option. In the
future, the proposed active learning scheme could be coupled with other recent
machine-learning potentials, e.g., NequIP [72] or MACE [73], to further
optimize the balance of accuracy and efficiency.

The GB energies obtained with the MTPs and the annealing-and-quenching (a$+$q)
approach for three structurally distinct GBs in Li6PS5Cl are within
$719\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{2}$,
which translates to $0.10.3\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$. The
range is consistent with that calculated by DFT for GBs in Li6PS5Cl with
smaller simulation cells
($0.090.54\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$) [74]. A close value of
$0.26\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$ was reported in Ref. [19]
for Li6PS5Br, a similar material system. The calculated GB energies for
Li6PS5Cl are mostly lower than those for other recently explored solid
electrolytes, e.g., Li3PS4
($0.21.2\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$) [42], Li0.16La0.62TiO3
($0.31.3\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$) [75] or Li3OCl
($0.31.1\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$) [40, 76].

In the present study, we have defined the GB width to be six times the
standard deviation (covering 99.7%) of the Gaussian fitted to the relative
atomic energy profiles of Li. In this way, the GB regions cover most of the Li
layers showing GB-induced concentration deviations with respect to the bulk,
resulting in a negligible Li segregation. The GB widths determined in this way
are within $1530\text{\,}\mathrm{\text{Å}}$, and thus larger than GB width
estimates for other solid electrolytes in previous computational studies,
e.g., $\beta$-Li3PS4 ($710\text{\,}\mathrm{\text{Å}}$) [40] or
Li0.375Sr0.4375Ta0.75Zr0.25O3 ($510\text{\,}\mathrm{\text{Å}}$) [41]. Our
results indicate that GBs have a relatively large impact on the local
structural modification of complex solid electrolytes. The computed GB
diffusion coefficients show sensitivity to the GB width employed in the
calculations, as demonstrated in Fig. 12. For Li6PS5Cl, where GBs enhance Li
diffusion compared to the anion-ordered bulk, smaller estimates of GB widths
result in larger GB diffusion coefficients and vice versa. With a GB width
larger than $7\text{\,}\mathrm{\text{Å}}$, the diffusion coefficients at
$300\text{\,}\mathrm{K}$ computed for all three GBs are between that computed
for the anion-ordered and the anion-disordered bulk structures.

Figure 12: Calculated GB diffusion coefficients at $300\text{\,}\mathrm{K}$ as
a function of the utilized GB width. Diffusion coefficients of the anion-
ordered and the anion-disordered bulk structures at $300\text{\,}\mathrm{K}$
are shown for comparison. In the present study, the GB widths of the three
investigated GBs were determined to be between
$1530\text{\,}\mathrm{\text{Å}}$.

Our simulations show that all three GBs enhance the Li self-diffusion in
Li6PS5Cl compared to the anion-ordered bulk, which is consistent with a
previous study for Li6PS5Br [19]. The effect of GBs on Li diffusion depends on
the characteristic atomic arrangements induced by the GBs. The $\Sigma
3(1\bar{1}2)[110]$ GB breaks the Li-coordinated cages and thus enhances Li
diffusion, which we refer to as the “cage-opening effect”. A similar
enhancement of Li diffusion was found due to the effect of anion disorder [9,
17]. For $\Sigma 3(\bar{1}11)[110]$, fewer structural changes (distortions)
are seen in the GB regions, and correspondingly, less long-range Li diffusion
is observed in the MD simulations compared to $\Sigma 3(1\bar{1}2)[110]$. The
amorphous-like structure formed in the GB $\Sigma 5(001)[001]$ increases the
Li diffusion. Our findings are in line with the Li-diffusion enhancement due
to amorphization reported for another sulfide-type solid electrolyte, i.e.,
Li3PS4 [77, 42].

At temperatures below $1\,000\text{\,}\mathrm{K}$, the calculated GB diffusion
coefficients fall within the interval set by the Arrhenius fits for the
diffusion coefficients of the anion-ordered and anion-disordered bulk
structures. Similar phenomena were also seen for other solid electrolytes: For
Li6PS5Br the MSDs of Li atoms in the GB regions are in between those for the
anion-ordered and the anion-disordered structures [19], and for $\beta$-Li3PS4
GB diffusion coefficients are in between those for the bulk and the amorphous
structures [42]. We expect that Li diffusion coefficients of other (high-
angle) GBs in Li6PS5Cl will also fall into the here-established Arrhenius
interval. For the anion-disordered structures of argyrodites Li6PS${}_{5}X$
($X$ $\in\\{\text{Cl, Br, I}\\}$), both experiments [5, 78] and simulations
[17, 19] indicate that maximum diffusivity is likely achieved at about
$50\text{\,}\mathrm{\char 37\relax}$ $X$/S-anion disorder. Therefore, we
expect that the Li diffusion coefficient of the bulk structure with a
different anion-disorder ratio may also fall into the obtained order-disorder
interval.

The Borisov relation [79, 80] and its derivatives [81] are often applied to
establish correlations between GB energy and diffusivity for metallic
materials [82, 83, 81, 84]. For Li6PS5Cl, the $0\text{\,}\mathrm{K}$ GB energy
of $\Sigma 3(1\bar{1}2)[110]$ is larger than that of $\Sigma
3(\bar{1}11)[110]$. The diffusion coefficient of $\Sigma 3(1\bar{1}2)[110]$ is
also larger than that of $\Sigma 3(\bar{1}11)[110]$ over a wide temperature
range. These results indicate that the Borisov relation may be applied for
Li6PS5Cl, but further research on more GBs is required to validate the
hypothesis.

Since the available experimental diffusivity data for Li6PS5Cl correspond to
polycrystalline samples [7, 8, 61, 31], it is meaningful to analyze the
diffusion coefficients of polycrystalline Li6PS5Cl. Figure 13 shows an
extrapolation of the diffusion coefficients obtained from the atomistic
simulations at $300\text{\,}\mathrm{K}$ to the macro-scale based on the Wiener
bounds [85, 86] (see Appendix C). The color-shaded regions in Fig. 13(a)
indicate the spread of the effective macroscopic diffusion coefficient
$D_{\textrm{macro}}$ for a polycrystalline microstructure with GBs, assuming
either the anion-ordered or anion-disordered bulk structure. Both color-shaded
regions narrow down with increasing grain size, and the impact of GBs becomes
nearly negligible above a grain size of $1\,000\text{\,}\mathrm{nm}$. This
finding is consistent with experiments that reported a minor influence of GBs
on Li transport [8, 61] for grain sizes of about $1\,500\text{\,}\mathrm{nm}$
[8]. To analyze a joint contribution of the anion-ordered and the anion-
disordered bulk regions in a macroscopic sample, Fig. 13(b) shows the Wiener
bounds for a mixture of the anion-ordered and anion-disordered regions
(without an additional GB impact). Compared with experimental data, a
significant interval of possible $D_{\textrm{macro}}$ values is observed over
a large volume fraction range of the two regions (about
$1095\text{\,}\mathrm{\char 37\relax}$). To systematically examine the
experimental and simulation results and to investigate the combined effect of
different bulk structures and GBs, continuum simulations with explicit
consideration of the microstructure are needed [87], which is planned for a
subsequent study. Additionally, the impact of other defects, e.g., dislocation
[88] and micro-pores [89], and strain [90] may also be considered.

Figure 13: Extension of diffusion coefficients at $300\text{\,}\mathrm{K}$
from atomistic (Fig. 11 and Table 5) to macro-scale based on the Wiener bounds
[85, 86] [Eq. (10) for the upper bound and Eq. (11) for the lower bound]. (a)
Grain size ($d$) effect for a polycrystal with anion-ordered or
$50\text{\,}\mathrm{\char 37\relax}$ anion-disordered bulk structures. The
shaded areas are the intervals set by the Wiener bounds. (b) Effect of mixing
the anion-ordered and the $50\text{\,}\mathrm{\char 37\relax}$ anion-
disordered bulk structures for a single crystal without GBs. The range
covering diffusion coefficients measured in experiments at about
$300\text{\,}\mathrm{K}$ [7, 8] is shown for comparison. It was evidenced from
experiments that the sample had an averaged grain size of less than
$1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ [8].

##  V Conclusions

Diffusion in solid electrolytes is a complex, multiscale phenomenon influenced
by microstructure and chemistry, particularly grain boundaries (GBs) and anion
disorder. This multiscale challenge pushes simulations to their limits as
classical force fields lack accuracy and large-scale ab initio methods are
computationally out of reach.

We have tackled this challenge and developed a quality-level-based active
learning scheme to efficiently and systematically train accurate machine-
learning interatomic potentials for complex atomic structures. These
potentials enable, for example, the acquisition of diffusivity data for GBs
and anion-disordered structures in solid electrolytes through accelerated
atomistic simulations with near ab initio accuracy.

Utilizing the proposed scheme, we have investigated Li-ion diffusion for three
structurally distinct GBs in Li6PS5Cl. These GBs exhibit low formation
energies, indicating their high stability in polycrystalline Li6PS5Cl. The GBs
enhance Li-ion diffusion compared to the anion-ordered bulk structure, with
the degree of enhancement varying according to the specific GB structure. The
underlying reason for the enhancement is traced back to the Li-cage opening
effect in the GB region. GBs may, thus, generally affect the macroscopic
diffusivity of polycrystalline Li6PS5Cl. Based on the present data, the
limiting case for diffusion enhancement is the 50% anion-disordered bulk.
Experimental diffusion data for Li6PS5Cl around room temperature fall into the
wide Arrhenius-extrapolated interval of diffusion coefficients for the
investigated atomic structures.

The proposed scheme facilitates the computation of energies and diffusivities
for various complex atomic structures in solid electrolytes. These data can,
for example, be subsequently integrated into continuum simulations to model
large-scale microstructures comprising multiple complex atomic structures.
Consequently, high-precision simulations of macroscale diffusion in solid
electrolytes are within reach, with potential applications in microstructure
engineering.

###### Acknowledgements.

This project is funded by Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC 2075 –
390740016. We acknowledge the support by the Stuttgart Center for Simulation
Science (SimTech), the state of Baden-Württemberg through bwHPC, and the
German Research Foundation (DFG) through grant no. INST 40/575-1 FUGG (JUSTUS
2 cluster) – 405998092. Yuji Ikeda is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – IK 125/1-1 –
519607530. Felix Fritzen is funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) in the Heisenberg program – FR 2702/10 –
517847245. Sergiy Divinski thanks the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) for partial support via the grant DI 1419/18-1.
Blazej Grabowski is funded by the European Research Council (ERC) under the
European Union’s Horizon 2020 Research and Innovation Program (grant agreement
no. 865855).

##  Appendix A Details on training and validation of MTP

A detailed flowchart of the proposed quality-level-based active learning
scheme is shown in Fig. 14. This scheme consists of three core parts:
Initialization, pre-training, and standard active learning.

Figure 14: Detailed flowchart of the quality-level-based active learning
scheme proposed and utilized in the present study. “Sub-sample” and “select-
add” are commands in the mlip package [54]. The MTP level [54] is denoted by
$m$, and $i$ enumerates the different, retrained MTPs of the same level. The
“pre-training” and “standard active learning” stages are performed at each MTP
level, as shown in Fig. 2. The number of “pre-training” cycles is controlled
by the hyperparameter $i_{\textrm{pre}}$. When the desired accuracy of the MTP
is reached, the proposed scheme finishes. In the present study, we used the
anion-ordered bulk structure as the input structure in the “initialization”
and set $i_{\textrm{pre}}=1$. MTP levels up to 20 ($M=20$) were tested
separately for different target structures in Li6PS5Cl (anion-ordered and
anion-disordered bulk, and GBs). The large-scale production runs were done
with level 18 MTPs to balance the accuracy against the computational costs
(Fig. 4).

In the initialization, ab initio MD simulations are performed for a small
supercell of the ideal (ordered) bulk structure of the target material to
provide the basis for the training set. To reduce the size of the training set
(which enables faster MTP fitting at the initial stage), a small number of
(several hundred) uncorrelated configurations are extracted (the so-called
“sub-sample”) from the obtained ab initio MD trajectories. These
configurations, together with the DFT-computed energies, forces, and stresses,
are then included in the training set. Note that, in the proposed active
learning scheme, ab initio MD simulations are performed only during the
initialization stage and that they are comparably efficient because of the
utilized ideal bulk structure and the small supercell size.

The labeled configurations obtained in the initialization are propagated to
the pre-training stage. The key idea behind the pre-training stage is to pre-
train the MTP on the ideal bulk structure before entering the standard active
learning stage for the target structures. The ideal bulk structure already
contains a lot of information about the atomic interactions, and thus, the
following standard active learning cycle for the target structure is shortened
(i.e., fewer DFT calculations and fewer standard active learning cycles are
needed). The pre-training stage can be run once or more times (depending on
$i_{\mathrm{pre}}$; we chose one) to select the most informative
configurations (so-called “select-add” based on the D-optimality criterion
[91, 92]) from the initial ab initio molecular dynamics runs. The selected
configurations are added to the training set, and the MTP is retrained. When
the pre-training stage is entered for the first time after initialization, an
MTP of level 4 ($m=4$) is chosen for the fitting. In later stages, higher-
level MTPs ($m=6,8,10,$ etc.) are chosen.

The pre-trained MTP and the accumulated training set (from initialization and
pre-training) are propagated to the next stage, which is a standard active
learning cycle carried out for the target structure (e.g., a structure with a
grain boundary). The pre-trained MTP is used to perform MD simulations for the
target structure, and configurations with large extrapolation grades are
selected (so-called “sample”). The energies, forces, and stresses of the
sampled configurations are calculated by DFT and added to the training set.
The MTP is then retrained using the training set enlarged with the sampled
configurations. The standard active learning for one target structure is
finished when no more configurations are sampled. In case the MTP is required
to describe multiple target structures, e.g., different anion-disordered
structures of Li6PS5Cl or different types of GBs, the standard active learning
is repeated for each target structure with the MTP of the same level ($i$
increases while $m$ is unchanged).

After running the standard active learning cycle for the first time, the
trained level 4 MTP ($m=4$) is outputted, and the accuracy of the MTP is
evaluated according to the fitting RMSEs in energies and forces. If higher
accuracy is required, an untrained MTP of the next level (level 6, $m=6$) is
utilized for another training round ($i$ reset to 1) with the pre-training and
standard active learning stages. The new configurations generated during each
stage are added to and accumulated in the training set. The MTP level in the
pre-training and standard active learning can continue to increase in this way
until the accuracy requirement of the MTP is satisfied. After that, the
highest level of the utilized MTPs is denoted by $M$, and the whole active
learning scheme finishes by outputting the finally-trained level $M$-MTP. With
this quality-level-based scheme, the configuration space of the target
structures is explored efficiently and systematically. It is important to note
that the computational cost of high-level MTPs increases rapidly,
necessitating a balance between accuracy and efficiency (see Sec. III.1). An
automatic code to perform the proposed active learning scheme is available
from the corresponding author upon reasonable request.

In the present study, conventional unit cells of the anion-ordered bulk
structure of Li6PS5Cl [Fig. 1(a)] in three volumes with atomic densities
17.73; 20.68; 23.943/atom were used in the initialization stage. Ab initio MD
simulations were performed at $1\,500\text{\,}\mathrm{K}$ for
$8\text{\,}\mathrm{ps}$ with the Nosé–Hoover thermostat and the canonical
($NVT$) ensemble implemented in vasp. A time step of $2\text{\,}\mathrm{fs}$
and a Nosé mass of
$3\text{\,}\mathrm{u}\,\mathrm{\text{Å}}\textsuperscript{2}$ were utilized. At
each ionic step, the energy was converged to within
${10}^{-3}\text{\,}\mathrm{eV}$ per simulation cell. The MD trajectories were
“sub-sampled” every $1\text{\,}\mathrm{ps}$. A set of untrained MTPs
configured at different levels [54] and with 8 radial basis functions was
utilized. The minimum and maximum cutoffs were set to 1.5;5, respectively. The
fitting weights for the energies, the atomic forces, and the stresses were set
to 1, $0.1\text{\,}\mathrm{\text{Å}}\textsuperscript{2}$, and
$0.001\text{\,}\mathrm{\text{Å}}\textsuperscript{6}$, respectively. In the
standard active learning stage, “sample” was performed by MTP MD simulations
with lammps at a temperature of $1\,000\text{\,}\mathrm{K}$ and for a maximum
of $0.4\text{\,}\mathrm{ns}$. The NPT ensemble with the same parameters as
described in Sec. II.4 was used. The extrapolation grade values of the
selecting and breaking thresholds were chosen to be 2 and 5, as suggested in
Ref. [55].

Target structures for the anion-ordered and anion-disordered bulk, as well as
the three GBs of Li6PS5Cl were constructed. A $1\times 1\times 2$ supercell
was used for the ordered bulk structure. The disordered bulk structure was
constructed based on a $1\times 2\times 2$ supercell of the ordered bulk
structure (see Sec. II.2). Three different anion-disordered structures with
random anion disorder were used to train MTPs of the anion-disordered bulk at
each level. In total, 2210 configurations were accumulated in the training set
up to level 18 (Table 3). GB target structures were constructed following Sec.
II.1 and with supercells made of the coincidence-site lattice unit cell. The
validation sets were constructed separately for the anion-ordered and the
anion-disordered bulk structures. For each target structure, the final MTP
trained up to level 20 was used to perform an MD simulation at
$600\text{\,}\mathrm{K}$. The validation set of each structure contained 500
configurations “sub-sampled” every $1\text{\,}\mathrm{ps}$ from the MD
simulation. The energies, forces, and stresses of those configurations were
calculated with DFT. Note that for the anion-disordered bulk structure, the
anion disorder of the configurations in the validation set differed from those
in the training set.

##  Appendix B Conversion of self-diffusion coefficients to ionic
conductivities

Based on the Nernst–Einstein equation [93, 94], a relation between ionic
conductivity, $\sigma$, and the self-diffusion coefficient, $D$, at
temperature, $T$, can be written as,

$\sigma(T)=\Lambda\frac{D(T)e^{2}}{k_{\mathrm{B}}T},$ (9)

where $e$ is the elementary charge. The coefficient $\Lambda$ is related to
the concentration of mobile ions, i.e., the Li-ions in the present case, and
to the Haven ratio [93, 94]. Table 6 lists the experimentally measured values
for $\sigma$ and $D$ at about $300\text{\,}\mathrm{K}$ from different studies,
and the resulting $\Lambda$’s [Eq. (9)]. As reflected by the strong scatter in
the $\Lambda$ values from the different studies, the determination of
$\Lambda$ is not trivial [95, 7]. We utilized
$\Lambda=$21.808\text{\times}{10}^{-3}\text{\,}\mathrm{L}\mathrm{i}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{3}$$
from the most recent experimental study [8] to convert the calculated self-
diffusion coefficients to ionic conductivity (Table 5).

Table 6: Comparison of experimentally measured ionic conductivity $\sigma$ and self-diffusion coefficient $D$ at about $300\text{\,}\mathrm{K}$. The coefficient $\Lambda$ is calculated based on the measured values according to Eq. (9).  Year | $T$ | $\sigma$ | $D$ | $\Lambda$ | Ref.
---|---|---|---|---|---
| ($\text{\,}\mathrm{K}$) | ($\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}$) | (${10}^{-7}\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$) | (${10}^{-3}\text{\,}\mathrm{L}\mathrm{i}\mathrm{/}\mathrm{\text{Å}}\textsuperscript{3}$) |
2011 | 313 | $0.001\,3$ | $0.77$ | $0.003$ | [31]
2019 | 313 | $9$ | $0.67$ | $22.625$ | [61]
2019 | 298 | $2.5$ | $0.387$ | $10.359$ | [7]
2020 | 298 | $3.4$ | $0.25$ | $21.808$ | [8]

##  Appendix C Extension from atomistic to macro-scale

Models with different geometrical assumptions of grains have been used to
obtain effective diffusion coefficients of polycrystals, e.g., the
Maxwell–Garnet model [96] and its modification [97], the Belova–Murch [44]
model, or the Chen–Schuh [98] model. In continuum mechanics, the Voigt–Reuss
bounds [99, 100, 101] are commonly used to give an admissible range of
effective properties rather than a geometry-dependent guess. They are
equivalent to the Wiener bounds [85, 86] used for permittivity. In the present
study, we used the Wiener bounds to obtain diffusion coefficients at the
macro-scale ($D_{\textrm{macro}}$) from the atomistic diffusion coefficients
(Table 5) at an exemplary temperature of $300\text{\,}\mathrm{K}$.

For a polycrystal with a set of components $J$, e.g., the ordered bulk, the
disordered bulk, and the GBs, etc., Chen et al. [86] have shown that, in the
absence of segregation, if diffusivity ratios of any two components are less
than ${10}^{2}\text{\,}$ (low-contrast system), the Wiener bounds remain valid
for evaluating the effective diffusion coefficients. Specifically, the upper
and lower bounds of the diffusion coefficients of the polycrystal are given by
the weighted arithmetic mean,

$D^{\textrm{upper}}_{\textrm{macro}}=\sum\limits_{j\in J}f^{j}D^{j},$ (10)

and the weighted harmonic mean,

$D^{\textrm{lower}}_{\textrm{macro}}=\left(\sum\limits_{j\in
J}\frac{f^{j}}{D^{j}}\right)^{-1},$ (11)

respectively. Here, $f^{j}$ is the volume fraction of the structure $s$ in the
polycrystal which satisfies,

$\sum\limits_{j\in J}f^{j}=1.$ (12)

Equations (10) and (11) are equivalent to the equations of the Hart [43, 44]
and the one-dimensional Maxwell–Garnet models [96], respectively. As was
validated in Ref. [102], effective diffusion coefficients of a polycrystal
calculated from the above-mentioned classical models fall within the Wiener
bounds [Eqs. (10) and (11)].

Assuming for simplicity an isotropic distribution of GB surface orientations,
the averaged GB diffusion coefficients can be estimated as an arithmetic mean
of the diffusion coefficients at the atomistic scale in three dimensions (Fig.
11 and Table 5),

$D^{\textrm{GB}_{k}}\approx\frac{2D_{\parallel}^{\mathrm{GB}_{k}}+D_{\perp}^{\mathrm{GB}_{k}}}{3},$
(13)

where $k$ indicates the GB type. With an average grain size of $d$, the volume
fraction $f^{\mathrm{GB}_{k}}$ of the type-$k$ GB can be approximated [103,
98],

$f^{\textrm{GB}_{k}}\approx
p^{\mathrm{GB}_{k}}H_{\mathrm{GB}}\frac{\delta^{\mathrm{GB}_{k}}}{d},$ (14)

where $\delta^{\textrm{GB}_{k}}$ is the width of the type-$k$ GB (Table 2).
Further, $H_{\textrm{GB}}$ is a dimensionless numerical factor accounting for
the shape and size distribution of the grains in a polycrystal. Here,
$H_{\textrm{GB}}=$2.910\,5\text{\,}$$ was used, corresponding to the grain
shape of Voronoi polyhedra and a nearly log-normal grain-size distribution
[103, 98]. The formation probability of the type-$k$ GB is reflected by
$p^{\mathrm{GB}_{k}}$. Assuming for simplicity all GB types to appear with
equal probability gives, for every $k$, $p^{\mathrm{GB}_{k}}\approx
1/n^{\mathrm{GB}}$, where $n^{\mathrm{GB}}$ is the total number of GB types.
The analysis has to be modified if Li ions are not equipartitionally
distributed between different components, i.e., in the case of segregation
[104].

We note that the analysis above corresponds to the so-called A-type kinetic
regime of GB diffusion [105], which holds if the diffusion length,
$\sqrt{D_{\textrm{macro}}t}$, is larger than the characteristic size of all
microstructure elements in the model [106]. In the absence of segregation, the
ratio of $\sqrt{D_{\textrm{macro}}t}$ to the characteristic size of any
microstructure element is typically required to exceed three [106]. Otherwise,
one has to deal with diffusion in a complex microstructure with a hierarchy of
the kinetic properties from different constituting elements (e.g., GBs of
different types, dislocations, pores, etc., and their specific geometric
arrangement), which also results in a hierarchy of diffusion regimes [107,
108].

## References

  * Famprikis  _et al._ [2019] T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam, and C. Masquelier, Fundamentals of inorganic solid-state electrolytes for batteries, Nat. Mater. 18, 1278 (2019).
  * Zhao  _et al._ [2020] Q. Zhao, S. Stalin, C.-Z. Zhao, and L. A. Archer, Designing solid-state electrolytes for safe, energy-dense batteries, Nat. Rev. Mater. 5, 229 (2020).
  * Deiseroth  _et al._ [2008] H.-J. Deiseroth, S.-T. Kong, H. Eckert, J. Vannahme, C. Reiner, T. Zaiß, and M. Schlosser, Li6PS${}_{5}X$: A class of crystalline Li-rich solids with an unusually high Li+ mobility, Angew. Chem. Int. Ed. 47, 755 (2008).
  * Boulineau  _et al._ [2013] S. Boulineau, J.-M. Tarascon, J.-B. Leriche, and V. Viallet, Electrochemical properties of all-solid-state lithium secondary batteries using Li-argyrodite Li6PS5Cl as solid electrolyte, Solid State Ion. 242, 45 (2013).
  * Kraft  _et al._ [2017] M. A. Kraft, S. P. Culver, M. Calderon, F. Böcher, T. Krauskopf, A. Senyshyn, C. Dietrich, A. Zevalkink, J. Janek, and W. G. Zeier, Influence of lattice polarizability on the ionic conductivity in the lithium superionic argyrodites Li6PS${}_{5}X$ ($X$ = Cl, Br, I), J. Am. Chem. Soc. 139, 10909 (2017).
  * Ruhl  _et al._ [2021] J. Ruhl, L. M. Riegger, M. Ghidiu, and W. G. Zeier, Impact of solvent treatment of the superionic argyrodite Li6PS5Cl on solid-state battery performance, Adv. Energy Sustain. Res. 2, 2000077 (2021).
  * Adeli  _et al._ [2019] P. Adeli, J. D. Bazak, K. H. Park, I. Kochetkov, A. Huq, G. R. Goward, and L. F. Nazar, Boosting solid-state diffusivity and conductivity in lithium superionic argyrodites by halide substitution, Angew. Chem. Int. Ed. 58, 8681 (2019).
  * Schlenker  _et al._ [2020] R. Schlenker, A.-L. Hansen, A. Senyshyn, T. Zinkevich, M. Knapp, T. Hupfer, H. Ehrenberg, and S. Indris, Structure and diffusion pathways in Li6PS5Cl argyrodite from neutron diffraction, pair-distribution function analysis, and NMR, Chem. Mater. 32, 8420 (2020).
  * Minafra  _et al._ [2020] N. Minafra, M. A. Kraft, T. Bernges, C. Li, R. Schlem, B. J. Morgan, and W. G. Zeier, Local charge inhomogeneity and lithium distribution in the superionic argyrodites Li6PS${}_{5}X$ ($X$ = Cl, Br, I), Inorg. Chem. 59, 11009 (2020).
  * Ganapathy  _et al._ [2019] S. Ganapathy, C. Yu, E. R. H. van Eck, and M. Wagemaker, Peeking across grain boundaries in a solid-state ionic conductor, ACS Energy Lett. 4, 1092 (2019).
  * Milan and Pasta [2022] E. Milan and M. Pasta, The role of grain boundaries in solid-state Li-metal batteries, Mater. Futures 2, 013501 (2022).
  * Deng  _et al._ [2017] Z. Deng, Z. Zhu, I.-H. Chu, and S. P. Ong, Data-driven first-principles methods for the study and design of alkali superionic conductors, Chem. Mater. 29, 281 (2017).
  * Stamminger  _et al._ [2019] A. R. Stamminger, B. Ziebarth, M. Mrovec, T. Hammerschmidt, and R. Drautz, Ionic conductivity and its dependence on structural disorder in halogenated argyrodites Li6PS${}_{5}X$ ($X$ = Br, Cl, I), Chem. Mater. 31, 8673 (2019).
  * Jiang  _et al._ [2022] M. Jiang, Z.-W. Chen, A. Rao, L.-X. Chen, X.-T. Zu, and C. V. Singh, Se-doped Li6PS5Cl and Li5.5PS4.5Cl1.5 with improved ionic conductivity and interfacial compatibility: A high-throughput DFT study, J. Mater. Chem. C 10, 18294 (2022).
  * Jeon  _et al._ [2024] T. Jeon, G. H. Cha, and S. C. Jung, Understanding the anion disorder governing lithium distribution and diffusion in an argyrodite Li6PS5Cl solid electrolyte, J. Mater. Chem. A 12, 993 (2024).
  * de Klerk  _et al._ [2016] N. J. J. de Klerk, I. Rosłoń, and M. Wagemaker, Diffusion mechanism of Li argyrodite solid electrolytes for Li-ion batteries and prediction of optimized halogen doping: The effect of Li vacancies, halogens, and halogen disorder, Chem. Mater. 28, 7955 (2016).
  * Morgan [2021] B. J. Morgan, Mechanistic origin of superionic lithium diffusion in anion-disordered Li6PS${}_{5}X$ argyrodites, Chem. Mater. 33, 2004 (2021).
  * Lee  _et al._ [2022] B. D. Lee, J.-W. Lee, J. Park, M. Y. Cho, W. B. Park, and K.-S. Sohn, Argyrodite configuration determination for DFT and AIMD calculations using an integrated optimization strategy, RSC Adv. 12, 31156 (2022).
  * Sadowski [2023] M. Sadowski,  _Properties of sulfide solid electrolytes studied by electronic structure calculations_ , Ph.D. thesis, Technische Universität Darmstadt, Darmstadt (2023).
  * He  _et al._ [2018] X. He, Y. Zhu, A. Epstein, and Y. Mo, Statistical variances of diffusional properties from ab initio molecular dynamics simulations, Npj Comput. Mater. 4, 1 (2018).
  * Das  _et al._ [2022] T. Das, B. V. Merinov, M. Y. Yang, and W. A. G. Iii, Structural, dynamic, and diffusion properties of a Li6(PS4)SCl superionic conductor from molecular dynamics simulations; prediction of a dramatically improved conductor, J. Mater. Chem. A 10, 16319 (2022).
  * Stegmaier  _et al._ [2021] S. Stegmaier, R. Schierholz, I. Povstugar, J. Barthel, S. P. Rittmeyer, S. Yu, S. Wengert, S. Rostami, H. Kungl, K. Reuter, R.-A. Eichel, and C. Scheurer, Nano-scale complexions facilitate Li dendrite-free operation in LATP solid-state electrolyte, Adv. Energy Mater. 11, 2100707 (2021).
  * Stegmaier  _et al._ [2022] S. Stegmaier, K. Reuter, and C. Scheurer, Exploiting nanoscale complexion in LATP solid-state electrolyte via interfacial Mg2+ doping, Nanomater. 12, 2912 (2022).
  * Gubaev  _et al._ [2023] K. Gubaev, V. Zaverkin, P. Srinivasan, A. I. Duff, J. Kästner, and B. Grabowski, Performance of two complementary machine-learned potentials in modelling chemically complex systems, Npj Comput. Mater. 9, 1 (2023).
  * Erhard  _et al._ [2024] L. C. Erhard, J. Rohrer, K. Albe, and V. L. Deringer, Modelling atomic and nanoscale structure in the silicon–oxygen system through active machine learning, Nat. Commun. 15, 1927 (2024).
  * Novoselov  _et al._ [2019] I. I. Novoselov, A. V. Yanilkin, A. V. Shapeev, and E. V. Podryabinkin, Moment tensor potentials as a promising tool to study diffusion processes, Comput. Mater. Sci. 164, 46 (2019).
  * Winter and Gómez-Bombarelli [2023] G. Winter and R. Gómez-Bombarelli, Simulations with machine learning potentials identify the ion conduction mechanism mediating non-Arrhenius behavior in LGPS, J. Phys. Energy 5, 024004 (2023).
  * Musaelian  _et al._ [2023] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nat. Commun. 14, 579 (2023).
  * Shapeev [2016] A. V. Shapeev, Moment tensor potentials: A class of systematically improvable interatomic potentials, Multiscale Model. Simul. 14, 1153 (2016).
  * Kong  _et al._ [2010] S.-T. Kong, H.-J. Deiseroth, C. Reiner, Ö. Gün, E. Neumann, C. Ritter, and D. Zahn, Lithium argyrodites with phosphorus and arsenic: Order and disorder of lithium atoms, crystal chemistry, and phase transitions, Chem. Eur. J. 16, 2198 (2010).
  * Deiseroth  _et al._ [2011] H.-J. Deiseroth, J. Maier, K. Weichert, V. Nickel, S.-T. Kong, and C. Reiner, Li7PS6 and Li6PS${}_{5}X$ ($X$: Cl, Br, I): Possible three-dimensional diffusion pathways for lithium ions and temperature dependence of the ionic conductivity by impedance measurements, Z. Anorg. Allg. Chem. 637, 1287 (2011).
  * Wang and Shao [2017] Z. Wang and G. Shao, Theoretical design of solid electrolytes with superb ionic conductivity: Alloying effect on Li+ transportation in cubic Li6PA${}_{5}X$ chalcogenides, J. Mater. Chem. A 5, 21846 (2017).
  * Golov and Carrasco [2021] A. Golov and J. Carrasco, Molecular-level insight into the interfacial reactivity and ionic conductivity of a Li-argyrodite Li6PS5Cl solid electrolyte at bare and coated Li-metal anodes, ACS Appl. Mater. Interfaces 13, 43734 (2021).
  * D’Amore  _et al._ [2022] M. D’Amore, L. E. Daga, R. Rocca, M. F. Sgroi, N. L. Marana, S. M. Casassa, L. Maschio, and A. M. Ferrari, From symmetry breaking in the bulk to phase transitions at the surface: A quantum-mechanical exploration of Li6PS5Cl argyrodite superionic conductor, Phys. Chem. Chem. Phys. 24, 22978 (2022).
  * Ou  _et al._ [2022] Y. Ou, Y. Ikeda, O. Clemens, and B. Grabowski, Dynamic stabilization of perovskites at elevated temperatures: A comparison between cubic BaFeO3 and vacancy-ordered monoclinic BaFeO2.67, Phys. Rev. B 106, 064308 (2022).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. 44, 1272 (2011).
  * Cheng  _et al._ [2018] J. Cheng, J. Luo, and K. Yang, Aimsgb: An algorithm and open-source Python library to generate periodic grain boundary structures, Comput. Mater. Sci. 155, 92 (2018).
  * Metzler and Klafter [2000] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep. 339, 1 (2000).
  * Einstein [1905] A. Einstein, Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen, Ann. Phys. 322, 549 (1905).
  * Dawson  _et al._ [2018] J. A. Dawson, P. Canepa, T. Famprikis, C. Masquelier, and M. S. Islam, Atomic-scale influence of grain boundaries on Li-ion conduction in solid electrolytes for all-solid-state batteries, J. Am. Chem. Soc. 140, 362 (2018).
  * Lee  _et al._ [2023] T. Lee, J. Qi, C. A. Gadre, H. Huyan, S.-T. Ko, Y. Zuo, C. Du, J. Li, T. Aoki, R. Wu, J. Luo, S. P. Ong, and X. Pan, Atomic-scale origin of the low grain-boundary resistance in perovskite solid electrolyte Li0.375Sr0.4375Ta0.75Zr0.25O3, Nat. Commun. 14, 1940 (2023).
  * Jalem  _et al._ [2023] R. Jalem, M. L. H. Chandrappa, J. Qi, Y. Tateyama, and S. P. Ong, Lithium dynamics at grain boundaries of $\beta$-Li3PS4 solid electrolyte, Energy Adv. 2, 2029 (2023).
  * Hart [1957] E. W. Hart, On the role of dislocations in bulk diffusion, Acta Metall. 5, 597 (1957).
  * Belova and Murch [2003] I. V. Belova and G. E. Murch, Diffusion in nanocrystalline materials, J. Phys. Chem. Solids 64, 873 (2003).
  * Kaur  _et al._ [1995] I. Kaur, Y. Mishin, and W. Gust,  _Fundamentals of grain and interface boundary diffusion_ (Wiley & Sons LTD, 1995).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse [1995] G. Kresse, Ab initio molecular dynamics for liquid metals, J. Non-Cryst. Solids 192-193, 222 (1995).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Blöchl  _et al._ [1994] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B 49, 16223 (1994).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS – a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271, 108171 (2022).
  * Wagih and Schuh [2023] M. Wagih and C. A. Schuh, Viewpoint: Can symmetric tilt grain boundaries represent polycrystals?, Scr. Mater. 237, 115716 (2023).
  * Novikov  _et al._ [2020] I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, The MLIP package: Moment tensor potentials with MPI and active learning, Mach. Learn.: Sci. Technol. 2, 025002 (2020).
  * Gubaev  _et al._ [2021] K. Gubaev, Y. Ikeda, F. Tasnádi, J. Neugebauer, A. V. Shapeev, B. Grabowski, and F. Körmann, Finite-temperature interplay of structural stability, chemical complexity, and elastic properties of bcc multicomponent alloys from ab initio trained machine-learning potentials, Phys. Rev. Mater. 5, 073801 (2021).
  * Stukowski [2009] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO – the open visualization tool, Model. Simul. Mater. Sci. Eng. 18, 015012 (2009).
  * Larsen  _et al._ [2017] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, The atomic simulation environment – a Python library for working with atoms, J. Phys. Condens. Matter. 29, 273002 (2017).
  * Glass  _et al._ [2006] C. W. Glass, A. R. Oganov, and N. Hansen, USPEX – Evolutionary crystal structure prediction, Comput. Phys. Commun. 175, 713 (2006).
  * Guo  _et al._ [2017] Y. Guo, Q. Wang, and W. A. Saidi, Structural stabilities and electronic properties of high-angle grain boundaries in perovskite cesium lead halides, J. Phys. Chem. C 121, 1715 (2017).
  * Koju and Mishin [2020] R. K. Koju and Y. Mishin, Relationship between grain boundary segregation and grain boundary diffusion in Cu-Ag alloys, Phys. Rev. Mater. 4, 073403 (2020).
  * Hanghofer  _et al._ [2019] I. Hanghofer, M. Brinek, S. L. Eisbacher, B. Bitschnau, M. Volck, V. Hennige, I. Hanzu, D. Rettenwander, and H. M. R. Wilkening, Substitutional disorder: Structure and ion dynamics of the argyrodites Li6PS5Cl, Li6PS5Br and Li6PS5I, Phys. Chem. Chem. Phys. 21, 8489 (2019).
  * Luo  _et al._ [2023] Y. Luo, J. A. Meziere, G. D. Samolyuk, G. L. W. Hart, M. R. Daymond, and L. K. Béland, A set of moment tensor potentials for zirconium with increasing complexity, J. Chem. Theory Comput. 19, 6848 (2023).
  * Bock  _et al._ [2024] F. Bock, F. Tasnádi, and I. A. Abrikosov, Active learning with moment tensor potentials to predict material properties: Ti0.5Al0.5N at elevated temperature, J. Vac. Sci. Technol. A 42, 013412 (2024).
  * Xu  _et al._ [2023] X. Xu, X. Zhang, A. Ruban, S. Schmauder, and B. Grabowski, Strong impact of spin fluctuations on the antiphase boundaries of weak itinerant ferromagnetic Ni3Al, Acta Mater. 255, 118986 (2023).
  * Xu  _et al._ [2024] X. Xu, X. Zhang, A. Ruban, S. Schmauder, and B. Grabowski, Accurate complex-stacking-fault Gibbs energy in Ni3Al at high temperatures, Scr. Mater. 242, 115934 (2024).
  * Zotov  _et al._ [2024] N. Zotov, K. Gubaev, J. Wörner, and B. Grabowski, Moment tensor potential for static and dynamic investigations of screw dislocations in bcc Nb, Model. Simul. Mat. Sci. Eng. 32, 035032 (2024).
  * Grabowski  _et al._ [2019] B. Grabowski, Y. Ikeda, P. Srinivasan, F. Körmann, C. Freysoldt, A. I. Duff, A. Shapeev, and J. Neugebauer, Ab initio vibrational free energies including anharmonicity for multicomponent alloys, Npj Comput. Mater. 5, 1 (2019).
  * Jung  _et al._ [2023a] J. H. Jung, P. Srinivasan, A. Forslund, and B. Grabowski, High-accuracy thermodynamic properties to the melting point from ab initio calculations aided by machine-learning potentials, Npj Comput. Mater. 9, 1 (2023a).
  * Jung  _et al._ [2023b] J. H. Jung, A. Forslund, P. Srinivasan, and B. Grabowski, Dynamically stabilized phases with full ab initio accuracy: Thermodynamics of Ti, Zr, Hf with a focus on the hcp-bcc transition, Phys. Rev. B 108, 184107 (2023b).
  * Forslund  _et al._ [2023] A. Forslund, J. H. Jung, P. Srinivasan, and B. Grabowski, Thermodynamic properties on the homologous temperature scale from direct upsampling: Understanding electron-vibration coupling and thermal vacancies in bcc refractory metals, Phys. Rev. B 107, 174309 (2023).
  * Zhou  _et al._ [2022] Y. Zhou, P. Srinivasan, F. Körmann, B. Grabowski, R. Smith, P. Goddard, and A. I. Duff, Thermodynamics up to the melting point in a TaVCrW high entropy alloy: Systematic ab initio study aided by machine learning potentials, Phys. Rev. B 105, 214302 (2022).
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun. 13, 2453 (2022).
  * Batatia  _et al._ [2022] I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, and G. Csanyi, MACE: Higher order equivariant message passing neural networks for fast and accurate force fields, in  _Adv. Neural Inf. Process Syst._, Vol. 35, edited by S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Curran Associates, Inc., 2022) pp. 11423–11436.
  * Xie  _et al._ [2024] W. Xie, Z. Deng, Z. Liu, T. Famprikis, K. T. Butler, and P. Canepa, Effects of grain boundaries and surfaces on electronic and mechanical properties of solid electrolytes, Adv. Energy Mater. 14, 2304230 (2024).
  * Symington  _et al._ [2021] A. R. Symington, M. Molinari, J. A. Dawson, J. M. Statham, J. Purton, P. Canepa, and S. C. Parker, Elucidating the nature of grain boundary resistance in lithium lanthanum titanate, J. Mater. Chem. A 9, 6487 (2021).
  * Quirk and Dawson [2023] J. A. Quirk and J. A. Dawson, Design principles for grain boundaries in solid-state lithium-ion conductors, Adv. Energy Mater. 13, 2301114 (2023).
  * Hayamizu  _et al._ [2016] K. Hayamizu, Y. Aihara, T. Watanabe, T. Yamada, S. Ito, and N. Machida, NMR studies on lithium ion migration in sulfide-based conductors, amorphous and crystalline Li3PS4, Solid State Ion. 285, 51 (2016).
  * Zhou  _et al._ [2021] L. Zhou, N. Minafra, W. G. Zeier, and L. F. Nazar, Innovative approaches to Li-argyrodite solid electrolytes for all-solid-state lithium batteries, Acc. Chem. Res. 54, 2717 (2021).
  * Borisov  _et al._ [1964] V. Borisov, V. Golikov, and G. Scherbedinskiy, Relation between diffusion coefficients and grain boundary energy, Phys. Met. Metall. 17, 881 (1964).
  * Pelleg [1966] J. Pelleg, On the relation between diffusion coefficients and grain boundary energy, Philos. Mag.: J. Theor. Exp. Appl. Phys. 14, 595 (1966).
  * Page  _et al._ [2021] D. E. Page, K. F. Varela, O. K. Johnson, D. T. Fullwood, and E. R. Homer, Measuring simulated hydrogen diffusion in symmetric tilt nickel grain boundaries and examining the relevance of the Borisov relationship for individual boundary diffusion, Acta Mater. 212, 116882 (2021).
  * Divinski  _et al._ [2010] S. V. Divinski, G. Reglitz, and G. Wilde, Grain boundary self-diffusion in polycrystalline nickel of different purity levels, Acta Mater. 58, 386 (2010).
  * Prokoshkina  _et al._ [2013] D. Prokoshkina, V. A. Esin, G. Wilde, and S. V. Divinski, Grain boundary width, energy and self-diffusion in nickel: Effect of material purity, Acta Mater. 61, 5188 (2013).
  * Li  _et al._ [2023] C. Li, S. Lu, S. Divinski, and L. Vitos, Theoretical and experimental grain boundary energies in body-centered cubic metals, Acta Mater. 255, 119074 (2023).
  * Karkkainen  _et al._ [2000] K. Karkkainen, A. Sihvola, and K. Nikoskinen, Effective permittivity of mixtures: Numerical validation by the FDTD method, IEEE Trans. Geosci. Remote Sens. 38, 1303 (2000).
  * Chen and Schuh [2006] Y. Chen and C. A. Schuh, Diffusion on grain boundary networks: Percolation theory and effective medium approximations, Acta Mater. 54, 4709 (2006).
  * Heo  _et al._ [2021] T. W. Heo, A. Grieder, B. Wang, M. Wood, T. Hsu, S. A. Akhade, L. F. Wan, L.-Q. Chen, N. Adelstein, and B. C. Wood, Microstructural impacts on ionic conductivity of oxide solid electrolytes from a combined atomistic-mesoscale approach, Npj Comput. Mater. 7, 1 (2021).
  * Faka  _et al._ [2024] V. Faka, M. T. Agne, M. A. Lange, D. Daisenberger, B. Wankmiller, S. Schwarzmüller, H. Huppertz, O. Maus, B. Helm, T. Böger, J. Hartel, J. M. Gerdes, J. J. Molaison, G. Kieslich, M. R. Hansen, and W. G. Zeier, Pressure-induced dislocations and their influence on ionic transport in Li+-conducting argyrodites, J. Am. Chem. Soc. 146, 1710 (2024).
  * Janek and Zeier [2023] J. Janek and W. G. Zeier, Challenges in speeding up solid-state battery development, Nat. Energy 8, 230 (2023).
  * Žguns and Yildiz [2022] P. Žguns and B. Yildiz, Strain sensitivity of Li-ion conductivity in $\beta$-Li3PS4 solid electrolyte, PRX Energy 1, 023003 (2022).
  * Podryabinkin and Shapeev [2017] E. V. Podryabinkin and A. V. Shapeev, Active learning of linearly parametrized interatomic potentials, Comput. Mater. Sci. 140, 171 (2017).
  * Gubaev  _et al._ [2018] K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, Machine learning of molecular properties: Locality and active learning, J. Chem. Phys. 148, 241727 (2018).
  * Uitz  _et al._ [2017] M. Uitz, V. Epp, P. Bottke, and M. Wilkening, Ion dynamics in solid electrolytes for lithium batteries, J. Electroceram. 38, 142 (2017).
  * Kizilyalli  _et al._ [1999] M. Kizilyalli, J. Corish, and R. Metselaar, Definitions of terms for diffusion in the solid state, Pure Appl. Chem. 71, 1307 (1999).
  * Marcolongo and Marzari [2017] A. Marcolongo and N. Marzari, Ionic correlations and failure of Nernst-Einstein relation in solid-state electrolytes, Phys. Rev. Mater. 1, 025402 (2017).
  * Garnett and Larmor [1997] J. C. M. Garnett and J. Larmor, XII. Colours in metal glasses and in metallic films, Philos. Trans. R. Soc. A 203, 385 (1997).
  * Kalnin  _et al._ [2002] J. R. Kalnin, E. A. Kotomin, and J. Maier, Calculations of the effective diffusion coefficient for inhomogeneous media, J. Phys. Chem. Solids 63, 449 (2002).
  * Chen and Schuh [2007] Y. Chen and C. A. Schuh, Geometric considerations for diffusion in polycrystalline solids, J. Appl. Phys. 101, 063524 (2007).
  * Voigt [1889] W. Voigt, Über die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper, Ann. Phys. 274, 573 (1889).
  * Reuss [1929] A. Reuss, Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle, J. Appl. Math. Mech. 9, 49 (1929).
  * Hill [1963] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids 11, 357 (1963).
  * Yoon  _et al._ [2023] J. S. Yoon, H. Sulaimon, and D. J. Siegel, Exploiting grain boundary diffusion to minimize dendrite formation in lithium metal-solid state batteries, J. Mater. Chem. A 11, 23288 (2023).
  * Thorvaldsen [1997] A. Thorvaldsen, The intercept method – 1. Evaluation of grain shape, Acta Mater. 45, 587 (1997).
  * Zhang and Liu [2012] Y. Zhang and L. Liu, On diffusion in heterogeneous media, Am. J. Sci. 312, 1028 (2012).
  * Harrison [1961] L. G. Harrison, Influence of dislocations on diffusion kinetics in solids with particular reference to the alkali halides, Trans. Faraday Soc. 57, 1191 (1961).
  * Paul  _et al._ [2014] A. Paul, T. Laurila, V. Vuorinen, and S. V. Divinski,  _Thermodynamics, diffusion and the Kirkendall effect in solids_ (Springer International Publishing, Switzerland, 2014).
  * Divinski  _et al._ [2004a] S. V. Divinski, J. S. Lee, and C. Herzig, Grain boundary diffusion and segregation in compacted and sintered nanocrystalline alloys, J. Metastable Nanocryst. Mater. 19, 55 (2004a).
  * Divinski  _et al._ [2004b] S. Divinski, F. Hisker, Y.-S. Kang, J.-S. Lee, and C. Herzig, Ag diffusion and interface segregation in nanocrystalline $\gamma$-FeNi alloy with a two-scale microstructure, Acta Mater. 52, 631 (2004b).