Nonequilibrium molecular dynamics for accelerated computation of ion–ion correlated conductivity beyond Nernst–Einstein limitation

Abstract

Condensed matters with high ionic conductivities are crucial in various solid devices such as solid-state batteries. The conduction is characterized by the cooperative ionic motion associated with the high carrier density. However, the high cost of computing correlated ionic conductivities has forced almost all ab initio molecular dynamics (MD) to rely on the Nernst–Einstein dilute-solution approximation, which ignores the cross-correlation effect. Here we develop a chemical color-diffusion nonequilibrium MD (CCD-NEMD) method, which enables to calculate the correlated conductivities with fewer sampling steps than the conventional MD. This CCD-NEMD is demonstrated to well evaluate the conductivities in the representative solid electrolyte bulk Li₁₀GeP₂S₁₂ and Li₇La₃Zr₂O₁₂. We also applied CCD-NEMD to the grain boundary of Li₇La₃Zr₂O₁₂ and demonstrated its applicability for calculating interfacial local conductivities, which is essential for investigating grain boundaries and composite electrolytes. CCD-NEMD can provide further accurate understanding of ionics with ionic correlations and promote developing solid devices.

Introduction

Understanding ionic transport phenomena in solid electrolytes (SEs) is fundamentally important for improving the performance of various solid devices, for example, all-solid-state batteries, solid oxide fuel cells, and sensors¹. In particular, Li-ion all-solid-state batteries have attracted considerable attention as next-generation batteries owing to their safety concerns. Finding SEs with a high ionic conductivity is crucial for further improvement, and atomistic simulations are required to enable accurate and fast estimation of ionic conductivities for the material design.

Currently, ab initio molecular dynamics (AIMD) is the most effective tool for calculating the conductivity as they are free from fitting parameters. Conductivity can be calculated via equilibrium molecular dynamics (EMD) using the following equation^2,3 (the left-hand side in Fig. 1a):

$$\begin{array}{lll}\\ \sigma = \frac{{e^2}}{{dVk_{{{\mathrm{B}}}}T}}\mathop {\int}\limits_0^\infty {\mathop {\sum}\limits_{i,j} {z_iz_j\left\langle {{{{\mathbf{v}}}}_i\left( t \right) \cdot {{{\mathbf{v}}}}_j\left( 0 \right)} \right\rangle {{{\mathrm{d}}}}t} } \cr \\ = \frac{{e^2}}{{Vk_{{{\mathrm{B}}}}T}}\mathop {\sum}\limits_i {z_i^2D_{{{{\mathrm{self}}}},i}} + \frac{{e^2}}{{dVk_{{{\mathrm{B}}}}T}}\mathop {\sum}\limits_{i,j \ne i} {\mathop {\int}\limits_0^\infty {z_iz_j\left\langle {{{{\mathbf{v}}}}_i\left( t \right) \cdot {{{\mathbf{v}}}}_j\left( 0 \right)} \right\rangle {{{\mathrm{d}}}}t} } \cr \\ = \mathop {{{{{\mathrm{lim}}}}}}\limits_{t \to \infty } \frac{{e^2}}{{2dtVk_{{{\mathrm{B}}}}T}}\left[ {\mathop {\sum}\limits_i {z_i^2\left\langle {\left[ {{{{\mathbf{x}}}}_i\left( t \right) - {{{\mathbf{x}}}}_i\left( 0 \right)} \right]^2} \right\rangle } + \mathop {\sum}\limits_{i,j \ne i} {z_iz_j\left\langle {\left[ {{{{\mathbf{x}}}}_i\left( t \right) - {{{\mathbf{x}}}}_i\left( 0 \right)} \right] \cdot \left[ {{{{\mathbf{x}}}}_j\left( t \right) - {{{\mathbf{x}}}}_j\left( 0 \right)} \right]} \right\rangle } } \right]\cr \\ = \sigma _ + ^{{{{\mathrm{self}}}}} + \sigma _ - ^{{{{\mathrm{self}}}}} + \sigma _{ + + }^{d{{{\mathrm{istinct}}}}} + \sigma _{ - - }^{{{{\mathrm{distinct}}}}} + 2\sigma _{ + - }^{{{{\mathrm{distinct}}}}}\\ \end{array}$$

(1)

where e, V, k_B, and T are the elementary charge, system volume, Boltzmann constant, and temperature, respectively. The z_i, x_i, and v_i are the charge valency, position, and velocity of particle i, respectively; d is the dimension of x_i and v_i, respectively. Detailed definition of D_self, $\sigma _ + ^{{{{\mathrm{self}}}}}$, $\sigma _ - ^{{{{\mathrm{self}}}}}$, $\sigma _{ + + }^{{{{\mathrm{distinct}}}}}$, $\sigma _{ - - }^{{{{\mathrm{distinct}}}}}$, and $\sigma _{ + - }^{{{{\mathrm{distinct}}}}}$ are provided in Supplementary Note 1. Equation (1) requires a long convergence time in terms of the AIMD timescale. To avoid this uncertainty, Nernst–Einstein (N–E) dilute-solution approximation has been often adopted, which neglects the carrier’s distinct correlated terms (i ≠ j; $\sigma _{ + + }^{{{{\mathrm{distinct}}}}}$, $\sigma _{ - - }^{{{{\mathrm{distinct}}}}}$, and $\sigma _{ + - }^{{{{\mathrm{distinct}}}}}$). The N–E approximated conductivity σ_dilute is related to the self-diffusion coefficient D_self via the equation $\sigma _{{{{\mathrm{dilute}}}}} = \sigma _ + ^{{{{\mathrm{self}}}}} + \sigma _ - ^{{{{\mathrm{self}}}}}$. The merit of N–E approximation is that the error is scaled by a factor of the square root of the number of particles, N^1/2, because all carriers independently contribute to the conductivity. However, this assumption relies on the infinite dilution limit, which does not hold for SEs with their high carrier density. Therefore, methodologies that allow for both efficient sampling and inclusion of a distinct ion–ion correlation effect (〈v_i(t) · v_j(0)〉 in Eq. (1)) are required.

Fig. 1: Schematic representation of ionic conductivity calculations in conventional MD and CCD-NEMD.

a Conventional MD: equilibrium MD (left) and color-diffusion NEMD (right) and b CCD-NEMD: the color charge is defined on the basis of the charge valency of the chemical units. For example, the color charges c_i of Li⁺ and the X^x–-unit are denoted as +1 and –X, respectively. The conductivity considering the ion–ion distinct correlation effect can be calculated faster using CCD-NEMD than EMD.

We focused on nonequilibrium molecular dynamics (NEMD) to enhance the sampling efficiency. The color-diffusion algorithm-based nonequilibrium molecular dynamics (CD-NEMD)^4,5 has been applied^{6,7,8,9,10,11} to accelerate the calculations of N–E approximated conductivity. CD-NEMD defines a fictitious charge called color charge, which does not interact with the other color charges, and imposes a fictitious external field called color field on these charges. The color field accelerates rare-event ion hopping in SEs, which enhances the sampling efficiency. However, the assignment of color charges in the previous studies on SEs is problematic. For instance, in Li-ion SEs, the color charges of half Li-ions and the other half are assigned as +1 and –1, respectively, to satisfy the color-charge neutrality^6,7,10. This method uses N–E approximation to estimate the conductivity, which prevents the calculation of ion–ion correlated conductivity. While CD-NEMD can enhance the sampling efficiency, there is room for improvement to avoid the reliance on N–E approximation.

Here, we extended the CD-NEMD to calculate the conductivities, including the ion–ion distinct correlation beyond the N–E approximation, namely “chemical” CD-NEMD (CCD-NEMD). We demonstrated that CCD-NEMD, where the color charge is equal to the charge valency of the ionic species, provides an ionic conductivity including the distinct correlation effect. We chose Li₁₀GeP₂S₁₂ (LGPS) as a validation of CCD-NEMD because LGPS is one of the promising SEs with a high experimental bulk ionic conductivity^12,13,14,15. The calculations were performed by AIMD. The results are used to demonstrate that CCD-NEMD can calculate ionic conductivities effect with smaller statistical errors than conventional EMD. Furthermore, we proposed its applicability to interfacial local conduction using a grain-boundary structure of the promising oxide-type Li-ion SE, namely, Li₇La₃Zr₂O₁₂ (LLZO).

Results and discussion

Overview of NEMD

The CCD-NEMD method determines the color charges based on charge valency of chemical units, whereas the previous CD-NEMD assigns values of +1 or –1 to carrier ions (the right-hand side in Fig. 1a). For instance, in the case of LGPS comprising Li⁺, PS₄^3–, and GeS₄^4– chemical units, the color charges c_i of CCD-NEMD are assigned as c_i∈Li = +1, ∑_i∈PS4 c_i = –3, and ∑_i∈GeS4 c_i = –4 (Fig. 1b), whereas those of CD-NEMD are assigned as c_i∈Li = (–1)ⁱ and c_i∈thers = 0. It should be noted that color-charge neutrality Σc_i = 0 must be satisfied for total momentum conservation. The color charge c_i exhibits a force c_iF_e from the external color field F_e, which induces the steady flux of color charge. Thus, we can obtain the ionic conductivity from the flux and F_e based on the linear response theory in CCD-NEMD. Furthermore, CCD-NEMD overcomes the following critical drawbacks of CD-NEMD (the right-hand side in Fig. 1a): (i) dependence on N–E dilute-solution approximation, which neglects significant ion–ion correlation in SEs, (ii) unavailability for SEs with one-dimensional narrow conduction channels.

For the validation of computing bulk ion–ion correlated ionic conductivity, we performed AIMD on the LGPS and classical MD on the LLZO solid electrolyte. In addition, we constructed the Σ3(112) grain boundary of LLZO to discuss the applicability of CCD-NEMD to interfacial conduction.

Issues of conventional CD-NEMD

The equations of motion in CD-NEMD after applying the color field F_e along the z-axis to the color charge c_i of the tagged particles are as follows^4,5:

$${{{\dot{\mathbf q}}}}_i = \frac{{{{{\mathbf{p}}}}_i}}{{m_i}}$$

(2)

$${{{\dot{\mathbf p}}}}_i = \left\{ {\begin{array}{*{20}{c}} {F_{i,\{ x\,or\,y\} } - \alpha p_{i,\{ x\,or\,y\} }} \\ {F_{i,z} + c_iF_{{{\mathrm{e}}}}} \end{array}} \right.$$

(3)

where α is the coupling coefficient of the thermostat to remove the energy originating from the color field and q_i, p_i, and F_i are the position, momentum, force via DFT (density functional theory), respectively. The thermostat is not applied to the z-direction to keep the external force only F_e. If the thermostat is applied to the z-direction, the effective external force will become c_iF_e – αp_i. Color-charge neutrality Σc_i = 0 must be satisfied for total momentum conservation. In previous studies of SEs, half of the carrier ions were given +1 color charge and the other half were given −1 charge like c_i = (–1)^{i 6,7,10}. The color-charge flux driven by the color field F_e is as follows:

$$J_{{{\mathrm{c}}}} = \frac{1}{V}\mathop {\sum}\limits_i {c_iv_i\left( t \right)} = \frac{1}{V}\mathop {\sum}\limits_i {\left( { - 1} \right)^iv_i\left( t \right)}$$

(4)

CD-NEMD provides the self-diffusion coefficient from the J_c(t) and F_e^5,16:

$$D_{{{{\mathrm{self}}}}} = \frac{{k_{{{\mathrm{B}}}}T}}{{(N/V)}}\mathop {{{{{\mathrm{lim}}}}}}\limits_{t \to \infty } \mathop {{{{{\mathrm{lim}}}}}}\limits_{F_{{{\mathrm{e}}}} \to 0} \frac{{\left\langle {J_c\left( t \right)} \right\rangle _{{{{\mathrm{NEMD}}}}}}}{{F_{{{\mathrm{e}}}}}}$$

(5)

Therefore, CD-NEMD requires D_self of the carrier ion to estimate conductivity and N–E approximation must be used, which prevents computing ionic conductivities with the ion–ion correlation (Problem-1 on the right-hand side in Fig. 1a). In addition, CD-NEMD is expected to have limited application in SEs with a narrow one-dimensional channel. Problem-2 on the right-hand side in Fig. 1a shows an example of ill definition of color charges. When the sum of all color charges within each conduction channel is zero and ions cannot pass each other, a net color flux cannot be generated. To avoid the loss of conduction of carrier ions, some researchers set all the color charges to +1 and artificially removed the momentum of the non-diffusing crystalline framework to prevent the drift of the entire system^8,9,11. However, this procedure may underestimate the cooperative conduction with the framework like paddle-wheel motion^17,18, which is essential in some SEs.

CCD-NEMD for computing conductivities with ion–ion correlation

We reformulated the conventional CD-NEMD to avoid both the reliance on N–E approximation and restricted applicability against SEs with narrow one-dimensional conduction channel. According to the linear response theory, the color-charge flux is expressed as⁵:

$$\begin{array}{ll}\\ \left\langle {J_{{{\mathrm{c}}}}\left( t \right)} \right\rangle _{{{{\mathrm{NEMD}}}}}\, =\, \mathop {{\lim }}\limits_{F_e \to 0} \frac{{VF_{{{\mathrm{e}}}}}}{{k_{{{\mathrm{B}}}}T}}\mathop {\int}\limits_0^t {\left\langle {J_{{{\mathrm{c}}}}\left( s \right)J_{{{\mathrm{c}}}}\left( 0 \right)} \right\rangle _{F_e = 0}{{{\mathrm{d}}}}s} \cr \\ \qquad\qquad\qquad=\, \mathop {{{{{\mathrm{lim}}}}}}\limits_{F_e \to 0} \frac{{F_{{{\mathrm{e}}}}}}{{Vk_{{{\mathrm{B}}}}T}}\left[ {\mathop {\sum}\limits_i {c_i^2} \mathop {\int}\limits_0^t {\left\langle {v_i\left( {{{\mathrm{s}}}} \right)v_i\left( 0 \right)} \right\rangle _{F_e = 0}{{{\mathrm{d}}}}s} + \mathop {\sum}\limits_{i \ne j} {\mathop {\int}\limits_0^t {c_ic_j\left\langle {v_i\left( {{{\mathrm{s}}}} \right)v_j\left( 0 \right)} \right\rangle _{F_e = 0}{{{\mathrm{d}}}}s} } } \right]\\ \end{array}$$

(6)

Although the sum of the cross term ∑_i≠j c_ic_j〈v_i(t)v_j(0)〉 becomes zero in the conventional CD-NEMD (c_i = (–1)ⁱ)^5,16, calculations of ionic conductivities must include the cross term to be consistent with the Green–Kubo formulation (Eq. (1)). Comparing between Eqs. (1) and (6), we propose chemical color-diffusion NEMD (CCD-NEMD), whose color charges c_i are determined as the charge valency z_i of chemical units. The color-charge flux is as follows:

$$J_c = \frac{1}{V}\mathop {\sum}\limits_i {c_iv_i\left( t \right)} = \frac{1}{V}\mathop {\sum}\limits_i {z_iv_i\left( t \right)}$$

(7)

Using the Eqs. (1), (6), and (7), we can obtain the ionic conductivity including ion–ion correlation term as:

$$\sigma = \mathop {{\lim }}\limits_{t \to \infty } \mathop {{\lim }}\limits_{F_e \to 0} \frac{{e^2\left\langle {J_c\left( t \right)} \right\rangle _{{{{\mathrm{NEMD}}}}}}}{{F_e}}$$

(8)

It should be noted that this definition of the color charge in CCD-NEMD satisfies the total momentum and is free from the issues arising from narrow one-dimensional conduction channels in CD-NEMD, as shown on the right-hand side in Fig. 1a.

We then summarize the various advantages of CCD-NEMD. Compared to CD-NEMD, CCD-NEMD can calculate the exact conductivity including the ion–ion correlation effect, which cannot be considered in CD-NEMD (see Supplementary Figs. 4 and 5 and Table 2 in Supplementary Note 5). It should be noted that in principle, CCD-NEMD can also be employed for a liquid system. Highly concentrated electrolytes will be meaningful examples because they exhibit strong correlation motions¹⁹.

We also describe the advantages of CCD-NEMD compared with the electric field NEMD. In classical MD, the electric field interacts with the point charges of atoms. If the point charge is the same as the charge valency^20,21, the conductivity can be calculated using the electric field strength E and ion current J_ion as J_ion/E. However, many classical force-field methods reduce the point charge from the charge valency to effectively consider the polarization effect. In this case, we cannot directly calculate the conductivity using the relationship σ = J_ion/E. By contrast, CCD-NEMD uses a fictitious charge based on the charge valency, which allows for the use of any force-field to calculate conductivities directly. AIMD with an electric field requires Berry-phase treatment²², which is more difficult to implement it. In terms of MLP-based MD, which have recently attracted much interest^23,24, the electric field is often meaningless because a lot of MLPs has no electric interaction term for energy calculations. By contrast, CCD-NEMD using fictitious charge and field can also be applied to MLP-based MD. Although some state-of-the-art MLPs include the calculation of point charges, the type of point charges should correspond to the Born effective charges to precisely calculate the conductivity using electric field^20,21. CCD-NEMD can be applied to a wider range of computational schemes than MD with electric field. In addition, this method can be easily implemented to any currently existing MD codes with minor modifications (Supplementary Note 2).

Relationship between Onsager theory and CCD-NEMD

We also discuss the connection between Onsager theory and CCD-NEMD. The total current density J is expressed by the Onsager transport equation^2,25,26:

$$J = \frac{1}{V}\mathop {\sum}\limits_i {z_iev_i}$$

(9)

$$J = - \mathop {\sum}\limits_{i,j} {z_ieL^{ij}\frac{{\partial \mu _j}}{{\partial x}}} - \mathop {\sum}\limits_{i,j} {z_iz_je^2L^{ij}\frac{{\partial \phi }}{{\partial x}}}$$

(10)

$$L^{ij} = \frac{V}{{k_{{{\mathrm{B}}}}T}}\mathop {\int}\limits_0^\infty {\left\langle {v_i(t)v_j(0)} \right\rangle dt}$$

(11)

where μ_j is the chemical potential, ϕ is the electrical potential, ρ_k is the concentration of species k, and L^ij is the Onsager transport coefficient. In the CCD-NEMD situation, the gradient of the chemical potential (∂μ_j/∂x) = 0, and − (∂ϕ/∂x) = F_e, and z_i = c_i, which gives us

$$J\mathop { \longrightarrow }\limits^{{{{\mathrm{CCD}}}} - {{{\mathrm{NEMD}}}}} \mathop {\sum}\limits_{i,j} {c_ic_je^2L^{ij}F_{{{\mathrm{e}}}}}$$

(12)

The relationship between the ionic conductivity σ and L^ij leads to

$$\sigma = \mathop {\sum}\limits_{i,j} {z_iz_je^2L^{ij}\mathop { \longrightarrow }\limits^{{{{\mathrm{CCD}}}} - {{{\mathrm{NEMD}}}}} e^2\left\langle J \right\rangle _{{{{\mathrm{NEMD}}}}}/F_{{{\mathrm{e}}}}}$$

(13)

Ionic migration in LGPS

We demonstrate the CCD-NEMD method on a realistic system, LGPS (Fig. 2a, b). The Li-ion conduction in LGPS has been shown to be anisotropic and higher in the c-axis than in the ab-plane ([110] direction) by the previous experimental^27,28 and computational²⁹ studies. The LGPS consists of an immobile framework of GeS₄, PS₄ tetrahedra and LiS₆ octahedra. The immobile framework constructs a fast and moderate Li-ion conduction channel along the c-axis and in the ab-plane, respectively, as shown in the isosurface of the probability density distributions for Li-ions (Fig. 2c). The c-axis conduction involves cooperative hopping, which indicates the importance of including the correlation effect of the carrier as already shown in the previous EMD study³⁰. To determine the diffusion correlation, we analyzed the distinct part of the van Hove correlation function G_d³¹:

$$G_d = \frac{1}{{4\pi r^2\rho N}}\left\langle {\mathop {\sum}\limits_{i \ne j} {\delta \left( {r - \left| {{{{\mathbf{r}}}}_j(t) - {{{\mathbf{r}}}}_i(0)} \right|} \right)} } \right\rangle$$

(14)

where δ(·) is the Dirac delta function, r is the Li–Li distance, ρ is the average number density of Li-ions, and N is the number of Li-ions. Our calculated G_d exhibits a pronounced peak at r = 0, which clearly proves correlated diffusion of Li-ion because the peak indicates that the position of the reference Li-ions is occupied by other Li-ions over an extremely short duration. The ionic conduction in LGPS involves strong Li–Li correlation. Thus, it is a suitable system for testing our proposed CCD-NEMD method.

Fig. 2: Crystal structure and transport property of Li₁₀GeP₂S₁₂.

The equilibrated structure of 2 Li₁₀GeP₂S₁₂ (8.72 × 8.72 × 12.64 ų) obtained by NVT simulation at 1200 K: a side view and b top view. c Isosurface of Li-ion probability density distribution, and d the distinct part of the van Hove correlation function of the EMD as a function of Li–Li distance r and time t at 1200 K.

Validation and efficiency of CCD-NEMD in LGPS

To validate our method and evaluate its efficiency, we compared the N–E approximated (Eq. (30) in Supplementary Note 1) and non-approximated (Eq. (46) in Supplementary Note 3) conductivities using EMD and one of the CCD-NEMD, denoted as σ_dilute, σ_EMD, and σ_NEMD, respectively. We show the cumulative averages of σ_dilute, σ_EMD, and σ_NEMD at 1200 K as a function of time in Fig. 3a to validate the effectiveness of the CCD-NEMD method and compare the computational costs. High-temperature simulation is suitable for the comparison because conductivities with a smaller statistical error can be obtained. The lines and error bars in Fig. 3a denote the average and standard deviation of the five independent samples with different initial velocities, respectively. As expected, σ_dilute is smaller than σ_EMD and σ_NEMD because N–E approximated, σ_dilute, neglects the contribution of the distinct ion–ion correlation. σ_NEMD and σ_EMD are in good agreement, which demonstrates that CCD-NEMD can appropriately reflect the contribution of ion–ion correlation. The convergence time is approximately 3.3 ps for σ_NEMD at a relative error of 20%, whereas those for σ_EMD and σ_dilute are approximately 28 ps and 10 ps, as shown by the arrows in Fig. 3b. The efficiency of CCD-NEMD is approximately nine times higher than that of conventional EMD. In addition, CCD-NEMD can calculate the conductivity with a small error and requires a shorter computational time than the conventional N–E-approximated EMD method. We performed the same validation using the LLZO electrolyte, which demonstrated that our method successfully calculated the correlated conductivity (see Supplementary Fig. 8 in Supplementary Note 6). These findings strongly support the generality of our method.

Fig. 3: Validation and efficiency of CCD-NEMD.

a Cumulative averaged conductivities (σ_dilute, σ_EMD, and σ_NEMD) and b the relative errors along c-axis at 1200 K. Each colored arrow corresponds to a position with a relative error of 20%. The lines and error bars indicate the average and standard deviation of the five independent samples with different initial velocities, respectively.

The temperature dependences of these conductivities in the c-axis direction and ab-plane are plotted in Fig. 4a, b, respectively. It should be noted that the calculations conducted to investigate the temperature dependence were performed using a smaller system than is typically used, which underestimates the conductivities in the ab-plane by 30% (see Supplementary Note 3). σ_EMD and σ_NEMD are in good agreement in the entire temperature region both in the c-axis and ab-plane. It should be noted that CCD-NEMD can obtain conductivities with smaller standard deviation than EMD. The conductivity in the c-axis was higher than that in the ab-plane as shown in previous studies^27,28,29. The activation energies for c-axis conduction of σ_EMD and σ_NEMD are 0.21 ± 0.03 and 0.17 ± 0.02 eV, ab-plane of σ_EMD and σ_NEMD are 0.35 ± 0.08 and 0.35 ± 0.03, respectively. These activation energies were also in good agreement between the values obtained from EMD and NEMD.

Fig. 4: Temperature dependence and anisotropy of ionic conductivity.

Arrhenius plot of the σT of LGPS a along c-axis and b in ab-plane. The values and error bars are the average and standard deviation of the five independent samples with different initial velocities, respectively.

We also plotted the Haven ratio H_R = σ_dilute/σ_NEMD as a function of temperature in Fig. 5. The CCD-NEMD clarifies the difference in H_R between the c-axis and ab-plane conductions. The H_R in the c-axis is smaller than that in the ab-plane, which means that the c-axis conduction includes more correlated motions that contribute to the conductivity. On the other hand, the experimental H_R at room temperature was approximately 0.5–0.7, which is expected for simple hopping mechanism¹³. A comparison between the experimental and computational H_R suggests that the conduction pathway in the ab-plane is more dominant (H_R ~ 0.69) than that in c-axis (H_R ~ 0.35) in the experimental situation. This result supports that the ab-plane with slow conduction could act as the bottleneck for long-range conduction in an experimental condition due to the channel blocking defects in the c-axis pathway as suggested in previous research^13,29.

Fig. 5: Haven ratio H_R = σ_dilute/σ_NEMD of LGPS along c-axis and in ab-plane as a function of temperature, respectively.

The values and error bars indicate the average and standard deviation of the five independent samples with different initial velocities, respectively.

CCD-NEMD for interfacial conduction: LLZO grain boundary

We discuss the applicability of CCD-NEMD for interfacial conduction in grain-boundary systems and composite electrolytes. NEMD can be used to calculate interfacial conductivities by defining the local flux in a straightforward manner, as proposed in previous studies^11,16. We constructed the energetically favorable Σ3(112) grain-boundary structure^32,33. Figure 6a, b show the relaxed structure of the Σ3(112) grain boundary. To calculate the local ionic conductivity in slab k, we define the local color current, J_c,k:

$$J_{{{{\mathrm{c}}}},k}\left( t \right) = \frac{1}{{V_{{{{\mathrm{slab}}}} - {{{\mathrm{k}}}}}}}\mathop {\sum}\limits_{i \in {{{\mathrm{slab}}}} - k} {c_iv_i\left( t \right)}$$

(15)

where V_slab–k is the volume of slab k. The local conductivity σ_slab-k can be determined using Eqs. (8) and (15):

$$\sigma _{{{{\mathrm{slab}}}} - k} = \mathop {{\lim }}\limits_{t \to \infty } \mathop {{\lim }}\limits_{F_{{{\mathrm{e}}}} \to 0} \frac{{e^2\left\langle {J_{{{{\mathrm{c}}}},k}\left( t \right)} \right\rangle _{{{{\mathrm{NEMD}}}}}}}{{F_{{{\mathrm{e}}}}}}$$

(16)

Fig. 6: Structures and grain-boundary ionic conductivity of the Σ3(112) grain boundary of LLZO.

Structures of the Σ3(112) grain boundary of LLZO depicting a all atoms and b only the Zr atoms. c Position-dependent ionic conductivity of the Σ3(112) grain-boundary structure of LLZO obtained via CCD-NEMD (σ) and CD-NEMD (σ_dilute). d Position-dependent Haven ratio (H_R). The values and error bars indicate the average and standard deviation of the five independent samples with different initial velocities.

The position-dependent ionic conductivities along the in-plane direction to the interface at 1000 K are presented in Fig. 6c. The position-dependent N–E approximated values (σ_dilute) are also shown in Fig. 6c. These conductivities decrease near the grain-boundary region. The AIMD calculation has suggested that the σ_dilute in the grain boundary was almost same as that in the bulk³³. On the other hand, the decreasing tendency around the grain boundary observed in the present results is consistent with the position-dependent σ_dilute value determined by previous EMD research using a different force-field³². It should be noted that in the previous study^32,33, the Li-ions were assigned to a slab based on the initial configuration, and the local σ_dilute values based on the initial configuration were calculated without considering the movement of the Li-ions between the slabs. The Li-ions in the grain boundary can move to the adjacent bulk-like slab, which could make result in the obtained grain-boundary diffusion coefficient being closer to the bulk coefficient than the actual value. While the conventional EMD approach suffers from ambiguity concerning the position-dependent diffusion coefficient and conductivity, the NEMD methods can robustly calculate these coefficients. We also calculated the position-dependent H_R values (Fig. 6d). In contrast to the case of conductivity, H_R did not significantly vary near the grain boundary.

In summary, we developed CCD-NEMD method, and successfully applied to predict the ionic conductivities with the ion–ion correlation effect, which cannot be included in the typically used N–E approximation. We performed the CCD-NEMD simulations for the ionic conductivity in a representative solid electrolyte, LGPS. The results show that CCD-NEMD leads to low computational cost with high statistical accuracy and is more efficient than the EMD using N–E approximation. Furthermore, in terms of the obtained ionic correlation, the experimental Li-ion long-range conduction is suggested to be dominated by the ab-plane slow conduction because the c-axis conduction could be blocked by defects. CCD-NEMD can be used to predict local conductivities in the interfacial region. This indicates that it can also be employed to estimate interfacial ionic conduction using the local flux, which is essential for enhanced conductivity in composites and grain-boundary resistance. The present CCD-NEMD would open a door to the more accurate and rapid MD research of ionics with ionic correlations.

Methods

All the simulations were performed using the CP2K package³⁴. The structures were drawn using VESTA³⁵ and VMD³⁶. Some analyses were conducted using Pymatgen^37,38.

Details of AIMD of LGPS

The experimental crystallographic positions consisting of 2 Li₁₀GeP₂S₁₂ (50 atoms) in 8.72 × 8.72 × 12.64 ų and 16 Li₁₀GePS₁₂ (400 atoms) in 17.44 × 17.44 × 25.28 ų were used as an initial structure. The smaller structure (50 atoms) was used to validate the temperature dependence of the conductivity because low-temperature EMD calculations require long trajectories. It should be noted that although the 50-atoms structure is smaller than the typical sized used in DFT calculation, it is adequate to validate the agreement between the EMD and NEMD results. In addition, the serious size effects of using such a small structure on the diffusive properties have not been observed in LGPS (see ref. ³⁰ and Supplementary Note 3). The 400-atoms structure, which is enough large for AIMD was used to evaluate the practical efficiency of our method. The forces F_i in Eq. (3) were calculated at the level of PBE-D3/DZVP-MOLOPT-SR-GTH^39,40,41 (400 Ry cutoff, Γ point) with the GTH-PBE pseudopotential⁴² using the Gaussian and planewave scheme⁴³ and CP2K/QUICKSTEP⁴⁴. The production runs were in the NVT ensemble using Nosé–Hoover thermostat^45,46 with a time step of 2 fs. The chain number of the thermostat in EMD was 3, and the thermostat was applied in the xyz direction. The chain number in NEMD was 1⁴⁷, and the thermostat was applied in directions perpendicular to the color field. After the equilibration for 10 ps, statistical averages were computed from trajectories (140–560 ps; EMD, 100–160 ps; NEMD. Details are presented in Supplementary Note 3).

The definitions of the color charges in CCD-NEMD were determined based on Li = +1, PS₄ = − 3, and GeS₄ = −4. In terms of the charge valency, we should assign P = + 5, Ge = +4, and S = –2. However, large color charges exhibit a strong force, which may destabilize the simulation. To address this concern, we divided the total charge valency of the chemical units and assigned it to each atom equally as color charge Li = +1, PS₄ = − 3 (P = S = − 0.6), and GeS₄ = −4 (G = S = − 0.8). The approach for assigning color charge for large chemical units would remain a factor to be considered. The strength of the color field was determined as 0.02 eV Å^−1 based on the linear response region, as shown in Supplementary Fig. 3 and Supplementary Note 3.

Details of classical MD for LLZO

We used the force-field optimized by Klenk and Lai, which well reproduced the Li-ion diffusivity^48,49. In the bulk cubic LLZO simulation (Supplementary Note 6), after the equilibration in the isotropic NpT ensemble (1000 K, 1 bar) using a Nosé–Hoover thermostat^45,46 and a barostat⁵⁰ for 10 ps, the equilibrated density was obtained during the following 50 ps operation. The production simulations of EMD and CCD-NEMD (F_e = 0.02 eV Å^−1 along the z-axis in Supplementary Fig. 7) were performed for 200 and 100 ps in the NVT ensemble (1000 K), respectively, using the equilibrated density obtained after the equilibration. The Σ3(112) grain-boundary model was constructed based on the Zr sublattice to preserve the ZrO₆ octahedron in the initial structure. The grain-boundary structure was relaxed following the same procedure as in the bulk simulation. The production simulation of CCD-NEMD (F_e = 0.02 eV Å^−1 along the y- or z-directions) was performed for 100 ps at 1000 K.