Novel stable structure of Li₃PS₄ predicted by evolutionary algorithm under high-pressure

Computation-aided materials searches have produced far better results than expected. For example, virtual screening of materials databases by using machine-learning techniques enables us to narrow the list of candidates possessing useful physical properties, resulting in the discovery of novel materials.^1–4 Another powerful tool of computational techniques is the use of an evolutionary algorithm, which is a field of machine learning, to predict stable crystal structures without any information in advance. About the difference between two calculation techniques, former needs the database, but latter does not need it. Latter technique can find the most stable structure, which has not been found. There have been many reports on the prediction by an evolutionary algorithm and discovery of new materials.^5–8 For example, although it is well known that the stable structure of NaCl has a composition ratio of Na:Cl = 1:1, it has been suggested from both calculations and experiments that NaCl has other stable structures with different composition ratios under high pressures over 10 GPa.⁵ 

The evolutionary algorithm is a calculation technique used to implement mechanisms suggested from biological evolution, such as reproduction, mutation, genetic modification, natural selection, and survival of the fittest. In genetic optimization problems, which apply natural selection theory to the calculation, individuals with a high fitness score, which is calculated according to the characteristics of individuals, produce many offspring. On the other hand, individuals with a low fitness score produce few offspring and die out. This procedure enables us to reach a solution gradually. When applying this algorithm to stable-structure searches, crystal parameters and formation energy can be considered as characteristics of individuals and a measure of fitness, respectively. To find the most stable structure is to solve the optimization problem of finding the global minimum of energy surface in a wide parameter space described by the crystal parameters. It is well known that the evolutionary algorithm works well for this kind of problem.

In this paper, we report one example of a materials search using the evolutionary algorithm. The system taken up in this study is designed for the development of materials for a solid-state electrolyte of lithium ion rechargeable batteries. Lithium ion batteries have high energy densities and high charge/discharge efficiencies. However, these batteries have safety problems because of their use of organic liquid electrolytes. Therefore, the development of all-solid batteries is matter of interest. Recently, it was reported that a new solid electrolyte, Li₁₀GeP₂S₁₂, showed a high ionic conductivity of 1.2 × 10⁻² S cm⁻¹, which is comparable to those of organic liquid electrolytes at room temperature.⁹ This finding has promoted new materials searches focused on lithium sulfide. For example, Li₃PS₄ is a simple ternary compound with three types of polymorphs α, β, and γ,¹⁰ and nanostructured β-Li₃PS₄ is of interest because of its high Li ion conductivity. Although the successive phase transitions of Li₃PS₄ at ambient pressures have been investigated, the high-pressure state of this compound has not been studied.

We found the new phases of Li₃PS₄ under high pressure, which is consistent to computational prediction. The evolutionary algorithm was employed to search for stable crystal structures under high pressure, where the formation energies of the structures representing the fitness of individuals were calculated by first-principles calculations. One new structure, δ-type, were predicted to appear under high pressure. Furthermore, in situ X-ray diffraction (XRD) under high pressure revealed the existence of δ-Li₃PS₄. Present results insist the computational techniques is usefull tool to search for new materials searches.

To predict which structures become stable under high-pressure conditions, we employed the evolutionary algorithm implemented in USPEX.^11–13 For the initial population, 30 random structures were generated in the first generation from Li atoms and PS₄ molecules. In order to decrease the total energy, structural relaxation was performed on each structure using the Vienna Ab initio Simulation Package^14,15 (VASP). From the stable structures in the top 60%, variation operators such as heredity, mutation, permutation, and so on were used to generate offspring. These calculations were repeated until the simulation reached the convergence criteria, where the best structure is that which did not change for five generations. Throughout the simulation, the ratio of Li, P, and S was kept at 3:1:4. The details of structure optimization and the total energy calculations with the VASP code are as follows. Exchange and correlation functions were given by the generalized gradient approximation as proposed by Perdew et al.¹⁵ Electron-ion interactions were represented by the projector augmented wave method with plane waves up to an energy of 420 eV.¹⁶ The Brillouin zone sampling in a primitive cell was performed by a Γ point-centered k-mesh, which was limited to 0.15 Å⁻¹.

For the structural optimization, the convergence criterion for electronic self-consistency and ionic relaxation loop were set to 10⁻⁶ eV and 0.01 eV, respectively. Gaussian smearing with the width of the smearing in 0.1 eV was used. For the total energy calculation, the convergence criterion for electronic self-consistency, an ionic relaxation loop of 10⁻⁶ eV, and the tetrahedron smearing method with Blöchl corrections were used. To investigate the dynamical stability of Li₃PS₄, we carried out first-principles phonon calculations with the finite displacement method^17,18 using the phonopy code¹⁹ combined with the VASP code.

In situ observation of the predicted crystal structure under high-pressure (HP) using XRD were conducted. The starting material for the HP experiment, polycrystalline γ-Li₃PS₄, was synthesized from stoichiometric mixtures of Li₂S (99.9%, MITSUWA) and P₂S₅ (99%, Aldrich). The mixtures were formed into pellets 10 mm in diameter, which were then sealed in an evacuated quartz tube and heated at 500 °C for 8 h. In situ energy-dispersive XRD measurements were carried out using the Kawai-type multi-anvil high pressure apparatus (SPEED-Mk.II-D) installed on the BL04B1 beamline at Spring-8 in Japan under the proposal No. 2016B1322. The samples were installed in boron nitride (BN) capsules. Further, Magnesia (MgO) and graphite were used as a pressure medium and heater, respectively. Temperatures were measured by W₉₇Re₃–W₇₅Re₂₅ thermocouples inserted into the high pressure cell. The pressures were determined by the unit-cell volume of MgO using equations of state of Jamieson et al.,²⁰ which was placed inside the high pressure cell. In situ XRD measurements were carried out the pressure (P) between 0.1 MPa and 20 GPa and the temperature (T) between 25 (R.T.) and 1000 °C. The incident white beam was reduced in size 50 × 100 μm² using two slits. The diffraction angle was fixed at 2θ = 5.00057°.

Starting from random structures, the evolutionary algorithm enables us to generate more stable structures step by step with increasing generation number. During the simulation, the number of atoms in the cell were fixed to 6, 4, and 8 for Li, P, and S, respectively. Figure 1 shows the evolution of crystal structures of Li₃PS₄ under various pressures. At each generation, 30 structures were generated by several variation operators. The lowest enthalpies for each generation are represented in this figure. At P = 0, the enthalpy decreased monotonically and reached its lowest value at the third generation. This stable crystal structure is shown in Fig. 2(a). The obtained lattice parameters were a = 7.818, b = 6.61, c = 6.204 Å, α = β = γ = 90°, with space group Pmn2₁. As already mentioned, Li₃PS₄ has three types of polymorphs—α, β, and γ.^10,21 Among them, γ-Li₃PS₄ is a stable crystal structure at low temperature. The lattice parameters of the experimentally reported γ-Li₃PS₄, shown in Fig. 2(b), were a = 7.708, b = 6.535, and c = 6.136 Å, with space group Pmn2₁ at 297 K. In the Table I, experimental crystal parameters can be compared to those of predicted structure in this study. Except for the small amount of deviation in the lattice parameters and atomic positions, the predicted stable structure showed many common structural characteristics with γ-Li₃PS₄, such as space group and atomic sites. The results strongly indicate that it is possible to predict stable structures in Li₃PS₄ by using the evolutionary algorithm. As far as the Li position concerned, these two structures are slightly different. Obtained total energy as a ground state in our calculation was -4.3940 eV/atom. On the other hand, the estimated total energy for the experimentally reported structure is -4.3928 eV/atom, which is slightly higher than the predicted one. This discrepancy may occur due to the difficulty of the determination of the position of the light weight Li ion in the experiment.

Under finite pressures, the evolution process also seems to successfully progress because the enthalpy decreases monotonically. Especially at P = 10–20 GPa, stable crystal structures were found very quickly. At 5–20 GPa, the obtained stable crystal structures were basically identical except for the unit cell volume. The obtained lattice parameters were a = 5.807 and c = 6.892 Å, with space group P42₁c(114) at 10 GPa. This crystal structure is shown in Fig. 2(c) and crystal parameters are listed in Table I. In this paper, the obtained new structures at 5–20 GPa is denoted as delta (δ) phases.

Stability of δ-Li₃PS₄ was examined in terms of vibration behavior, the so-called phonon dispersion relation. First-principles phonon calculations with the finite displacement method^17,18 were carried out. For the present calculation, a 2 × 2 × 2 supercell containing 128 atoms was used and the atomic displacement distance was set to 0.01 Å. The supercell was created from the fully optimized primitive cell of δ-Li₃PS₄ at P = 10 and 20 GPa. Figures 3(a) and (b) show the phonon dispersion relation for δ-Li₃PS₄ at P = 10 GPa, and at P = 20 GPa, respectively. The vertical axis shows the frequency and the horizontal axis shows the wave vector. The wave vectors of the phonon dispersion of δ-Li₃PS₄ is shown along the high-symmetry directions Γ-M-X-Γ-Z-A-R-Z.

The phonon spectrum for all structures is distributed below 20 THz, and it shows three acoustic and many optical branches overlapping in a complex manner. We focused on a softening of the acoustic modes, which show an imaginary number of the frequency and incompatibility with the dynamic stability of structure. There is no part where the frequency is less than zero for δ-Li₃PS₄. Furthermore, the eigenstate of the phonons at P = 20 GPa is distributed in a more high-energy region as compared to the spectra at P = 10 GPa. This is a reasonable variation that is expected under high pressures because the lattice contraction caused by the high pressure strengthens the chemical bonding, which increases the phonon frequency. In more detail, there are three optical branches around 13, 15, and 17 THz at 10 GPa. For these branches, main contributions come from S, Li, and P atoms, respectively. Among of them, S-branch moves lower, on the other hand, Li- and P-branch moves to higher energy. As a result, S-branch seems to be integrated to the lower complex phonon branches. Although this behavior contains important information about the changes in bonding with pressure, more detailed analysis is necessary in future work. These calculation results show that δ-Li₃PS₄ is dynamically stable at both P = 10 GPa and 20 GPa.

The pressure-induced solid–solid phase transformation at T = 0 can be discussed with the enthalpy instead of the Gibbs free energy. Figure 4 shows the pressure dependence of the enthalpy for Li₃PS₄ structures, where reference state is γ-type structure. At the low-pressure region of 0 < P < 5 GPa, the order of the enthalpies of Li₃PS₄ is γ < β < δ; therefore, γ-Li₃PS₄ is the most stable structure. With increasing pressure, the enthalpy of δ-Li₃PS₄ becomes lower than those of γ-Li₃PS₄ and β-Li₃PS₄. At P = 10 GPa, the enthalpy of δ-Li₃PS₄ is the lowest value among of all the structures. From fig. 4, the transformation pressure is found to be approximately 5 GPa. Above this pressure, the δ-Li₃PS₄ structure becomes relatively stable as compared to the other structures. This pressure dependence of enthalpy, in addition to the structure prediction using the evolutionary algorithm and the phonon calculation, suggests that δ-Li₃PS₄ can be obtained at pressures greater than 5 GPa.

In situ observation of the predicted δ-Li₃PS₄ were conducted under high-pressure (HP) using X-ray diffraction (XRD). A polycrystalline sample of γ-Li₃PS₄ was used as the starting material for the HP experiment. Before the HP experiments, the quality of the powdered samples was examined from angle dispersive XRD patterns. Figure 5 shows the XRD patterns of the sintered sample, which showed that the main phase of the sample was γ-Li₃PS₄ and that a very small amount (∼1 %) of Li₂PS₃ was present as an impurity phase.

The above prediction about the existence of a new phase was confirmed by experiments. Figure 6 shows in situ XRD patterns for Li₃PS₄. The XRD pattern of the sample at ambient pressure, shown in Fig. 6(a). We note that several weak peaks can be ascribed to the BN container and unknown impurities, and other Bragg peaks could be indexed to γ-Li₃PS₄. However crystal parameters were slightly changed and relative intensity was different with the those of γ-Li₃PS₄. Even at the ambient pressure, the sample was affected from the small pressure during the sample preparation, and the crystal structure may be distorted even in this situation. Structural change under the low pressure region must be investigated in future study. Figure 6(b) shows the XRD pattern at P = 5 GPa at R.T., 100, 200, and 300 °C, where the sample was heated after applying pressure at R.T. The diffraction patterns totally changed at R.T. as compared to those shown Fig. 6(a). At P = 5 GPa, these Bragg peaks can be indexed to a δ-Li₃PS₄ unit cell. The main peak around d = 2.7 Å and several other peaks are matched the calculated patterns of δ-Li₃PS₄. With increasing temperature up to 300 °C, peaks for BN and Li₃PO₄ appeared in the XRD patterns. From the XRD data, we obtained the lattice parameters of the unit cell for the structures. The obtained lattice parameters were a = 5.770 and c = 6.947 Å at R.T. and a = 5.852 and c = 6.863 Å at 100 °C, respectively. There is good agreement between the experimental and theoretical values of the equation of a = 5.807 and c = 6.892 Å. Furthermore, we found that δ-Li₃PS₄ phases were transformed from γ-Li₃PS₄ at P = 5 GPa, which also shows good agreement with the theoretical predictions. These results confirm that δ-Li₃PS₄ exists at T < 300 °C and at P = 5 GPa.

For δ-Li₃PS₄, several physical properties, including Li ion conductivity, are also attractive. It was reported that bulk γ-Li₃PS₄ showed a low ionic conductivity of 3 × 10⁻⁷ S cm⁻¹ at R.T.; on the other hand, nanoporous β-Li₃PS₄ had an extremely high ionic conductivity of 1.6 × 10⁻⁴ S cm⁻¹ at R.T.⁴ For the conductivity measurements, a sufficient amount of the bulk sample is necessary. We tried to synthesize δ-Li₃PS₄ as a single phase several times, but the attempts ended in failure until now. The XRD patterns of a sample synthesized under high pressure are understood to be mixtures of δ-Li₃PS₄, β-Li₃PS₄, Li₂PS₃, and Li₃PO₄. Especially, we note that this sample was easily oxidized. The physical properties of δ-Li₃PS₄ are not clear. Wang et al.²² revealed a fundamental relationship between anion packing and ionic transport in fast Li-conducting materials. A body-centered cubic-like anion framework allows direct Li hops between adjacent tetrahedral sites, which is most desirable for achieving high ionic conductivity. On the basis of these findings, there is a serious disadvantage in δ-Li₃PS₄ because it is synthesized under pressure. However, in situ measurement of the physical properties of δ-Li₃PS₄ is necessary, which will be addressed in a future work.

Finally, a new solid electrolyte, Li₁₀GeP₂S₁₂, was discovered as the highest ionic conductor among all-solid electrolytes.² It was very difficult to explore other new solid electrolytes only from the experimental techniques because some dangerous chemicals are included in this system. As shown in this study, it is more efficient in terms of new materials searches to estimate some experimental parameters, such as temperature, pressure, and composition in advance. To develop a combination of computational and experimental techniques is very useful for materials searches.

By combining computational techniques and HP synthesis for Li₃PS₄, new phases were successfully found under high pressure. The evolutionary algorithm was employed to search for stable crystal structures, where the formation energy of the structures as representing the fitness of individuals was calculated by first-principles calculations. Following the theoretical prediction for δ-Li₃PS₄, in situ synchrotron XRD measurements revealed the existence of δ-Li₃PS₄. This new δ-Li₃PS₄ structure was stable at pressures over 5 GPa. The combination of computational prediction and experimental observation works well for the new materials searches.

The authors gratefully acknowledge Prof. Tsuneyuki, and Prof. Nakamura for the valuable discussion. The synchrotron radiation experiments were performed at the BL04B1 of Spring-8 with the approval of the Japan Synchrotron Radiation Research Institute (JASRI) (Proposal No. 2016B1322, and 2017A1273).