Introduction
Long-term storage and safe disposal of toxic wastes (nuclear, radioactive, or otherwise) is a pressing environmental concern that must be addressed for a sustainable future. Glasses and crystalline ceramic materials represent some of the common inorganic host materials that have been explored and have shown promise for the immobilization of waste materials.[Reference Weber, Ewing, Catlow, de la Rubia, Hobbs, Kinoshita, Matzke, Motta, Nastasi, Salje, Vance and Zinkle1,Reference Donald, Metcalfe and Taylor2] One of the waste forms that is generated from nuclear power plants is the toxic and radioactive iodine (I-129)[Reference Watanabe, Ikoma, Yamada, Suetsugu, Komatsu, Stevens, Moriyoshi and Tanaka3], and there is a growing priority to capture radioactive iodine, because unattended I-129 entering the ecosystem is a severe concern due to its long half-life, high solubility in water, and potential human ingestion.[Reference Masuda, Kato, Tanaka, Nakajima, Okamoto, Sakuragi and Yoshida4] More recently, materials that form an “apatite” crystal structure type have been considered as a potential host lattice for immobilizing the I-129 isotope.[Reference Watanabe, Ikoma, Yamada, Suetsugu, Komatsu, Stevens, Moriyoshi and Tanaka3,Reference Riley, Vienna, Strachan, McCloy and Jerden5]
Compounds with the apatite crystal structure type have a general formula A10(BO4)6X2, where the A-site is occupied by larger monovalent (Na+, K+, etc.), divalent (Ca2+, Sr2+, Ba2+, Pb2+, etc.), or trivalent (La3+, Y3+, Ce3+, Sm3+, etc.) cations, the B-site is occupied by smaller metals and metalloids (P5+, As5+, V5+, Si4+, etc.), and the X-site is filled by halide, hydroxide, or oxide anions (F−, Cl−, Br−, I−, OH−, O2−, etc.). The aristotype or the parent crystal structure of an apatite has a hexagonal symmetry in the space group, P63/m (# 176), with a one-dimensional channel parallel to the c-axis. Crystallographically, there are two distinct A-sites, AI and AII, occupying the Wyckoff sites 4f and 6h, respectively. The crystal structure is made up of a network of regular BO4 tetrahedra that is connected to a AIO6 metaprism and an irregular AIIO6X1,2 polyhedral unit. The crystal structure is shown in Fig. 1. Apatites also have been shown to undergo structural-phase transformation to lower symmetry structures as a function of chemical substitutions, defects (vacancies and interstitials), and temperature.[Reference White and ZhiLi6,Reference White, Ferraris, Kim and Madhavi7] Several desirable properties such as structural flexibility and chemical diversity, durability against radiation,[Reference Kim, Dong and White8] and accessibility to straightforward synthesis techniques justify the choice of the material.[Reference White and ZhiLi6,Reference White, Ferraris, Kim and Madhavi7,Reference Wopenka and Pasteris9,Reference Ardanova, Get'man, Loboda, Prisedsky, Tkachenko, Marchenko, Antonovich, Chivireva, Chebishev and Lyashenko10]
Crystal structures of A10(BO4)6X2 apatites in the aristotype P63/m space group projected onto the (0001) plane. The structure is made of BO4 tetrahedral units that are connected to a sixfold coordinated AIO6 metaprism unit. In addition, an irregular AIIO6X1,2 polyhedral unit circumscribes the one-dimensional channel parallel to the c-axis.
Despite the progress in improving the properties of iodine-containing apatites (iodoapatites) such as the loading capacity, leaching rate, thermal recovery, and mechanical properties via processing routes,[Reference Stennett, Pinnock and Hyatt11,Reference Yao, Lu, Sun, Wang, Ewing and Lian12] most of the experimental studies have focused only on nonstoichiometric forms of lead vanadate [Pb10(VO4)6I2] iodoapatites. The quest for iodoapatites with better properties is mainly limited by the lack of exploration. With potentially many choices for chemical substitutions in both A- and B-sites in the apatite structure, experimental investigation for all possible combinations is not realistic. A priori knowledge of the compositional limit of a particular chemical constituent in the apatite lattice can provide important insights toward tailoring new compounds with desired properties.
Here, we develop a computational approach built on the foundations of machine learning (ML) and density functional theory (DFT) methods to predict possible new iodoapatites. The idea of combining ML methods with DFT calculations (including the high-throughput variant) to predict the formability and thermodynamic stability is not entirely new. There is a growing body of literature in exploring the potential of ML methods for predicting compound formability. Some of the notable contributions include the prediction of formability of perovskite-based compounds, Heusler alloys, high entropy alloys, and metallic glasses (to name a few).[Reference Legrain, Carrete, van Roekeghem, Madsen and Mingo13–Reference Ren, Ward, Williams, Laws, Wolverton, Hattrick-Simpers and Mehta21]
In terms of prior ML studies on iodoapatites, Wang [Reference Wang22] built an artificial neural network model to predict the channel size of apatite structure as a function of iodine incorporation. The hypothesis was that large channels host large anions (such as I−) and small channels host small anions (such as F−). As a result, the goal was to identify A- and B-site elements that will lead to large channels so that iodine can be incorporated. Wang neither considered the criterion of formability nor thermodynamic stability in the published work. One of the plausible reasons could be because of the lack of databases in the literature where such formability data are compiled. We bridge the gap in this paper by constructing a database from surveying the experimental literature and demonstrate the utility of ML methods by predicting novel and formable iodoapatite compounds.
In this work, we employed ensemble ML models that were based on tree-based methods (random forests and gradient tree boosting) to establish a relationship between the crystal chemistry descriptors and apatite-phase formability. Ensemble models combine predictions from 100 or 1000 s (or more) of ML models and are shown to be more accurate than predictions from a single model.[Reference Optiz and Maclin23] Further, the apatite phase formability problem is compounded by the small nature of the training dataset sizes (336 data points). One of the challenges in the application of ML methods on a small dataset is its high variability. Ensemble methods provide a means to address the bias–variance tradeoff by reducing the effect of the variance term. In addition to ML methods, we also employed DFT calculations to calculate the thermodynamic stability of apatite compounds. We utilized the convex hull (CH) construction as implemented in the open quantum materials database (OQMD)[Reference Saal, Kirklin, Aykol, Meredig and Wolverton24,Reference Kirklin, Saal, Meredig, Thompson, Doak, Aykol, Rühl and Wolverton25] and calculated the decomposition energy (ΔE D) associated with each apatite compound. The ΔE D cutoff was determined from comparing the DFT data with experimental XRD (x-ray diffraction)-phase formation data. Finally, we used our trained ML models to predict the formability of previously unexplored iodoapatites. These results were augmented with thermodynamic stability predictions from DFT, especially for those compounds where the two ML models disagreed. The key outcome is the prediction of 18 novel iodoapatites as formable, two as thermodynamically stable, and one as metastable. We recommended these promising predictions for experimental synthesis and validation.
Methods and approach
Machine learning
Dataset
The first step in the ML approach is to develop a representative dataset of known observations. For developing the classification models, we have considered a total of 336 compositions taken from the literature.[Reference White, Ferraris, Kim and Madhavi7,Reference Ito26–Reference Silva, Vasconcelos, Sombra and Valente33] In Figure 2, a schematic overview of our dataset is given. The chemical compositions that contain fluoride, chloride, oxide, and hydroxide anion in the X-site are well represented in our dataset, whereas those containing bromide and iodide are relatively sparse. For the given chemical composition, we distinguished between stable apatites and nonapatites based on the peaks assigned to the apatite structure from the XRD pattern. For example, in our database, we assigned stable apatite structure to Ca10(PO4)6F2, whereas for Zn10(PO4)6F2 we assigned it to be a nonapatite, because the synthesis of fully stoichiometric Zn10(PO4)6F2 apatite was not shown to be experimentally feasible.[Reference Grisafe and Hummel28] In addition, when an experimental synthesis resulted in the formation of secondary phases, we assigned it to be nonapatites. Instead of using the entire XRD pattern in our analysis, we defined a single categorical variable called “apatite formability” and assigned a categorical label “Yes” to the chemical compositions that are experimentally determined to have a single-phase apatite structure and “No” to the compositions that do not form single-phase apatite structure. Of the 336 compositions in our database, 224 had formable apatite structure and 112 were nonapatites. It should also be noted that in our database, we do not explicitly differentiate between the different crystal symmetries of apatites. Apatites are known to exist in P63/m, P63/cm, P63, $P\bar {3}$
, P21/b, P21, $P\bar {6}$
, and $P\bar {1}$
[Reference White, Ferraris, Kim and Madhavi7,Reference Baikie, Mercier, Elcombe, Kim, Page, Mitchell, White and Whitfield34] space groups, and the label “Yes” can correspond to any of these space groups. The dataset is given in Supplementary Material.
Schematic overview of various A-site, B-site, and X-site elements considered in this study. The choice of these chemical constituents was based on the available data from the past experimental studies. The numbers in the pie chart indicate the total count of chemical compositions that contains the specific X-site element. For example, in our data set, there were 114 chemical compositions that had fluoride anion in the X-site. Overall, we considered a total of 336 chemical compositions, out of which 224 had state apatite structure and 112 were nonapatites. In addition, the figure also describes the specific combination of various A-site and B-site elements used for a particular X-site element.
Past experimental and theoretical studies have indicated the importance of three descriptors: ionic radii, electronegativity, and polarizing power, in describing the apatite-phase formation.[Reference Suzuki, Hatsushika and Miyake35–Reference Balachandran and Rajan39] While ionic radii can be related to the lattice energy of solids,[Reference Muller and Roy40] the physical meaning of electronegativity and polarizing power (a relative measure of the power of a cation to attract an anion) [Reference Ahrens41] takes into account the bond directionality and charge density information. The polarizing power of a cation is defined as the ratio of the charge of the ion to its ionic radius. Therefore, by combining the three parameters, we hypothesize that the electronic structure and bonding of complex apatites can be well approximated. The full list of variables used in this work to quantify each stoichiometry is given in Table I. The descriptors of Shannon' ionic radii,[Reference Shannon42] Pauling's electronegativity[Reference Pauling43], and polarizing power of the cation were used as the predictor (or independent) variables, and the phase formability information (“apatite formability” variable) extracted from XRD analysis was used as the dependent variable. We used the group electronegativity value of 3.7 for the OH anion that was scaled to be compatible with the Pauling scale.[Reference Anslyn, Dougherty and Murdzek44] In Supplementary Fig. 1, the numerical values of the predictor variables used in the analysis are provided. Using these values, the average ionic radii, average electronegativity, and average polarizing power of the cation at each crystallographic site were calculated by multiplying each variable by its weight, adding these values and dividing the total by the total weight of A-, B-, and X-sites separately. For example, the average ionic radii of the A-site (〈r A〉) in A$^1_p$
A$^2_q$
(BO4)6X2 is calculated using the following equation:
where, $r_{\rm A^1}$
and $r_{\rm A^2}$
are the ionic radii of A1 and A2 chemical elements that occupy the A-site of the lattice, respectively, and p and q are the relative weights of A1 and A2 chemistries, respectively.
List of descriptors used for training the model.
All descriptors must be interpreted in terms of their average values as calculated based on Eq. (1).
The logic explained above in Eq. (1) was repeated in calculating the average values of all the predictor variables listed in Table I. The 13 crystal chemistry descriptors in Table I served as the input for ML. Since the response (or output) is a categorical, binary variable (Yes or No), the type of ML analysis is referred to as classification learning. The goal of ML is to establish a relationship between the response and the inputs. We utilized random forests and gradient tree boosting for classification learning, which are discussed next.
Random forests
Random forest is an ensemble learning method that generates a specified number of unbiased decision trees then aggregates the decisions of all trees to make a final prediction of the class label.[Reference Breiman45,Reference Liaw and Wiener46] These trees are generated using a bootstrap method where the training data are sampled, with replacement, to grow the forest. Most random forest algorithms typically withhold 30% of the training data during this process, although the number may vary slightly depending on the algorithm. These withheld data points are referred to as the out-of-bag (OOB) samples and are used to assess the accuracy of the trained model. This is a straight forward process since the true labels of the OOB samples are contained in the training set.
The algorithm works by first sampling a pre-specified number of bootstrapped samples (ntree), and larger values of ntree generally increase the accuracy of the model but are more complex to train and often do not increase accuracy by a substantial amount beyond a certain point. The algorithm then randomly samples, with replacement, a pre-specified number of input parameters (mtry) and generates a decision tree for each bootstrapped sample based on those parameters. We utilized the random forest method as implemented in the randomForest R-package.[Reference Liaw and Wiener46]
Gradient tree boosting
Gradient tree boosting is another ensemble learning method that uses decision trees as base learners.[Reference Friedman47,Reference Friedman48] In contrast to bagging algorithm of the random forests, the boosting algorithm learns each decision tree sequentially. The first tree is used to coarsely segment the observations in the data. Then, the subsequent trees were fit in such a way that the error of the current model was further minimized. The loss function used in the model is a binomial log-likelihood function. The final prediction of the class is estimated as the sum of predictions of the individual trees. There are three tuning parameters: (i) the number of trees, (ii) the maximum number of splits in each tree, and (iii) the learning rate of the gradient of loss function. We fixed the maximum number of splits in each tree and the learning rate hyperparameters to the default values. We tuned the number of trees using the k-fold cross-validation method. We utilized the gradient tree boosting method as implemented in the gbm R-package.[Reference Greenwell, Boehmke, Cunningham and Developers49]
Density functional theory
DFT calculations were performed within the generalized gradient approximation (GGA) as implemented in Quantum ESPRESSO.[Reference Giannozzi, Baroni, Bonini, Calandra, Car, Cavazzoni, Ceresoli, Chiarotti, Cococcioni, Dabo, Dal Corso, de Gironcoli, Fabris, Fratesi, Gebauer, Gerstmann, Gougoussis, Kokalj, Lazzeri, Martin-Samos, Marzari, Mauri, Mazzarello, Paolini, Pasquarello, Paulatto, Sbraccia, Scandolo, Sclauzero, Seitsonen, Smogunov, Umari and Wentzcovitch50] The PBEsol exchange-correlation functional[Reference Perdew, Ruzsinszky, Csonka, Vydrov, Scuseria, Constantin, Zhou and Burke51] was used, and the core and valence electrons were treated with ultrasoft pseudopotentials. The Brillouin zone integration was performed using a Monkhorst-Pack k-point mesh centered at Γ and a 60-Ry plane-wave cutoff for wavefunctions (600 Ry kinetic energy cutoff for charge density and potential). The scalar relativistic pseudopotentials were taken from the PSLibrary.[Reference Corso52] The atomic positions and the cell volume were allowed to relax until an energy convergence threshold of 10−8 eV and Hellmann–Feynman forces <2 meV/Å, respectively, were achieved. We considered the following space groups to determine the lowest energy apatite structure: (i) in the case of fluorapatites (X=F), hexagonal (P63/m), and triclinic $(P\bar {1})$
space groups, (ii) in the case of chlorapatites (X=Cl), bromapatites (X=Br), and iodoapatites (X=I), hexagonal (P63/m), monoclinic (P21/b), and triclinic $(P\bar {1})$
space groups, and (iii) in the case of hydroxyapatites (X=OH), monoclinic (P21/b and P21) and hexagonal (P63) crystal symmetries. The choice of these symmetries was motivated by the earlier work in the literature.[Reference White, Ferraris, Kim and Madhavi7,Reference Balachandran, Rajan and Rondinelli53] We did not consider substitutional defects, vacancies, or interstitials in our calculations.
To evaluate the thermodynamic stability and decomposition reaction of each apatite compound, we considered the large pool of compounds assembled in the OQMD.[Reference Saal, Kirklin, Aykol, Meredig and Wolverton24,Reference Kirklin, Saal, Meredig, Thompson, Doak, Aykol, Rühl and Wolverton25] To construct the CH for a given A–B–O–X element set, the OQMD uses the total energy data computed by DFT of all (both existing and hypothetical) elements and compounds in the OQMD. Due to computational challenges involved in treating solid solutions, DFT calculations were only run on a subset of the total training data that represented simple compounds (Supplementary Table 1). For simplicity, all magnetic materials and those containing partially filled d-orbitals in the divalent (Fe2+, Mn2+, and Ni2+) or pentavalent charge states (Fe5+, Mn5+, Cr5+, Re5+, etc.) were not considered in this analysis. In total, the ΔE D values were computed for 38 apatite compounds in the training set. We extracted the reactants from the OQMD and performed spot DFT calculations. Our product is the apatite phase in its lowest energy structure. We then calculated the distance from the CH as $\Delta E^{{\rm D}}= E ^{{\rm product}}-\sum _{i=1}^{n} E^{{\rm reactants}}_i$
(where i = 1, …, n represent the distinct reacting species and E is the total energy from DFT calculations), which can take either a positive [the hypothetical A10(BO4)6X2 compound is metastable or unstable] or a negative value [the hypothetical A10(BO4)6X2 compound is more stable than other compounds that lie on the CH, and thus, the apatite is more likely to form].
Results and discussion
The hyperparameters for the random forest models were optimized by iterating over a grid space of ntree and mtry. We explored a range of 1–2000 trees for ntree and 1–13 for mtry. The optimal OOB error rate was calculated to be 13.1% for the combined ntree and mtry values of 955 and 3, respectively. As a result, the final predictions about apatite formability were made by combining predictions from 955 models. Similarly, the hyperparameters for the gradient tree boosting (number of trees) were optimized using a k-fold cross-validation, where we explored several values of k-fold. The best result was obtained for a k-value of 15 (number of trees = 500). In this case, the final predictions were made by combining predictions from 500 models. In order to reliably compare random forest and gradient tree boosting methods, we randomly split our dataset into training and test sets. A 80%–20% training set–test set split was used to build and test the trained ML models, respectively. We also repeated the process five times with five different random seeds, and the average test set error rate from the five random trials was calculated. The random forest and gradient tree boosting models performed with an average test set error rate of 15% and 18%, respectively. These results indicate that the random forest models perform relatively better compared with the gradient tree boosting models on this dataset.
In Supplementary Table 1, the results from DFT-CH analysis are given for 38 compounds for which apatite-phase formation data are known from experimentally measured XRD data. Ideally, we expect the DFT-CH data to consist of two independent probability distributions, one corresponding to experimentally formable and the other to nonformable compounds. The violin plots of this data appear to show that the distributions are independent (Fig. 3). To verify independence statistically, we conducted a nonparametric Mann–Whitney U-test (also referred to as the Wilcoxon–Mann–Whitney test or the rank-sum test). This nonparametric test is commonly used to determine the independence of any two distributions,[Reference Fay and Proschan54] where a P-value of <0.05 generally means that two samples are from distributions with separate median values. The P-value for the test was 2.04 × 10−5, confirming that these were indeed independent distributions.
Distributions of formable and nonformable compounds based on their DFT-CH predicted ΔE D. One of the nonformable compounds, Ga10(PO4)6F2, is identified as an outlier in the violin plot because of an uncharacteristically large ΔE D in the DFT-CH analysis. The degree of metastability (DOM) is set at +16 meV/atom, and any compound with ΔE D ≤ DOM are predicted to be potentially formable.
Comparing several different cutoff energies, there is some ambiguity in the area of formability between the values of −10 to +16 meV/atom. General thermodynamical stability intuition would lead us to believe that compounds with a positive ΔE D would be energetically unfavorable. For this reason, we defined a “degree of metastability (DOM)” and referred to compounds whose ΔE D falls in this range (≥0 and ≤ +16 meV/atom) to be metastable and potentially formable using nonequilibrium synthesis methods. Additional details can be found in Supplementary Section III.
On the basis of the DOM analysis, we recommend re-evaluation of Pb10(PO4)6I2 compound. Cao et al. [Reference Cao, Chong, Thirion, Mauro, McCloy and Goel55] used wet chemical synthesis methods and led to the conclusion that the incorporation of (PO4) in Pb10(VO4)6I2 had a compositional limit (4.77 wt%), beyond which the single-phase structure was not experimentally observed. Therefore, we labeled Pb10(PO4)6I2 as nonformable in our training set. However, the DFT-CH analysis predicts a ΔE D value of +5.14 meV/atom (only weakly metastable) for Pb10(PO4)6I2. Therefore, nonequilibrium processing methods can potentially stabilize an apatite phase in Pb10(PO4)6I2. Further, random forest and gradient tree boosting methods predicted that Pb10(PO4)6I2 is nonformable with a probability of 51% and 87%, respectively. Therefore, the random forest models are uncertain about assigning a label for Pb10(PO4)6I2, which further indicates that the compound should be re-evaluated. On the other hand, gradient tree boosting is fairly confident (87%) to assign a nonformable label to Pb10(PO4)6I2.
Having trained the ML models using formable and nonformable experimental data and having established a thermodynamic stability and DOM criteria (ΔE D) using known data from DFT-CH analysis, we now predict the formability of previously unexplored iodoapatites. Our strategy is shown in Fig. 4. A dataset was constructed by including a list of potential iodoapatites considering possible permutations of several different atoms in the A- and B-sites. We included Ba, Ca, Cd, Ga, Pb, and Sr as possible A-site atoms and As, Mo, Ga, Nb, P, Re, Os, Ta, V, and W as B-site atoms. The combinations of these A- and B-site atoms gave 60 potential iodoapatite stoichiometries, and six have been already experimentally synthesized. Out of the remaining 54 compounds, 42 agreed between random forest and gradient tree boosting, which included predictions of 18 compositions as formable and 24 as nonformable. The list of 18 formable novel apatites is given in Table II. For the 12 compounds where random forest and gradient tree boosting predictions disagreed (i.e., gradient tree boosting predicted them to be formable, whereas random forest predicted them to be nonformable), we performed DFT-CH analysis to provide additional insights into the thermodynamic stability. The results are given in Table III. From the data, it can be seen that the thermodynamic stability analysis inclines toward the random forest predictions.
The strategy to predict the formability and thermodynamic stability of previously unexplored or hypothetical 54 A10(BO4)6I2 using two ML methods (random forest and gradient tree boosting) and DFT-CH thermodynamical stability analysis. The ML models predict 18 as formable and 24 as nonformable. There were 12 compounds for which the outcomes from ML models were uncertain (gradient tree boosting predicted them to be formable, but random forest predicted them to be nonformable). We performed DFT-CH analysis on those 12 uncertain compounds to gain insights into their thermodynamic stability, which predicted two to be stable, one to be metastable, and nine to be unstable.
List of 18 novel and formable iodoapatites predicted by random forest and gradient tree boosting.
The ΔE D (meV/atom) from DFT-CH analysis are also given for a few compounds in parentheses. Any compound with ΔE D ≤ +16 meV/atom will be predicted as formable in the apatite structure, whereas compounds with ΔE D > +16 meV/atom will predicted as nonformable from the DFT-CH analysis.
DFT-CH predictions for 12 novel iodoapatites, where random forest and gradient tree boosting methods disagree in their predictions.
ΔE D (in meV/atom) is the decomposition energy or the distance from the CH. When ΔE D ≤ 0, then the compound is predicted to be thermodynamically stable. When ΔE D > +16, then the compound is predicted to be unstable. The inequality condition, 0 < ΔE D ≤ +16, represents the degree of metastability (DOM) region.
We predicted two compounds to be thermodynamically stable with ΔE D < 0, namely [Ba10(NbO4)6I2 and Ba10(TaO4)6I2]. However, the calculated ΔE D values were uncharacteristically small. Plausible reasons include mismatch in the choice of DFT codes. For instance, our DFT calculations were performed using an ultrasoft pseudopotential and PBEsol-GGA functional in Quantum ESPRESSO, whereas the decomposition energies were taken from the OQMD that uses projector-augmented wave pseudopotential and PBE-GGA functional in VASP (The Vienna Ab initio simulation package). Other source for discrepancy can come from the consideration of only relatively fewer or nonrepresentative compounds in the construction of CH for the Ba–Nb–O–I and Ba–Ta–O–I element sets. More discussions on this topic can be found in Ref. [Reference Balachandran, Emery, Gubernatis, Lookman, Wolverton and Zunger15]. For consistency, we also checked whether Materials Project[Reference Jain, Ong, Hautier, Chen, Richards, Dacek, Cholia, Gunter, Skinner, Ceder and Persson56] can predict the decomposition reactions for Ba10(NbO4)6I2 and Ba10(TaO4)6I2. However, our query returned with “No data available.” This is likely because there is no experimental data on Ba10(NbO4)6I2 and Ba10(TaO4)6I2 anywhere in the published literature. Therefore, we recommend the readers to interpret the ΔE D for Ba10(NbO4)6I2 and Ba10(TaO4)6I2 in Table III with cautious optimism. Furthermore, the random forest models predicted that there is only 32% and 28% probability of forming an apatite phase in Ba10(NbO4)6I2 and Ba10(TaO4)6I2, respectively.
We now turn our attention to Sr10(PO4)6I2 in Table III that has a ΔE D value of +8.03 meV/atom. The ΔE D is below the DOM threshold of +16 meV/atom. Thus, we predict Sr10(PO4)6I2 to be metastable and potentially formable. The random forest models, on the other hand, predicted that Sr10(PO4)6I2 has a 44.3% probability of forming an apatite phase. These results can be used as guidelines for the synthesis of novel iodoapatite compounds.
We also performed DFT-CH calculations on six hypothetical compounds that were predicted to be formable by both random forest and gradient tree boosting methods (see Table II): Ca10(AsO4)6I2, Sr10(VO4)6I2, Ba10(PO4)6I2, Sr10(AsO4)6I2, Ba10(AsO4)6I2, and Ba10(VO4)6I2. Five of the compounds have ΔE D ≤ +16 meV/atom, indicating that they are only weakly metastable and potentially formable. The exception is Ba10(AsO4)6I2 that has a calculated ΔE D of + 21.27 meV/atom. Since the value of ΔE D for Ba10(AsO4)6I2 is > +16 meV/atom (DOM), the prediction from DFT-CH analysis is that Ba10(AsO4)6I2 is not formable in the apatite crystal structure. This result is in disagreement with the ML predictions. All compounds have the monoclinic space group (P21/b) as the lowest-energy crystal structure.
Finally, we compare our predictions with the work of Wang[Reference Wang22] who predicted that when the A-site is Sr2+, Pb2+, Cd2+, and Ba2+ and when the B-site is V5+, As5+, Cr5+, and Mo5+, then the structure will have the suitable channel size for incorporating iodine in the lattice. Our ML predictions generally agree with Wang's work, when A={Sr2+ and Pb2+} and B={V5+, As5+, and Mo5+} (see Table II). Our work also predicts new and formable Ca2+-based iodine-containing apatites that were missed by Wang. However, there were disagreements when considering the Cd-containing compounds, where DFT-CH analyses predict some of the Cd-containing iodoapatites to be thermodynamically unstable (see Table III). Experimental synthesis and characterization will be needed to provide validation and feedback.
Conclusion
Predicting the formability and thermodynamic stability is one of the important first steps for assessing the potential of realizing any compound as a waste form for immobilizing toxic elements, such as radioactive iodine. We have made progress by exploring two complementary computational methods, such as data-driven ML and physics-based DFT, to address the problem. ML methods predict 18 new iodoapatites as formable. In addition, there were 12 compounds that posed a challenge to ML methods. We performed DFT-CH analysis on those 12 compounds and predicted two to be thermodynamically stable and one to be metastable (DOM ≤ +16 meV/atom). The remaining nine were predicted to be unstable. We recommend the promising formable new compounds for experimental synthesis and validation. In addition to line compounds, the ML models also have the capability to rapidly predict the formability of complex solid solutions, which are beyond the scope of the current state-of-the-art DFT methods. Although the analysis and results presented focused on the formability and stability of iodine-containing apatites, our ML models are capable of exploring other chemical substitutions, including toxic Hg.[Reference Wang, Kim, Dionysiou, Sorial and Timberlake57] Furthermore, our computational approach is generic and can seamlessly be adapted to other crystal systems (e.g., lanthanide phosphates[Reference Toro-González, Clifford, Copping, Mirzadeh and Rojas58]) to address problems in environmental remediation and related applications, where phase formability and stability are important for design considerations.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1557/mrc.2019.103.
Acknowledgments
All authors acknowledge support from the Advanced Research Computing Services at the University of Virginia (UVa), which maintains the Rivanna supercomputing cluster where we performed the DFT calculations. M.V.A and P.V.B acknowledge support from the UVa start-up funds. T.Q.H acknowledges support from the DARPA-Topological Excitations in Electronics program (grant# D18AP00009). The content of the information does not necessarily reflect the position or the policy of the Government, and no official endorsement should be inferred. Approved for public release; distribution is unlimited. All authors thank Dr Shruba Gangopadhyay for providing insightful comments on the manuscript.
