Robust and efficient calculation of activation energy by automated path search and density functional theory

Because inorganic solid electrolytes are one of the key components for application to all-solid-state batteries, high-ionic-conductivity materials must be developed. Therefore we propose a method of efficiently evaluating the activation energy of ionic diffusion by calculating a potential-energy surface (PES), searching for the optimal diffusion path by an algorithm developed using dynamic programming (DP), and calculating the corresponding activation energy by the nudged elastic band (NEB) method. Taking $\beta \text{−}{\mathrm{Li}}_3{\mathrm{PS}}_4$ as an example, the activation energy of Li-ion diffusion was calculated as 0.43, 0.25, and 0.40 eV in the $a$-, $b$-, and $c$-axis directions, respectively, which is in good agreement with previously reported values. By comprehensively searching for the lowest energy path by PES-DP, the arbitrariness of the path selection can be eliminated, and the activation energy must only be calculated using the NEB method a few times, which greatly reduces the computational cost required for evaluating activation energy and enables the high-throughput screening of solid-state electrolytes.

Figure 1

(a) Crystal structure of $\beta \text{−}{\mathrm{Li}}_3{\mathrm{PS}}_4$ rendered by vesta [28]. Li, P, and S atoms and ${\mathrm{PS}}_4$ polyhedra are green, purple, yellow, and purple, respectively. (b) Asymmetric unit in unit cell. Grid points in asymmetric unit are represented by black spheres.

Figure 2

Reduction in number of grid points. Crosses represent [P or S] (orange) and Li (green) sites. Gray circles represent grid points from which energy is calculated. White circles are located within ionic radii of P or S sites and are removed from calculation target.

Figure 3

Voronoi tessellation of Li sublattice and PES in Voronoi region. Li sites at Voronoi center correspond to local minimum of PES. When additional Li ion is introduced at grid point in Voronoi region, one Li ion at center of Voronoi region is removed. (a, b) Ground state, (c)–(f) hopping state. After crossing Voronoi region, vacancy moves to adjacent Li site because Li site nearest to grid point changes.

Figure 4

Lowest-activation-energy diffusion paths calculated by DP in each crystal-axis direction. Diffusion paths to (a) $b$, (b) $c$, and (c) $a$ axes are shown by black spheres, and white spheres represent PES global minimum.

Figure 5

PES along $b$-axis diffusion path. Dotted line corresponds to grid point nearest to Li site.

Figure 6

Energy profiles of diffusion path along each crystal axis as calculated by NEB method. Reaction coordinates are normalized from 0 to 1, and red circles correspond to each NEB image. Dashed line indicates energy curve obtained by cubic spline interpolation.

Figure 7

Diffusion paths calculated by AIMD and NEB methods. Yellow surface indicates region showing high probability density of Li ions, as calculated by AIMD, and black spheres show NEB-optimized diffusion paths.

Figure 8

(a) Two diffusion paths between $8d$ Li sites. Path obtained by DP is $8d\text{−}4c\text{−}8d$ (path 1), and path connecting shortest distance is $8d\text{−}8d$ (path 2). (b) Energy profiles for Li-ion diffusions of paths 1 and 2.

Figure 9

Relationship between number of calculation points required and ionic-radius scale factor. Grid points within scaled-ionic-radius (i.e., ionic radius multiplied by abscissa) distance from atomic coordinates of P or S ions are omitted from calculations.

Figure 10

(a) Diffusion path of ${\mathrm{Li}}_3{\mathrm{YBr}}_6$ showing lowest activation energy and (b) corresponding PES. (c) Energy profile of diffusion path calculated using NEB method. O and T indicate octahedral and tetrahedral sites, respectively. Reaction coordinates from 0 to 1 in panel (c) correspond to path coordinates from 0 to 0.5 in panel (b).