Synthesis challenges, thermodynamic stability, and growth kinetics of La–Si–P ternary compounds

Abstract

Although many new compounds have been recently predicted with the help of machine learning, the successful experimental synthesis of these compounds remains challenging. Computational insights about the thermodynamic stability and phase formation kinetics among the ground state and competing metastable phases are highly desirable to rationalize and attempt to overcome synthesis challenges experimentally. In this work, we explore synthetic challenges within ternary La–Si–P compounds through feedback between experimental and computational studies. We discuss the experimental challenges in forming three computationally predicted ternary phases (La₂SiP, La₅SiP₃, and La₂SiP₃). To understand the synthetic challenges, we performed molecular dynamics (MD) simulations using an accurate and efficient artificial neural network machine learning (ANN-ML) interatomic potential. We study the phase stability and formation kinetics of these ternary phases in relation to the reported and synthesized La₂SiP₄ phase. While the growth of the La₂SiP₄ phase can be reproduced by our MD simulation, our results indicate that the rapid formation of a Si-substituted LaP crystalline phase is a major barrier to the synthesis of the predicted La₂SiP, La₅SiP₃, and La₂SiP₃ ternary compounds, agreeing well with experimental observations. Our simulations also suggest that there is a narrow temperature window in which the La₂SiP₃ phase can be grown from the solid–liquid interface.

---

1. Introduction

Rapid advances in artificial intelligence/machine learning (AI/ML) algorithms, information infrastructures, and data sciences, as well as computer hardware/software offer promising opportunities to develop new, transformative strategies for materials design, discovery, and synthesis.^1–9 ML techniques have been demonstrated to be very successful in the rapid screening of vast composition–structure spaces to select promising candidates for first-principles calculations, greatly accelerating the pace of new materials prediction and discovery.^10–15

Currently, many computationally predicted new compounds have not been experimentally realized, largely due to the lack of prediction of their optimal synthetic conditions.¹⁶ Most computational predictions are performed at 0 K, while experimental syntheses are at finite temperatures. Knowledge about temperature effects on the relative thermodynamic stabilities among the predicted ground state structures and competing stable and metastable structures is essential for optimizing synthetic parameters. Moreover, as kinetics dominate during experimental synthesis, it is vital to understand the temperature dependence of stability and the phase formation and growth kinetics. Hence, it is highly desirable to combine an effective and timely computational–experimental feedback loop with AI/ML to enable efficient predictive synthesis.

The fundamental thermodynamic energy landscape of multi-element phase spaces and complex kinetics of crystal nucleation and growth make predictive synthesis an outstanding challenge. A full computational exploration of predictive synthesis is not realistic at present. However, advances in computational modeling/simulation and AI/ML provide opportunities to make relevant estimates for a narrow subset of parameters to guide experimental efforts. For example, understanding the potential decomposition pathways and products of various synthesis conditions is useful in eliminating a significant set of dead-end routes.^17,18 Simulations can then guide experimentalists through synthetic parameter selection and protocol design to enhance the efficacy of synthesis. Such efforts will greatly improve the synthetic predictability and rational design of emergent materials.

While molecular dynamics (MD) simulation is a powerful tool to study the thermodynamics and kinetics of real materials at the atomic level, the bottleneck of such simulations is the lack of accurate and efficient interatomic potentials for systems of interest, especially those with 3 or more chemical elements. Although ab initio MD can accurately describe interatomic interactions, it is computationally expensive. Larger-size and/or longer-time MD simulations can be carried out at low computational costs using empirical interatomic potentials in analytic functional forms with fitted parameters, but it is difficult to obtain accurate and transferable potentials for systems of 3+ elements. Recently, machine learning has emerged as a promising strategy to tackle challenges in constructing interatomic potentials that balance the accuracy and efficiency requirements of MD simulations.^19–40 In this paper, we use ternary La–Si–P compounds to demonstrate that reliable MD simulations can be realized using an artificial neural network (ANN) ML interatomic potential,⁴¹ yielding a greater understanding of synthetic challenges in complex chemical systems.

The La–Si–P phase space is relatively sparse, with only 5 experimentally reported phases: LaSi₂P₆, La₂SiP₄, and 3 polymorphs of LaSiP₃ (ref. 42–44). All five of these compounds are semiconductors featuring [SiP₄] tetrahedral units and homonuclear P bonding ranging from P–P dumbbells to one-dimensional (1D) P chains to two-dimensional (2D) disordered nets. Noncentrosymmetric LaSi₂P₆ and Pna2₁ LaSiP₃ are potential infrared (IR) nonlinear optical (NLO) materials,⁴⁵ and both LaSiP₃ layered phases (Pna2₁ and Aea2) exhibit low thermal conductivities, a key trait in candidate thermoelectric materials.⁴³ Thus, ternary La–Si–P compounds hold promise in energy-related applications and have prompted further computational investigations in this complex phase space. Recently, La-rich ternary phases La₂SiP and La₅SiP₃ were predicted by ML-guided first principles calculation studies to be low-energy metastable phases.¹⁴ ML-guided exploration of P-rich La–Si–P ternary compounds also revealed 16 stable and metastable compounds, including stable La₂SiP₃.¹⁵ In this paper, we address the synthetic difficulties in realizing these predicted materials. We perform MD simulations to study the phase stability and formation of predicted novel ternaries, La₂SiP, La₅SiP₃, and La₂SiP₃,^14,15 using an accurate and efficient ANN-ML interatomic potential developed by us.⁴¹ For comparison, we also study La₂SiP₄, which was recently synthesized by our experimental group.⁴² These simulations, coupled with experimental synthesis feedback, provide valuable information for understanding the inherent challenges in the synthesis of La–Si–P ternary compounds.

2. Experimental synthesis challenges

There are inherent challenges to the high-temperature solid state synthesis of La–Si–P phases. A direct reaction of the refractory elements, La and Si, with volatile and reactive red P promotes the formation of thermodynamic dead-ends, LaP and SiP, rather than the target phase.⁴³ However, employing an arc melted La–Si precursor in lieu of elements can aid in overcoming some synthetic barriers. Atomic mixing of the refractory elements allows them to be introduced into the reaction medium simultaneously with P in close spatial and temporal proximity. These LaSi + P reactions encourage prompt formation of ternary compounds upon heating and minimize binary admixtures. Another method to promote the reactivity of the elements is to employ a flux with a lower melting temperature, as demonstrated for the synthesis of LaSi₂P₆ in Sn and NaCl fluxes.^44,45 Yet, despite these synthetic tricks, the synthesis of ternary La–Si–P compounds at elevated temperatures (>1450 K) remains arduous due to reaction container limitations.

It should also be noted that while the chemical potential of P provides an advantage in the exploration and discovery of P-rich ternary compounds,¹⁴ P reactivity also brings about additional challenges. At sufficient temperatures, reaction contents may reach a melt which can further attack the reaction vessel, resulting in unwanted admixtures and potentially damaging the container. Additionally, most solid-state reactions are conducted in silica ampoules that begin to deform at T ∼ 1500 K, so synthesis at high temperatures is not trivial. Metal ampoules, such as Nb or Ta containers, are often used as an alternative to silica but these metals react with P at elevated temperatures to form Nb or Ta metal phosphides. A future avenue for the synthesis of La–Si–P phases at high temperature may be to employ carbon crucibles, but these containers still need to be sealed in a secondary silica jacket. Anticipating these synthetic hurdles, we use the La–Si–P phase space as a model for experimental–computational feedback to explore competition between phase stability and phase selection kinetics of complex ternary compounds.

The crystal structures of the low-energy ML-predicted La₂SiP (further indicated by the short notation 211), La₅SiP₃ (513), and La₂SiP₃ (213) compounds^14,15 and experimentally synthesized La₂SiP₄ (214) phase⁴² are shown in Fig. 1. In the La-rich 211 and 513 compounds, [P@La₆] octahedra connect via edges and share corners with [Si@La₆] trigonal prisms. However, in the P-rich compounds, neither of these motifs is observed; instead, [SiP₄] tetrahedra are present. In La₂SiP₃, these tetrahedra share corners to form 1D chains propagating down the [010] direction with [La@P₇] polyhedra filling space between these chains. In La₂SiP₄, additional P atoms are accommodated as P–P dumbbells connecting the [SiP₄] tetrahedra together in 1D chains along [001], and interstitial La atoms are coordinated by 9P atoms, [La@P₉], or Si + 8P atoms, [La@SiP₈]. The phonons of these ternary compounds have been studied by first principles calculations.^14,15 These compounds were shown to be dynamically stable without imaginary phonon modes being observed. The relative stability of the predicted 213 phase against the experimentally known 214, 113, and 126 phases as a function of temperature, including the contributions of vibrational entropy to the free energy, has been given in our previous paper.¹⁵

Fig. 1 The crystal structures of the 211, 513, 213, and 214 phases. In (a) 211 and (b) 513, [P@La₆] octahedra are shown in pink and [Si@La₆] trigonal prisms are shown in cyan. In (c) 213 and (d) 214, [SiP₄] tetrahedra are shown in blue. [La@P₇] polyhedra in (c) 213 are shown in yellow, while [La@P₉] or [La@SiP₈] polyhedra in (d) 214 are shown in green. Note that in 211 and 513 crystals, the [P@La₆] octahedra is the same building block of binary LaP. The half of [La@P₇] polyhedra in 213 has a pyramid structure (shown by pink bonds of P–P at the top of figure (c) as an example), which is just a half building block of the LaP crystal, while there is no octahedra-like motif in the 214 crystal. These LaP motifs in 211, 513, and 213 crystals are responsible for the rapid growth of the LaP phase during the crystal growth process, which will be demonstrated in the following MD simulations.

To corroborate ML-guided predictions, we previously attempted the synthesis of the predicted La-rich compounds, La₂SiP and La₅SiP₃.¹⁴ Despite our attempts to vary synthetic parameters (elemental vs. precursor starting materials, and reaction temperatures ranging from 773–1273 K), all experiments resulted in the formation of the stable binary LaP and unreacted Si. Because the synthesis of La-rich ternary phases was challenging, we instead focused our synthetic efforts on the P-rich compounds, hypothesizing that the chemical potential of P may help promote phase formation. We attempted several syntheses targeting the formation of La₂SiP₃ using arc-melted La–Si precursors, CsCl salt flux, and by varying synthesis temperature in the range of 1073–1423 K (Table S1). We were unsuccessful in forming the predicted La₂SiP₃ phase as revealed by powder X-ray diffraction (PXRD). However, these synthetic trials led to two significant observations: (i) LaP was present in all reactions, mostly as a stable major product, alluding that it is a barrier to the synthesis of ternary La–Si–P phases, and (ii) most reactions produced the known ternary phases Aea2 LaSiP₃ or La₂SiP₄ rather than La₂SiP₃. The reported melting point of Aea2 LaSiP₃ is ∼1308 K,⁴³ and this phase was present as a minor product in reactions conducted at 1073 K. However, by raising the reaction temperature to 1273 K while maintaining the same duration, unindexed diffraction peaks appeared that did not correspond to any known binaries and ternaries or predicted La–Si–P ternaries (Fig. S1(a) and S2(a)). Reactions conducted at 1273 K or 1423 K for longer durations yielded La₂Si₂O₇, indicating a partial reaction of the melt with the silica ampoule, and some La₂SiP₄ or unindexed peaks (Fig. S2(a)). Overall, this series of experiments revealed that melting La–Si–P phases may be a way to access new ternary phases in this system.

3. Melting temperature and its potential effects on synthesis

To better understand the results from the aforementioned synthetic attempts, we evaluate the thermal stabilities and melting temperatures of the predicted low-energy La₂SiP, La₅SiP₃, and La₂SiP₃ and experimentally synthesized La₂SiP₄ (ref. 42) phases by MD simulations using the developed ANN-ML La–Si–P interatomic potential. The performance of ANN-ML potential on crystal growth and melting can be found in Fig. S3. For more information about the ANN-ML interatomic potential used in the present MD simulations, we refer readers to ref. 41.

The MD simulations are performed using the LAMMPS package.⁴⁶ The initial simulation boxes are crystalline supercells containing 4608, 5184, 5184, and 6860 atoms for the 211, 513, 213, and 214 phases, respectively. An isothermal–isobaric (NPT) ensemble and a Nose–Hoover thermostat are used in the simulations.^47,48 The MD time step for the simulation is 2.5 fs. The four crystalline phases are first heated to 300 K for 500 ps, then continuously heated up to 3600 K or 2600 K, respectively, at a rate of 10¹² K s⁻¹. After thermalizing at high temperature, cooling simulations from the liquids are also performed with a continued cooling rate of 10¹² K s⁻¹. Fig. 2 shows the evolution of the instantaneous potential energy (E − 3k_BT) as a function of temperature upon heating and cooling for the 4 systems. We see a clear jump in the potential energy around 2300 K, 2200 K, 1900 K, and 1600 K, respectively, for the La₂SiP, La₅SiP₃, La₂SiP₃, and La₂SiP₄ phases, indicating crystal–liquid transitions are occuring. Snapshot structures of the simulation cell before and after melting are given as insets in Fig. 2 to illustrate changes in the atomic structures. We expect the temperatures where the potential energies rapidly increase should be well above the melting temperatures due to the finite system size and simulation time used. We then estimate the melting temperature T_m of each La–Si–P phase by checking the temperature where the potential energy (E − 3k_BT) starts to deviate from the linear temperature-dependence upon heating and cooling as shown in Fig. 2. The simulation results show the T_m estimated from heating and cooling in this way are very close with each other. The average T_m from heating and cooling simulations estimated in this manner also agree well with that obtained by solid/liquid coexistence method show below. Therefore, the temperature where the potential energy starts to deviate from the linear temperature-dependence upon both heating and cooling is a good indicator of the melting point. By comparison, the temperatures where the potential energies exhibit a rapid increase during the heating process are several hundred Kelvin (300–500 K) above the melting point of these compounds. Nevertheless, it is reasonable to believe that these compounds, if they can be synthesized, are stable at high temperatures of at least 1200 K.

Fig. 2 The potential energy E − 3k_BT as a function of temperature for (a) La₂SiP, (b) La₅SiP₃, (c) La₂SiP₃, and (d) La₂SiP₄ phases upon heating the crystalline phases (red) and cooling the liquid state (blue) in MD simulations using the ANN-ML interatomic potential from ref. 41. The insets are snapshot atomistic structures before and after melting during heating. The blue and red lines are linear fittings of E − 3k_BT in the low and high temperature regions, respectively. The two vertical dashed lines show the temperatures where the potential energy starts to deviate from its linear temperature-dependence upon heating and cooling, respectively.

We also calculate the melting temperatures of the 211, 513, 213, and 214 La–Si–P ternary phases using the solid–liquid phases coexistence method by MD simulations.⁴⁹ First, the supercells in the whole simulation boxes containing 3136, 3024, 4608, and 3584 atoms, respectively, for the 211, 513, 213, and 214 phases are performed with NPT ensemble simulations at temperatures below the melting points identified in Fig. 2. Then, the atoms in the right half of box are heated up to high temperature to ensure they are totally melted into a liquid state, while the atoms in the left half are fixed. After the atoms in the liquid half are cooled down to the initial temperature, the atoms in the whole box are simulated with NPT ensemble for a short time ∼10 ps to equilibrate the solid/liquid interface. Then, the MD simulations are switched to a NVE ensemble to search for the coexistence of crystalline and liquid phases by adjusting the total energy of the systems through scaling the velocities of the atoms. When the coexistence of crystalline and liquid phases is reached, the melting temperatures can be determined by averaging the kinetic energy of the system over a sufficient time in the MD simulation. The simulation results are shown in Fig. 3. For the 211 and 513 phases, we observe precipitation of a LaP binary phase, which is sandwiched between the parent (211 or 513) solid and liquid phases. The LaP phase precipitated from both the 211 and 513 liquids contains some Si atoms. The three phases (solid, liquid, and LaP) coexist during the MD simulation times of ∼0.9 and 0.6 ns respectively, and at similar temperatures, 1900 ± 30 K and 1930 ± 30 K respectively, for the 211 and 513 phases as shown in Fig. 3(a) and (b). The coexistence of the 213 and 214 phases with their corresponding liquids is also shown in Fig. 3(c) and (d), from which the melting temperatures are calculated to be 1340 ± 20 K and 1200 ± 20 K, respectively. More information about the crystal/liquid coexistence simulations is available in the SI (Fig. S4–S6).

Fig. 3 Crystal/liquid coexistence MD simulations for the (a) 211, (b) 513, (c) 213, and (d) 214 phases. The liquid regions are highlighted by the red dashed lines. While the 213 and 214 solid phases are in direct contact with their liquids, the LaP crystalline phase (highlighted by the blue lines) grows from the 211 and 513 liquids and separates the solid phases from the remaining liquid. Note that the precipitated LaP phases contain some Si atoms.

The MD simulation results reveal that 211 and 513 phases exhibit high melting temperatures above 1900 K. The high melting points of these compounds may make them difficult to synthesize under the current experimental conditions, given the potential reactivity of the La–Si–P melt and reaction container limitations. However, these phases could be accessed by other means, such as arc-melting. The predicted melting temperature, T_m, for LaP from our simulation⁴¹ is ∼2690 K, well beyond the temperature range accessible by standard solid state furnace reactions and the stability range of silica ampoules. It is possible LaP could be activated by heating above the melting point of the 211 phase, T_m ∼ 1900 K (Fig. 3). To probe this hypothesis, we first attempted to synthesize La₂SiP via arc melting (Table S2). A pre-reacted sample with nominal composition “La₂SiP” (containing LaP and a mixture of LaSi, La₅Si₄, and La₂Si₃) was arc-melted and the melted ingot was characterized by PXRD, revealing the presence of LaP, LaSi, and minor LaSi₂. Scanning electron microscopy (SEM) studies with energy-dispersive X-ray spectroscopy (EDS) elemental analysis also indicated the presence of minor amounts of a ternary La–Si–P phases in that sample. Poor homogenization of the pre-reacted sample and the inability to form 211 by arc melting may be due to the non-conductive nature of LaP. As such, we varied the ratio of LaSi : LaP by combining the LaSi precursor and LaP powder in a 0.8 : 0.2 ratio to ensure the conductive binary silicide was the majority phase. According to PXRD, this reaction produced binary phases only, with La₂Si₃ being the most prominent. Similar to the first arc-melted sample, SEM/EDS analysis identified only minor ternary La–Si–P phases present.

In contrast to the 211 and 513 phases, the melting temperature of the 213 and 214 phases are significantly lower. Compared to available experimental data for the melting temperatures of La–Si–P ternary compounds,⁴¹ the ANN-ML potential used in the present MD simulations systematically underestimates the melting temperature of LaSiP₃ and La₂SiP₄ by about 10%. We note that the melting point for Si calculated by GGA-DFT⁵⁰ is also about 10% lower than the experimental value. The melting point underestimation by our ANN-ML potential could be partially attributed to the training data, which are generated by GGA-DFT calculations. It is likely that the melting points for other ternary La–Si–P phases predicted by the ANN-ML potential are also about 10% smaller than their actual values. The predicted melting temperature for the 214 phase is about 1200 K from the simulation while the experimental value is ∼1330 K. According to this scaling, the experimental T_m of the newly predicted La₂SiP₃ phase should be about 1490 K. Employing arc-melted precursors and/or CsCl salt flux in experimental syntheses partially helped in overcoming synthetic barriers to form reported ternary La–Si–P compounds, albeit not the target La₂SiP₃ phase, as well as minor unidentified La–Si–P phases (Table S1). However, the primary hurdle in phase formation appears to be the presence of LaP as a major, stable product. The lower predicted melting temperatures for 213 and 214 phases indicate they may be more accessible by traditional synthetic methods. There is a ∼150 K temperature window between the predicted melting points of 214 and 213 (Fig. 3), which suggests it may be possible to transform one ternary phase to another and subsequently isolate 213. Based on the experimental T_m = 1330 K for La₂SiP₄ (ref. 41) determined by differential scanning calorimetry (DSC), we probed this hypothesis by heating a sample containing (LaP + La₂SiP₄) to 1373 K for 5 days and then quenching it in air from high temperatures. This process yielded LaP, minor amounts of La₂Si₂O₇, as well as a new crystalline phase (or phases) as indicated by the presence of unindexed diffraction peaks (Table S3 and Fig. S2(a)). Semi-quantitative elemental analysis by SEM/EDS revealed some possible ternary compositions (normalized to 1 La): LaSi_0.7(2)P_1.1(1) and LaSi_0.29(9)P_0.96(7) along with Si-doped LaP (Fig. S2(b) and (c)), indicating there may be additional unreported ternary La–Si–P phases. Attempts to improve the crystallinity of and isolate the unknown phase(s) are underway.

We attempted to resolve the identity of the unknown phases(s) through computational modeling and simulations. We performed adaptive genetic algorithm (AGA) searches^51,52 for low-energy structures with La₃SiP₃ (LaSi_0.33P) and La₄SiP₄ (LaSi_0.25P) compositions. We obtained a metastable La₄SiP₄ structure with a formation energy 67 meV per atom above the convex hull (Fig. S7). This structure features LaP slabs composed of [P@La₆] octahedra that are connected via Si–Si bonds between [P@La₄Si₂] distorted octahedra. The most intense reflections in the simulated XRD pattern of La₄SiP₄ match some of the unindexed diffraction peaks observed in experiments (Fig. S7(c)). We also constructed some LaSi_0.25P and LaSi_0.31P hypothetical structures by randomly inserting Si atoms into the interstitial sites of the LaP lattice, then performing structural relaxations to obtain optimized crystal structures. As shown in Fig. S8, the simulated PXRD patterns for these model structures are similar to that of pristine LaP with some additional diffraction peaks, much like the products observed in our experiments. These simulations suggest that the experimentally-observed ternary La–Si–P phase(s) may be structurally-related to LaP and the Si-substituted compounds shown in Fig. S7 and S8.

4. Crystallization at the solid–liquid interfaces

To gain more insight into the growth kinetics of the low-energy La–Si–P ternary phases, MD simulations were performed to investigate the growth of the crystalline phases at solid–liquid interfaces (SLI) for the four ternary La–Si–P compounds of interest. The SLI MD simulations are performed using rectangular boxes along certain crystallographic directions. The boxes contain 3136, 3024, 4608, and 3584 atoms, respectively, for the 211, 513, 213, and 214 phases – the same as those in the crystal/liquid coexistence MD simulations described above. To prepare the initial SLI configuration, we first run MD simulations for the crystalline La–Si–P structures at a given temperature T below the melting temperatures. Then, half of the atoms in the simulation box are kept fixed at a snapshot configuration, while another half of the atoms in the simulation box are heated (at a rate of 10¹² K s⁻¹) well above the melting point for 10 ps to ensure the atoms in this section are completely in a liquid state. Then, the atoms in the liquid are cooled down (at a rate of 10¹² K s⁻¹) to the given temperature T. Thus, an initial configuration of SLI with crystalline structure on the left and liquid on the right at a given T is prepared. Starting from the initial configuration of the SLI, MD simulations are performed at a constant temperature T and constant pressure (with pressure perpendicular to the interface) of 0 kbar to study the growth of the crystal from the SLI.

4.1. La₂SiP and La₅SiP₃ phases

MD simulations of crystallization at the SLI were performed at temperatures of T = 1520 K for the 211 phase and T = 1544 K for the 513 phase (which are about 0.80 of their melting temperature, T_m) using an NPT ensemble. The potential energies decrease with simulation time as shown in Fig. 4(a) and (b), suggesting continued crystallization from the interfaces. For the 211 phase, the potential energy decreases rapidly from 0–2 ns and around 6 ns, indicating rapid crystallization around these two intervals. For the 513 phase, the potential energy decreases rapidly in the first 2 ns and fluctuates near a constant value during the next 28 ns of simulation, suggesting very quick crystallization in this system. These crystallization behaviors are also confirmed by the snapshots of the atomic structures at different simulation times (Fig. 4(c) and (d)). However, the grown crystalline phases are mostly the binary LaP phase, as determined by comparison of the pair correlation functions of the grown structures to that LaP phase at the same temperature (see Fig. 4(e) and (g) for region (i) and (iii) respectively). The rapid growth of the LaP phase from the 211 and 513 interfaces are due to presence of [P@La₆] octahedra (Fig. 1), which is also the main motif of the LaP phase. It is interesting to see some of the 211 phase (the as-grown crystal in region (ii), where the g(r) is shown in Fig. 4(f)) forms between the grown LaP phase, presumably due to the richer La and Si compositions in the remaining liquid after the quick formation of LaP in the first 2 ns. The potential energy drop around the 6 ns is correlated with the formation of the 211 phase from the remaining liquid (highlighted by the red dashed line in Fig. 4(c)). These simulation results suggest that the 211 phase would be synthesizable starting with La- and Si-rich compositions. The persistence of some liquid (highlighted by the red dashed line in Fig. 4(d)) in the simulation of 513 phase would also be due to a richer La composition in the remaining liquid.

Fig. 4 MD simulation of crystal growth for (a) La₂SiP and (b) La₅SiP₃ phases at the crystal/liquid interface at T = 1520 K (211) and T = 1544 K (513). Figure (c) shows, in the 211/liquid MD simulation, 12 LaP crystal layers (yellow region) grown at two interfaces in the first 6 ns simulation time. In the as-grown LaP layers, about 40% P sites are substituted by Si atoms, and some Si atoms also occupy interstitial sites. The composition of the liquid (highlighted by red dashed lines) is La₂Si_1.84P_0.16. This liquid crystallizes into a part of the 211 structure beginning around 6 ns of simulation time. The blue region (ii) shows the whole 211 crystal structure, where the composition is La₂Si_1.3P_0.7, suggesting that the P atoms are substituted by Si atoms when 211 crystal grows. Similarly, figure the (d) shows the as-grown crystal in the 513/liquid simulation also adopts the LaP structure with some P atoms substituted by Si atoms. Comparison of the total pair-correlation functions between the perfect and as-grown crystals are shown in figures (e–g).

It is also interesting to note that the LaP phase grown from both the 211 and 513 interfaces are well-ordered with partial substitution of P by Si. In the 211 phase simulations, the grown LaP phase has 40% of P being substituted by Si, while only about 10% of P is substituted in the LaP phase grown from the 513 liquid (Fig. 4). This difference is attributed to the different Si content (25% vs. 11%) in the parent liquid. Substitution of Si in the LaP matrix was also confirmed experimentally by SEM/EDS on selected samples, ranging from ∼10% to ∼50% (Fig. S1 and S2). More details about the Si atomic positions in the grown LaP phase from the 211 SLI can be found in the SI (Fig. S9). These crystallization simulation studies provide further explanation for the observed experimental results wherein Si-substituted LaP was the major phase formed in the arc-melting experiments (Table S2). As arc-melting involves rapid heating and cooling rates at high temperatures, crystallization of thermodynamically stabilized products may be favored.

4.2. La₂SiP₃ and La₂SiP₄ phases

MD simulations of the crystallization at the La₂SiP₃ solid–liquid interface were performed at a temperature T = 1187 K (which is about 0.89 of the melting T) using an NPT ensemble. The Fig. 5(a) shows the potential energy per atom as a function of simulation time. The potential energy exhibits a quick drop in several steps (around 1, 5, and 15 ns), suggesting that crystallization taking place in steps. It is interesting to see both 213 and LaP grown after 15 ns of simulation (indicated by (ii) and (i) in the plots). The crystallization snapshots at different simulation times are given in Fig. 5(c). We verify that the grown crystalline phases are LaP and 213 by comparing the pair correlation functions of the as-grown structure at the same temperature, as shown in the Fig. 5(e) and (f). We note that, when the simulations are performed at a temperature of 1225 K, only the growth of the 213 phase is observed as shown in Fig. S10(b). However, simulations at other temperatures (1200 K, 1250 K, and 1275 K) show that both LaP and 213 phases grow at the 213/liquid interface (see Fig. S10(a), (c), and (d)). These results can explain why LaP phase is constantly obtained in the experimental synthesis attempts discussed above. There may be a narrow temperature window of about 75 K where the 213 can grow above the melting point of 214 phase. However, realizing such a narrow temperature window in experimental synthesis would be a challenge. In comparison, the growth of the 214 phase at T = 1050 K (which is about 0.88 of melting T_m) takes about 30 ns of the MD simulations to complete over the simulation cell, as shown in Fig. 5(b) and (d). The pair correlation function of the as-grown structure confirms that it is indeed the 214 phase (see Fig. 5(g)). The simulated growth of the 214 phase, without interference by the LaP phase, also agrees with our experimental results.

Fig. 5 MD simulation of crystal growth of (a) La₂SiP₃ and (b) La₂SiP₄ phases at the crystal/liquid interface at T = 1187 K (213) and T = 1050 K (214). Figure (c) shows both LaP and 213 grow at the 213/liquid interface (indicated by (i) and (ii) in the plots). In the as-grown LaP structure, some P atoms are substituted by Si atoms, similar to the 211/liquid and 513/liquid simulations. Figure (d) shows 214 growing at the 214/liquid interface without interference by the LaP phase. Comparison of the total pair-correlation functions between the perfect and as-grown crystals are shown in figure (e–g).

5. Summary

In this paper, we explore synthetic challenges of La–Si–P ternary compounds through feedback between experiments and computational simulations. Varying experimental attempts, including arc-melting and high-temperature quenching, failed to produce the predicted La₂SiP, La₅SiP₃, and La₂SiP₃ compounds. Instead, the results of these syntheses showed that LaP was present in all reactions, mostly as a stable major product. This revealed that the LaP phase is indeed a significant barrier to the synthesis of ternary La–Si–P phases.

Using an ANN-ML interatomic potential for La–Si–P ternaries, we performed MD simulations to investigate the thermal stabilities and melting points of the La₂SiP, La₅SiP₃, and La₂SiP₃ compounds, as well as known binary LaP and ternary La₂SiP₄ phases. Our simulations show that the binary LaP phase is very stable with a melting temperature of 2690 K. All four ternary La–Si–P phases are also thermodynamically stable upon heating above 1200 K. The melting point of the La₂SiP₄ phase is the lowest among the studied phases. From the crystal/liquid coexistence MD simulations, the melting temperature of the La₂SiP₄ phase is calculated to be about 1200 K, which is about 10% lower than the experimental value. The agreement of the melting temperature of this ternary compound between the simulation and experiment is reasonably good, considering that the melting temperature is not explicitly included in the interatomic potential training.⁴¹ As the ANN-ML potential MD simulations revealed a ∼150 K temperature window between the predicted melting points of the 214 and 213 phases, we aimed to synthetically transform the 214 phase to the 213 phase by annealing a sample containing La₂SiP₄ to 1373 K (above the melting point of 214). The process yielded LaP, minor amounts of La₂Si₂O₇, as well as a new crystalline phase (or phases) as indicated by the presence of unindexed diffraction peaks.

The accuracy and efficiency of the ANN-ML interatomic potential for La–Si–P ternaries enabled us to perform MD simulations to study the growth kinetics of the crystalline phases at the solid–liquid interfaces for La₂SiP, La₅SiP₃, La₂SiP₃ and La₂SiP₄. In the 211/513 crystal growth simulations, MD studies correctly identify that the formation of Si-substituted LaP is a major kinetic barrier to ternary formation, which was consistent with experimental observations. Our simulations also suggest that the 211 phase could be formed under La- and Si-rich conditions. In the 214/213 crystal growth simulations, while the 213 liquid was easily transformed into Si-substituted LaP, growth of 214 was obtained below the melting temperature. This result also explains that, while the 214 phase can be made experimentally, strong competition from the LaP phase is a major barrier to synthesizing the 213 phase. Although the simulations indicate that there is a narrow temperature window below melting to facilitate the growth of 213 phase, realizing such a narrow temperature window in experimental synthesis would be difficult.

Although DFT calculations suggest that the 211, 513 and 213 phases have favorable formation energy, the convex hull was calculated at 0 K and the stability of these phases might be altered significantly at high temperatures. The experimental synthetic approach was semi-successful for the studied La–Si–P system, given the challenges related to highly reactive La in the melt and high P vapor pressure which inhibited certain synthetic endeavors. The proposed methodology should be more applicable for systems requiring less extreme synthetic temperatures, i.e. T_m < 1000 K, and/or for systems lacking such volatile components, such as borides and silicides. The melting temperatures of the predicted 211, 513 and 213 phases are all well above 1300 K, especially for the La-rich compositions, so experimental challenges are expected. Our MD simulations using an ANN-ML interatomic potential provided valuable information and insights to understand the observed experimental synthesis challenges and also suggested some other potential avenues for synthesis in future investigations. A comprehensive phase diagram establishing the relative stability of competing solid and liquid phases under varying conditions (such as chemical compositions and temperatures, pressures) would be critical for determining the synthesis pathways for complex materials. Unfortunately, to the best of our knowledge, no such ternary phase diagram for the La–Si–P system is currently available. Constructing a comprehensive and reliable ternary phase diagram computationally remains challenging and intensive, primarily due to the need for precise free energy calculations for both liquid and solid phases across vast chemical and thermal ranges. This area warrants further study, given that phase diagrams are foundational for guiding the synthesis of complex materials.

Conflicts of interest

The authors declare no competing interests.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Supplementary information: the movies of MD simulations and the structure file of new metastable La₄SiP₄ crystal found by AGA searching. See DOI: https://doi.org/10.1039/d5ta05069c.

Acknowledgements

Work at Ames National Laboratory was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Science and Engineering Division, including a grant of computer time at the National Energy Research Scientific Computing Center (NERSC) in Berkeley. Ames National Laboratory is operated for the U.S. DOE by Iowa State University under contract # DE-AC02-07CH11358. We thank Dr Yaroslav Mudryk (Ames National Laboratory) for access to the arc-melting setup. L. Tang acknowledges the support by the National Natural Science Foundation of China (Grant No. 11304279).

References

1. J. E. Gubernatis and T. Lookman, Machine learning in materials design and discovery: examples from the present and suggestions for the future, Phys. Rev. Mater., 2018, 2, 120301.
2. R. Vasudevan, G. Pilania and P. V. Balachandran, Machine learning for materials design and discovery, J. Appl. Phys., 2021, 129, 070401.
3. A. G. Kusne, T. Gao, A. Mehta, L. Q. Ke, M. C. Nguyen, K. M. Ho, V. Antropov, C. Z. Wang, M. J. Kramer, C. Long and I. Takeuchi, On-the-fly machine-learning for high-throughput experiments: search for rare-earth-free permanent magnets, Sci. Rep., 2014, 4, 6367.
4. A. Kabiraj, M. Kumar and S. Mahapatra, High-throughput discovery of high curie point two-dimensional ferromagnetic materials, npj Comput. Mater., 2020, 6, 35.
5. J. Cai, X. Chu, K. Xu, H. Li and J. Wei, Machine learning-driven new material discovery, Nanoscale Adv., 2020, 2, 3115–3130.
6. G. Katsikas, S. Charalampos and K. Joseph, Machine learning in magnetic materials, Phys. Status Solidi B, 2021, 258, 2000600.
7. T. D. Rhone, W. Chen, S. Desai, S. B. Torrisi, D. T. Larson, A. Yacoby and E. Kaxiras, Data-driven studies of magnetic two-dimensional materials, Sci. Rep., 2020, 10, 15795.
8. G. A. Landrum and H. Genin, Application of machine-learning methods to solid-state chemistry: Ferromagnetism in transition metal alloys, J. Solid State Chem., 2003, 176, 587–593.
9. I. Miyazato, T. Tanaka and K. Takahashi, Accelerating the discovery of hidden two-dimensional magnets using machine learning and first principle calculations, J. Phys.: Condens. Matter, 2018, 30, 06LT01.
10. T. Xie and J. C. Grossman, Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties, Phys. Rev. Lett., 2018, 120, 145301.
11. C. W. Park and C. Wolverton, Developing an improved crystal graph convolutional neural network framework for accelerated materials discovery, Phys. Rev. Mater., 2020, 4, 063801.
12. W. Xia, M. Sakurai, B. Balasubramanian, T. Liao, R. Wang, C. Zhang, H. Sun, K.-M. Ho, J. R. Chelikowsky, D. J. Sellmyer and C.-Z. Wang, Accelerating novel magnetic materials discovery using a machine learning guided adaptive feedback, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2204485119.
13. R. H. Wang, W. Y. Xia, T. J. Slade, X. Y. Fan, H. F. Dong, K. M. Ho, P. C. Canfield and C. Z. Wang, ML-guided discovery of ternary compounds involving La and immiscible Co and Pb elements, npj Comput. Mater., 2022, 8, 258.
14. H. J. Sun, C. Zhang, W. Y. Xia, L. Tang, G. Akopov, R. H. Wang, K. M. Ho, K. Kovnir and C. Z. Wang, Machine learning guided discovery of ternary compounds containing La, P and group IV elements, Inorg. Chem., 2022, 61, 16699–16706.
15. W. Y Xia, L. Tang, H. J. Sun, C. Zhang, K. M. Ho, G. Viswanathan, K. Kovnir and C. Z. Wang, Accelerating materials discovery using integrated deep machine learning approaches, J. Mater. Chem. A, 2023, 11, 25973–25982.
16. K. Kovnir, Predictive Synthesis, Chem. Mater., 2021, 33, 4835–4841.
17. J. R. Neilson, M. J. McDermott and K. A. Persson, Modernist materials synthesis: Finding thermodynamic shortcuts with hyperdimensional chemistry, J. Mater. Res., 2023, 38, 2885–2893.
18. A. J. Martinolich and J. R. Neilson, Toward Reaction-by-Design: Achieving Kinetic Control of Solid State Chemistry with Metathesis, Chem. Mater., 2017, 29, 479–489.
19. J. Behler and M. Parrinello, Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces, Phys. Rev. Lett., 2007, 98, 146401.
20. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons, Phys. Rev. Lett., 2010, 104, 136403.
21. S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt and K.-R. Müller, Machine learning of accurate energy-conserving molecular force fields, Sci. Adv., 2017, 3, e1603015.
22. A. V. Shapeev, Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials, Multiscale Model. Simul., 2016, 14, 1153–1173.
23. 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., 2020, 2, 025002.
24. Z. Fan, W. Chen, V. Vierimaa and A. Harju, Efficient molecular dynamics simulations with many-body potentials on graphics processing units, Comput. Phys. Commun., 2017, 218, 10–16.
25. Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen and T. Ala-Nissila, Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport, Phys. Rev. B, 2021, 104, 104309.
26. Y. Zhang, J. Xia and B. Jiang, Physically Motivated Recursively Embedded Atom Neural Networks: Incorporating Local Completeness and Nonlocality, Phys. Rev. Lett., 2021, 127, 156002.
27. Y. Zhang, J. Xia and B. Jiang, REANN: A PyTorch-based end-to-end multi-functional deep neural network package for molecular, reactive, and periodic systems, J. Chem. Phys., 2022, 156, 114801.
28. J. Han, L. Zhang, R. Car and E. Weinan, Deep potential: a general representation of a many-body potential energy surface, Commun. Comput. Phys., 2018, 23, 629–639.
29. H. Wang, L. Zhang, J. Han and E. Weinan, DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics, Comput. Phys. Commun., 2018, 228, 178–184.
30. L. Zhang, J. Han, H. Wang, R. Car and E. Weinan, Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics, Phys. Rev. Lett., 2018, 120, 143001.
31. T. Wen, L. Zhang, H. Wang, E. Weinan and D. J. Srolovitz, Deep potentials for materials science, Mater. Futures, 2022, 1, 022601.
32. L. Zhang, J. Han, H. Wang, W. Saidi and R. Car, End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems, Advances in Neural Information Processing Systems, 2018, vol. 31.
33. L. Tang, K. M. Ho and C. Z. Wang, Molecular dynamics simulation of metallic Al-Ce liquids using a neural network machine learning interatomic potential, J. Chem. Phys., 2021, 155, 194503.
34. L. Tang, Z. J. Yang, T. Q. Wen, K. M. Ho, M. J. Kramer and C. Z. Wang, Development of interatomic potential for Al–Tb alloys using a deep neural network learning method, Phys. Chem. Chem. Phys., 2020, 22, 18467–18479.
35. L. Tang, Z. J. Yang, T. Q. Wen, K. M. Ho, M. J. Kramer and C. Z. Wang, Short- and medium-range orders in Al₉₀Tb₁₀ glass and their relation to the structures of competing crystalline phases, Acta Mater., 2021, 204, 116513.
36. J. Wang, H. Shen, R. Yang, K. Xie, C. Zhang, L. Chen, K.-M. Ho, C.-Z. Wang and S. Wang, A deep learning interatomic potential developed for atomistic simulation of carbon materials, Carbon, 2022, 186, 1–8.
37. J. Wu, Y. Zhang, L. Zhang and S. Liu, Deep learning of accurate force field of ferroelectric HfO₂, Phys. Rev. B, 2021, 103, 024108.
38. C. Zhang, Y. Sun, H.-D. Wang, F. Zhang, T.-Q. Wen, K.-M. Ho and C.-Z. Wang, Crystallization of the P₃Sn₄ Phase upon Cooling P₂Sn₅ Liquid by Molecular Dynamics Simulation Using a Machine Learning Interatomic Potential, J. Phys. Chem. C, 2021, 125, 3127–3133.
39. C. Zhang, L. Tang, Y. Sun, K.-M. Ho, R. M. Wentzcovitch and C.-Z. Wang, Deep machine learning potential for atomistic simulation of Fe-Si-O systems under Earth's outer core conditions, Phys. Rev. Mater., 2022, 6, 063802.
40. T. Q. Wen, C. Z. Wang, M. J. Kramer, Y. Sun, B. L. Ye, H. D. Wang, X. Y. Liu, C. Zhang, F. Zhang, K.-M. Ho and N. Wang, Development of a deep machine learning interatomic potential for metalloid-containing Pd-Si compounds, Phys. Rev. B, 2019, 100, 174101.
41. L. Tang, W. Xia, G. Viswanathan, E. Soto, K. Kovnir and C.-Z. Wang, Developing a Neural Network Machine Learning Interatomic Potential for Molecular Dynamics Simulations of La-Si-P Systems, J. Chem. Phys., 2025, 163, 084109.
42. G. Akopov, G. Viswanathan and K. Kovnir, Synthesis, Crystal and Electronic Structure of La₂SiP₄, Z. Anorg. Allg. Chem., 2021, 647, 91.
43. G. Akopov, J. Mark, G. Viswanathan, S. J. Lee, B. C. McBride, J. Won, F. A. Perras, A. L. Paterson, B. Yuan, S. Sen, A. N. Adeyemi, F. Zhang, C.-Z. Wang, K.-M. Ho, G. J. Miller and K. Kovnir, Third Time's the Charm: Intricate Non-centrosymmetric Polymorphism in LnSiP₃ (Ln = La and Ce) Induced by Distortions of Phosphorus Square Layers, Dalton Trans., 2021, 50, 6463–6476.
44. P. Kaiser and W. Jeitschko, The Rare Earth Silicon Phosphides LnSi₂P₆ (Ln = La, Ce, Pr, and Nd), J. Solid State Chem., 1996, 124, 346–352.
45. Y. S. Sun, J. D. Chen, S. D. Yang, B. X. Li, G. L. Chai, C. S. Lin, M. Luo and N. Ye, LaSiP₃ and LaSi₂P₆: Two Excellent Rare-Earth Pnictides with Strong SHG Responses as Mid- and Far-Infrared Nonlinear Optical Crystals, Adv. Opt. Mater., 2021, 9, 2002176.
46. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1.
47. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
48. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
49. J. R. Morris, C. Z. Wang, K. M. Ho and C. T. Chan, Melting line of aluminum from simulations of coexisting phases, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 3109.
50. D. Alfe and M. J. Gillan, Exchange-correlation energy and the phase diagram of Si, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 205212.
51. X. Zhao, M. C. Nguyen, W. Y. Zhang, C. Z. Wang, M. J. Kramer, D. J. Sellmyer, X. Z. Li, F. Zhang, L. Q. Ke, V. P. Antropov and K. M. Ho, Exploring the structural complexity of intermetallic compounds by an adaptive genetic algorithm, Phys. Rev. Lett., 2014, 112, 045502.
52. S. Q. Wu, M. Ji, C. Z. Wang, M. C. Nguyen, X. Zhao, K. Umemoto, R. M. Wentzcovitch and K. M. Ho, An adaptive genetic algorithm for crystal structure prediction, J. Phys.: Condens. Matter, 2014, 26, 035402.

---