Weiyi
Xia
a,
Ling
Tang
b,
Huaijun
Sun
c,
Chao
Zhang
d,
Kai-Ming
Ho
e,
Gayatri
Viswanathan
af,
Kirill
Kovnir
af and
Cai-Zhuang
Wang
*ae
aAmes National Laboratory, U.S. Department of Energy, Ames, IA 50011, USA
bDepartment of Applied Physics, College of Science, Zhejiang University of Technology, Hangzhou, 310023, China
cJiyang College of Zhejiang Agriculture and Forestry University, Zhuji 311800, China
dDepartment of Physics, Yantai University, Yantai 264005, China
eDepartment of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA
fDepartment of Chemistry, Iowa State University, Ames, IA 50011, USA
First published on 29th September 2023
We present an integrated deep machine learning (ML) approach that combines crystal graph convolutional neural networks (CGCNN) for predicting formation energies and artificial neural networks (ANN) for constructing interatomic potentials. Using the La–Si–P ternary system as a proof-of-concept, we achieve a remarkable speed-up of at least 100 times compared to high-throughput first-principles calculations. The ML approach successfully identifies known compounds and uncovers 16 new P-rich compounds with formation energies within 100 meV per atom above the convex hull, including a stable La2SiP3 phase. We also employ the developed ML interatomic potential in a genetic algorithm for efficient structure search, leading to the discovery of more metastable compounds. Moreover, substitution of La atoms with Ba reveals a new stable quaternary compound, BaLaSiP3. Our generic and robust approach holds great promise for accelerating materials discovery in various compounds, enabling more efficient exploration of complex chemical spaces and enhancing the prediction of material properties.
Rapid advances in AI/ML algorithms, information infrastructures and data science, as well as computer hardware/software, offer great opportunities to develop new transformative strategies for materials design and discovery.15–22 Recently, machine learning (ML) techniques have been employed to assist in rapidly screening vast composition-structure spaces to select promising candidates for first-principles calculations. Successful application of such ML-guided approaches in accelerating materials discovery has been demonstrated.23–27 However, the efficiency of these ML-guided approaches relies heavily on the accuracy of the ML model predictions. ML models with poor accuracy can either substantially overestimate the number of candidate structures, thus putting an unnecessary burden on subsequent first-principles calculations, or badly miss promising target structures. Thus, it is highly desirable to have a robust approach to provide accurate ML screening, selecting only a minimal number of promising candidate structures for first-principles calculations, and the success rate of the ML-guided predictions is boosted to dramatically accelerate the pace of materials discovery.
In this paper, we show that a crystal graph convolutional neural network (CGCNN) ML model, trained by first-principles calculation data specifically focus on the material system of interest, gives much more accurate prediction of ternary compound formation energies to efficiently select a small fraction of candidate structures (<1.5%) from a large hypothetical structure pool. Subsequent structure relaxation for the structures selected by CGCNN using interatomic potentials trained by an artificial neural network ML (ANN-ML) method, further reduces the candidate structures to a very small number for final refinement by fist-principles calculations, thus substantially boosting the efficiency of the ML-guided approach. Moreover, the developed ANN-ML interatomic potentials can be used to perform efficient and reliable structure searches based on GA. Such an integrated ML-guided approach is illustrated in Fig. 1.
Fig. 1 Flowchart of the proposed integrated ML-guided approach for accelerating materials discovery. |
We use La–Si–P ternary system as a test bed to demonstrate the accuracy and efficiency of the proposed integrated ML approach for the discovery of novel ternary compounds. Many reported silicon phosphides of transition and/or rare-earth metals crystallize in noncentrosymmetric (NCS) crystal structures. The absence of inversion symmetry and significant hybridization of d-, f-, and p-orbitals promotes a plethora of emergent properties such as unconventional superconductivity, topologically non-trivial quantum properties over large energy windows, and quasiparticle behavior.28–31
Two experimentally-observed ternary compounds (La(SiP3)2 and LaSiP3) have been reported,32,33 and another compound (La2SiP4) has also been recently synthesized by experiment34 along with a high-temperature disordered polymorph of LaSiP3. The discovery of new ternary phases in this system is both intriguing and challenging due to the complex interatomic interactions involved. Our proposed integrated ML approach not only correctly captures the three known ternary phases, but also efficiently predicts a novel energetically stable La2SiP3 phase and 15 other low-energy metastable P-rich ternary phases (P composition ≥ 50%) with formation energies within 0.1 eV per atom above the ternary convex hull. Among the 15 low-energy metastable P-rich ternary phases (P composition ≥ 50%), 5 of them are NCS crystals.
In CGCNN, a crystal structure is represented by a crystal graph that encodes both atomic information and bonding interactions between atoms. A convolutional neural network is added on top of the crystal graph to construct the proper descriptors, which are optimized for predicting target properties.23 In this way, composition–structure–property relationships can be efficiently learned and predicted by CGCNN. The training data in CGCNN are primarily generated by first-principles calculations, which enables a sufficient volume of data for the supervision training. The first CGCNN motel for formation energy predictions of compounds was proposed and trained by Xie and Grossman23 using the structures and formation energies of 28046 structures from first-principles calculations from the Materials Project (MP) database.35 The training structures in this model cover a wide range of chemical elements and composition ratios, and we refer to this model as a general CGCNN (g-CGCNN) model. While the g-CGCNN model has been shown to be very useful for accelerating the discovery of novel complex ternary compounds,25–27 the accuracy and efficiency of the CGCNN model can be significantly improved further if the model is retrained using the data specifically targeting the materials being studied.25,36 For La–Si–P ternary system, we retrained the CGCNN using 228284 structures and formation energies from the dataset that was prepared to train an ANN-ML interatomic potential for the La–Si–P system (see below and ESI†). We refer to this later CGCNN model as a specific CGCMM (s-CGCNN) model. In the training of both g-CGCNN and s-CGCNN, we use the exact same descriptors, size of neural network, and all other hyperparameters. The only difference is the training data. While the training data for g-CGCNN cover many chemical elements in the periodic table, the training data for the s-CGCNN model focuses on the structures containing only the three relevant elements La, Si, and P. To compare the prediction accuracy of these two CGCNN models for La–Si–P ternary compounds, we applied the models to a set of 806 La–Si–P ternary compounds studied in ref. 27. The first-principles calculation results for formation energies of these ternary compounds were studied in ref. 27, but the first-principles results were not included in the training dataset for either CGCNN model. The comparison of the formation energy from the CGCNN models and first-principles calculation results are shown in Fig. 2, in which mean absolute error (MAE) are provided. We can see that the g-CGCNN model substantially underestimates the formation energies of these compounds (Fig. 2a), while the agreement with the first-principles results is much better by the s-CGCNN predictions (Fig. 2b).
We then applied the two CGCNN models to perform a fast screening of the formation energies of ternary La–Si–P compounds in a hypothetical structure pool, which was generated based on known ternary structure lattices. We collected 28472 known ternary structures from the MP database and replaced the three elements with La, Si, and P. For each ternary structure from the MP database, five hypothetical lattices were generated by uniformly scaling the bond lengths of the structure such that the volume of the unit cell varied from the original by a factor of 0.96 to 1.04 in increments of 0.02. The use of the scaling factor for the crystal unit cell volume helps the CGCNN model differentiate the energetic stability of the same structure with different bond lengths. There are also six ways to shuffle the three elements in a given template ternary structure. Therefore, by multiplying reported ternary structures by 6 compositional combinations and 5 volume steps, 854070 hypothetical ternary La–Si–P structures were generated.
The distributions of the formation energies (Ef) from the g-CGCNN and s-CGCNN predictions for this set of structures are shown in Fig. 3, where Ef is defined with respect to the elemental ground-state bulk phases of the constituent elements La, Si, and P, i.e., , with E(LamSinPp) as the total energy of the LamSinPp compound, E(La), E(Si) and E(P) as the per-atom energy of the ground state of La, Si and P crystals, respectively. Negative formation energies indicates that the formation of the ternary compounds is energetically favorable with respect to the three elemental phases.
Fig. 3 The formation energies (Ef) of 854070 hypothetical ternary La–Si–P compounds as predicted by (a) the g-CGCNN model and (b) the s-CGCNN model. |
From Fig. 3a, we can see that more than 80% of the ternary La–Si–P structures from the hypothetical structure pool are predicted to have negative formation energies by the g-CGCNN model. In particular, the three known ternary compounds lie in the 42nd percentile low-energy region, meaning that about the 42% of the structures from the hypothetical structures must be selected as candidate structures in order to capture these known phases. Thus the g-CGCNN model is not efficient in accelerating materials discovery for this system. In comparison, the prediction from the s-CGCNN model, as shown in Fig. 3b, shows dramatic improvement in the efficiency of candidate structure selection. The fraction of the structures that have negative formation energy is much less (33%), and the three known ternary phases lie in the 5.6th percentile low-energy region, meaning that only 5.6% of the structures must be selected as candidates to capture the three experimentally-observed phases.
After removing very similar structures, our s-CGCNN screening revealed 12383 La–Si–P structures with negative formation energies that could be selected for further evaluations. While using first-principles calculation to optimize the unit cell shape and the atomic positions of the 12383 structures is manageable, we found that using an interatomic potential trained by an artificial neural network ML (ANN-ML) method can further boost the efficiency of the ML-guided approach. More details of the ANN-ML interatomic potential for the La–Si–P system are given in the ESI.† Structural optimization with the ANN-ML potential using the LAMMPs code is 105–106 times faster than that using regular DFT potentials,37 and enables the optimization of the 12383 structures to be done much quickly. The accuracy and efficiency of the ANN-ML interatomic potentials also enable efficient exploration of the temperature effects on the thermodynamic stability and the phase formation kinetics of the predicted phases, as will be discussed below.
The artificial neural network (ANN) machine learning (ML) interatomic potential for La–Si–P system is developed by the deep learning software package DeePMD-kit.44 The training data including the energies and forces are generated by ab initio MD (AIMD) simulations and ab initio calculations using the VASP package39,40 with the projector-augmented-wave (PAW) method.41 The exchange and correlation energy functional with non-spin-polarized generalized gradient approximation (GGA) in the Perdew–Burke–Ernzerhof (PBE) form41 is used. The energy cutoff for the plane wave basis set is 270 eV. The time step is 3 fs and NVT ensemble with Nose–Hoover thermostat45,46 are employed in AIMD simulations. The Brillouin zone is sampled by the gamma point. The training data set consists of 228284 structures and their formation energies, including liquid snapshot structures of La20Si20P60, La10Si10P80, La10Si45P45, La10Si80P10, La5Si90P5, La5Si5P90, La, Si, P, and the distorted crystalline structures of LaSi, LaP, SiP, LaSiP3, etc. from Materials Project database.35 The local structure configuration information within a cutoff radial of 7.0 Å is sampled to train the ANN-ML model with four hidden layers and 120 nodes per layer. The tanh function is used as activation function in neural network model.
The binding energies of 12383 structures and the known ternary and binary compounds obtained from the ANN-ML potential optimization enabled us to calculate the formation energy above convex hull (denoted as Ehull) for these structures, which gives a better quantitative description of the thermodynamic stability in comparison with using Ef. Ehull of any given phase on the ternary convex hull can be calculated by comparing the phase's formation energy with respect to three nearby known phases on the convex hull. These three known phases can be ternaries, binaries, or elemental crystalline phases, and the chemical compositions of these three phases correspond to the vertexes of a triangle (called a Gibbs triangle) that encloses the composition of the phase of interest for which Ehull is calculated. Therefore, Ehull determines the thermodynamic stability of the given phase against decomposition into the nearby three known phases.
The distribution of formation energies above convex hull (Ehull) of the 12383 structures predicted by the ANN-ML potentials is shown in Fig. 4. We can see that only a very small fraction of structures has Ehull below or close to the convex hull from the ANN-ML potential prediction. The known LaSiP3, La(SiP3)2, and La2SiP4 phases are predicted to have Ehull values of −0.036 eV per atom, 0.007 eV per atom, and 0.041 eV per atom respectively – all very close to the convex hull. Applying an Ehull cutoff value of 0.05 eV per atom to select candidate structures with composition ratio P ≥ 50%, we obtain 29 candidates including the three known ternary phases. These results indicate the ANN-ML potentials has sufficient accuracy to efficiently select promising candidate structures for the discovery of low-energy ternary compounds. To search for more low-energy metastable structures and to take into the consideration that the ANN-ML would have some errors in the Ehull predictions, we selected structures with Ehull ≤ 0.3 eV per atom and composition ratio P ≥ 50% predicted by the ANN-ML potential for further evaluation by first-principles calculations. These were 353 structures to be studied by first-principles calculations which could easily be performed within a few days using commonly available cluster computers, since most of the structures after the ANN-ML potential relaxation are already very close to the final first-principles calculation results. Therefore, using ANN-ML potentials to relax the structure can substantially reduce the number of candidate structures for first-principles calculations and provides another significant speed up to materials discovery on the top of s-CGCNN screening.
Fig. 4 The distribution of formation energies above convex hull (Ehull) of the 12383 structures predicted by the ANN-ML potential. |
The first-principles calculations were performed using the projector augmented wave (PAW) method38 within DFT as implemented in VASP code.39,40 The exchange and correlation energy is treated by the generalized gradient approximation (GGA) and parameterized by the Perdew–Burke–Ernzerhof formula (PBE).41 A plane-wave basis set with a kinetic energy cutoff of 520 eV was used to expand the electronic wave functions, ensuring high accuracy in our DFT calculations for the candidate structures identified by the ML process. The convergence criterion for the total energy was set to 10−5 eV. Monkhorst–Pack's sampling scheme42 was adopted for the Brillouin zone sampling with a k-point grid of 2π × 0.033 Å−1, and the unit cell lattice vectors (both the unit cell shape and size) were fully relaxed together with the atomic coordinates until the forces on each atom was less than 0.01 eV Å−1.
We also performed phonon calculations to investigate the dynamic stabilities of the two predicted La2SiP3 phases. The Phonopy package was used with the finite displacement method. A supercell of 2 × 2 × 2 is created, and the forces were calculated with the supercell with displacement along different directions. These DFT calculations adopted the same parameters used previously in the structural optimization as well as the formation energy calculations. Force constants were then calculated from the set of forces. Dynamic matrices were built from the force constants and thus phonon frequencies and eigenvectors were obtained with the specified q points. The phonon dispersion along with the phonon density of states are shown in Fig. 5c and d. No imaginary vibrational frequencies were found, indicating that these two phases are both dynamically stable. La atoms dominate the contribution of low frequency vibrational modes up to 4 THz in both structures. Both Si and P atoms contribute mostly to high frequency modes from 4 to 13 THz. The calculated band structures at the DFT-PBE level for these two structures are shown in Fig. S5a and b respectively in the ESI.† They both show metallic character but with small density of states around Fermi level. The bands across the Fermi level indicate a layered-like character formed. By investigating the density of states around the Fermi level, we found that the contribution mostly comes from La atoms, which is consistent with the arrangement of La atoms as layers in the two structures.
To investigate the effect of temperature on thermodynamics stability of the two newly discovered La2SiP3 compounds with respect to the three reported ternary compounds, we evaluated the Gibbs free energy as a function of temperature for the two La2SiP3 phases, the ordered LaSiP3 polymorph, LaSi2P6, and La2SiP4. Without pressure and at a given temperature T and volume V, the Gibbs free energy G of a crystalline compound is given by:
The results of the Gibbs free energy as the function of temperature for the five ternary La–Si–P phases from our calculations are shown in Fig. S6 in the ESI.† To characterize the thermodynamic stability of these ternary structures with different compositions, we calculated the formation Gibbs energies Ghull of each ternary phase with respect to the nearby competing binary and ternary phases in the ternary convex hull at each temperature. To this end, we have also calculated the temperature-dependent Gibbs free energies of relevant binary (LaP, LaP2, LaSi, LaSi2, SiP, SiP2) as shown in Fig. S6 in the ESI.† As long as the Gibbs free energies of the relevant phases at each temperature are obtained, Ghull of each ternary phase with respect to the ternary convex hull at each temperature can be calculated in the same way as Ehull calculation, by simply replacing the total energies of the competing phases with the corresponding Gibbs free energies G of the phases.
The calculated Ghull as the function of temperature for the five La–Si–P ternary phases are plotted in Fig. 6. We can see that two of the known phases, LaSiP3 and LaSi2P6 are the most stable ones at low temperature below 1100 K. However, La2SiP3 becomes more stable than the LaSi2P6 above 1200 K, and same for La2SiP4 above 1800 K. Since the La2SiP4 has been recently successfully synthesized,34 we expect that the newly discovered La2SiP3 phases, which is more thermodynamically favorable than La2SiP4 phase at all temperatures below 2000 K, can be synthesized as well. The results from this Gibbs energy analysis are important for understanding the thermodynamic stability of the compounds to guide experimental synthesis.
Fig. 6 Calculated Ghull as a function of temperature for the four competitive La–Si–P ternary phases. |
Moreover, GA searches using the developed ANN-ML interatomic potential also reveal additional new metastable La–Si–P phases as compared to the results from the CGCNN-ML. In the right panel of Fig. 8a and b, we show new metastable structures of the LaSiP3 and La2SiP4 phases, respectively, obtained from our GA search. In comparison with the lowest-energy structures of the two phases (which were also captured by our GA search) in the left panels of the figure, the new metastable structures have formation energies very close to those of lowest-energy structures (within 50 meV per atom). These new metastable structures could be accessible in experiment at higher temperatures. It is interesting to note that the metastable LaSiP3 structure is very similar to the lowest-energy one with the same noncentrosymmetric Pna21 space group, except the subtle difference in the orientation of the 1D P cis–trans chains. In the lowest-energy structure, all P-chains run parallel to the axis a direction, while in the metastable structure, the P-chains run along the a and b directions alternatively. In the La2SiP4 structures, the difference between the lowest energy and metastable structures is the relative orientation of the SiP4 tetrahedra (shown in blue in the figure).
Similarly, we calculated the Gibbs free energies as the function of temperature for the BaLaSiP3 quaternary phase, as well as the relevant binary and ternary phases (LaP, LaSiP3, Ba3P2, Ba3(Si2P3)2). The results are shown in Fig. S6 and S8 respectively in the ESI.† Using these free energy data, we calculated Ghull as the function of temperature for the BaLaSiP3 quaternary phase and plotted it in Fig. S9.† We can see that BaLaSiP3 phase stays below the convex hull at all temperatures below 2000 K, which we expect it can be synthesized as well.
Using the La–Si–P ternary system as a proof-of-concept, we show that, needing only 29 structures to be checked by first-principles calculations, the integrated ML approach correctly captures the three known ternary phases (La(SiP3)2, and ordered LaSiP3, La2SiP4) and predicted a novel stable La2SiP3 phase in the P-rich (P ≥ 50%) regime. The GA search using the developed ANN-ML interatomic potential quickly captured the three stable La–Si–P ternary phases and reveal several new low-energy metastable phases. The described approach is generic and robust which can be readily applied to any compounds of interest. Finally, combining these algorithms with chemical intuition, we expanded this work to explore quaternary Ba–La–Si–P phases by substituting half of the La atoms in the newly predicted La2SiP3 phase to yield an electron-balanced phase, BaLaSiP3, and probe its stability.
Our band structure calculations as shown in Fig. S5 and S7† indicate these two compounds are metallic but with only small number of bands across the Fermi level. Especially, the bands across the Fermi level in the La2SiP3 compound exhibit linear-like dispersion, suggesting that electrons in this compound would have high mobility. In order to propose potential applications of the newly discovered compounds (La2SiP3, BaLaSiP3 and other metastable phases shown in the ESI†), more evaluation and understanding on the properties of these compounds will be needed. This will be an interesting topic of future investigations.
For ordered crystalline compounds, it has been widely shown that vibrational entropy plays an important role in determining the relative thermodynamical stability and phase transition between different competing phases as the temperature varies.47–49 Our calculation results on the temperature effects as shown in Fig. 6 also support this point of view. For most materials, including the La–Si–P ternary compounds studied in this paper, the vibrational entropy can be accurately calculated by first-principles DFT calculations. Therefore, the method for calculating temperature-dependent Ghull to determine the temperature effects on the relative stability of the competing compounds described in this paper can also be applied to other compounds.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ta03771a |
This journal is © The Royal Society of Chemistry 2023 |