Machine Learned Interatomic Potentials for Ternary Carbides trained on the AFLOW Database

Large density functional theory (DFT) databases are a treasure trove of energies, forces and stresses that can be used to train machine learned interatomic potentials for atomistic modeling. Herein, we employ structural relaxations from the AFLOW database to train moment tensor potentials (MTPs) for four carbide systems: HfTaC, HfZrC, MoWC and TaTiC. The resulting MTPs are used to relax ~6300 random symmetric structures, and are subsequently improved via active learning to generate robust potentials (RP) that can relax a wide variety of structures, and accurate potentials (AP) designed for the relaxation of low-energy systems. This protocol is shown to yield convex hulls that are indistinguishable from those predicted by AFLOW for the HfTaC, HfZrC and TaTiC systems, and in the case of the MoWC system to predict thermodynamically stable structures that are not found within AFLOW, highlighting the potential of the employed protocol within crystal structure prediction. Relaxation of over three hundred Mo$_{1-x}$W$_x$C stoichiometry crystals first with the RP then with the AP yields formation enthalpies that are in excellent agreement with those obtained via DFT.


INTRODUCTION
Machine learned interatomic potentials (ML-IAPs), which are trained on density functional theory (DFT) data, have irrevocably changed the way in which computational materials science is carried out.They have increased the complexity of the systems that can be studied via static calculations [1], the duration of molecular dynamics trajectories [2], and made it possible to routinely model effects, such as anharmonicity [3], which are commonly neglected due to the large computational cost associated with the requisite DFT calculations.Moreover, they are becoming increasingly important in the field of crystal structure prediction (CSP) where it may be necessary to relax hundreds or thousands of structures to find the global minimum and a handful of low-lying local minima [4].As CSP methods tackle systems with increased combinatorial complexity such as ternaries and quaternaries [5], it is becoming increasingly important to develop protocols that can be used to train and employ ML-IAPs for CSP.
Some of the most well-known ML-IAPs include neural networks (NNs) [6,7], the spectral neighbor analysis potential (SNAP) [8], moment tensor potentials (MTPs) [9], the Gaussian approximation potential (GAP) [10], and, more recently ultra-fast (UF) [11], and ephemeral data-derived (EDD) potentials [12].Provided an ML-IAP can accurately describe the potential energy surface (PES) of a multi-component system, it can significantly accelerate a CSP search while at the same time enabling exploration of wider regions of compositional space [4].Indeed, recent algorithmic developments have shown promising results: CSP methods including random and evolutionary searches, as well as particle swarm optimization have been coupled with various ML-IAPs and applied to predict crystalline structures of boron [13][14][15], carbon [14,16], sodium [14], phosphorus [17], lithium [18] as well as Mg-Ca [19] and various metal-tin [20,21] alloys at 1 atm and under pressure.For ternary crystalline systems, early work showed that adaptive classical potentials can accelerate CSP in the Mg-Si-O system [22], and recently ephemeral deep learning potentials combined with random searching have been used to predict candidate structures for a reported Lu-N-H superconducting phase [23], and Zn(CN) 2 metal-organic-frameworks [24].
A challenge in developing ML-IAPs for CSP is that they need to be able to predict the energies of unstable structures, as well as those that are close to the local minima in the PES.Moreover, because ML-IAPs are poor at extrapolation, they cannot accurately predict the energies (or forces or stresses) of configurations that differ significantly from those they have seen before.Traditionally, DFT data sets containing thousands to hundreds-of-thousands of structures have been generated for the training and testing of ML-IAPs [14,25].A wide range of procedures have been employed to create these DFT data sets including generating structures randomly [14], perturbing the geometries of such structures via 'shaking' [12], ab initio molecular dynamics runs at various temperatures [26], decorating predefined lattices with atoms of different types while simultaneously varying the chemical composition [25], relaxation of structures generated via constrained evolutionary searches [7], straining crystal lattices, creating defect structures, and more [26,27].
Unfortunately, even when large DFT data sets are employed for training, it is unlikely that the resulting ML-IAPs can predict, with sufficient accuracy, the energies of the various structures encountered in the course of a CSP run.One strategy that has been proposed for the generation of a multi-purpose ML-IAP, given a limited number of DFT calculations, relies on the automatic iterative building of the fitting database by selecting the most diverse structures [28].Additional techniques include various active learning [29,30] and learning-on-thefly [31] methods, where the potential is updated and improved during the course of the search.It has been suggested that active learning could be used to generate two ML-IAPs: a ro-bust one that is able to optimize any structure the CSP algorithm encounters and make rough predictions, and an accurate one trained on (and able to relax) only the low-energy structures [14,25].
ML-IAP based simulations where the potential is updated on-the-fly may be initialised using either an empty/untrained potential [25], or one that has been pre-trained.The former strategy may require just as many, if not more, DFT evaluations than the latter because the likelihood of encountering a configuration that is deemed extrapolative is high, necessitating a retraining of the ML-IAP [14].Thus, there appears to be 'no-free-lunch' since the construction of a reliable ML-IAP requires numerous expensive DFT calculations.At the same time, a large number of databases exist -AFLOW (Automatic FLOW) [32][33][34], the Materials Project (MP) [35], the Open Quantum Materials Database [36], etc. -each containing millions of DFT evaluations of the energies, forces and stresses of extended systems.Furthermore, it is becoming standard practice for researchers to deposit the DFT data generated during the course of a computational project in repositories such as NOMAD [37], OCELOT [38] and NIST Materials Data [39].One way this data is being used is to train ML-IAPs for the computational study and exploration of the vast PES of all possible chemistries.Some of the forefront examples of such 'universal' ML-IAPs, which can predict energies, forces and stresses using equivariant graph neural networks, include M3GNet [40], CHGNet [41], ALIGNN-FF [42], MACE-MP-0 [43], and GNoME [44].
The training of potentials on already existing DFT-data is illustrated here by combining outputs present within the AFLOW [32] database with MTPs [9], as implemented within the Machine Learned Interatomic Potentials (MLIP) program [45].Specifically, chemically sensible structures, which are randomly selected from the relaxation trajectories stored within the AFLOW database, are used to train an MTP that is subsequently employed to relax a large number of random symmetric structures spanning a wide composition range.In contrast to the recently developed universal ML-IAPs [40][41][42][43][44]46], here we only train on a subset of the data found within AFLOW, chosen for the application in-mind.Therefore, the ML-IAPs developed here are system-specific, and not universal.What distinguishes our study from prior works [7,12,14,26,27] is that rather than generating our own DFT-training set, we employ already existing data found within AFLOW.In a further step, the AFLOW-trained potentials are improved via active learning, generating ML-IAPs that can be robust (able to roughly optimize any configuration) or accurate (for the optimization of only low-energy structures near the convex hull).Thus, only a small number of supplementary DFT calculations are required to develop system-specific MLIPs that enable the computational exploration of evermore complex PESs towards the discovery of materials.
A utility package that automates this training process, the Plan for Robust and Accurate Potentials (PRAPs), is described.The method is used to determine the zero-Kelvin phase diagrams of four ternary metal carbides, chosen because they represent materials with superlative mechanical properties [47,48].The PRAPs generated convex hulls whose lowest energy structures are optimised with DFT are in good agreement with the DFT-relaxed AFLOW hulls.Moreover, in the case of the CMoW system, thermodynamically stable structures are found that are not present in the AFLOW data.Further calculations with the accurate potential find a variety of Mo 1−x W x C stoichiometry phases at/near the tie-line indicating the possibility of a solid-solution with a very low miscibility gap critical temperature.

RESULTS AND DISCUSSION
A. Plan for Robust and Accurate Potentials (PRAPs) We created MTPs of varying complexity for a number of ternary metal carbides and investigated their capabilities to predict structures outside of their training sets.Our choice of MTPs was motivated by their excellent balance between model accuracy and computational efficiency [26,49], their application towards multicomponent systems [25,27,[50][51][52], their ability to predict phonons and thermodynamic properties, [3,53,54] and the availability of a powerful active learning scheme (ALS) that has been interfaced with the MTP method [25,55].
The MLIP software package trains MTPs and uses them to relax the geometries and minimize the energies of a wide variety of chemical systems [45].The simplest form of training, basic training, employs the energy, force and stress (EFS) data of a set of configurations, as obtained from DFT or other quantum chemical calculations, to generate an MTP.The complexity of the MTP is described by a user-selected level, a notation containing information about the number of basis functions and parameters comprising the potential.The ALS employs a D-optimality criterion to calculate the extrapolation degree or grade, γ, for every structure that is generated throughout the course of the simulation (relaxation trajectory or molecular dynamics run) [25,55].MLIP automatically selects configurations to be added to the training set if their γ falls within a user-defined range; we choose the default 2 < γ ≤ 10.The calculation (relaxation or molecular dynamics run) is terminated if γ exceeds the upper bound, triggering retraining of the MTP, and the procedure repeats until the simulation finishes with γ ≤ 10.Further information about MTPs, including their functional form, the quantities included in the cost function, and details of the active learning procedure are provided in Section S1 of the Supplementary Information.
In what follows we give a brief overview of the PRAPs utility package employed in this study; a forthcoming manuscript will describe its composition and usage.PRAPs interfaces with MLIP and automates the creation of a Robust Potential (RP) and Accurate Potential (AP) using the aforementioned ALS, and employs basic training for supplemental tasks.From a given a set of configurations obtained from AFLOW, a subset of ∼800 configurations is chosen randomly and used to train an MTP (Figure 1, gray box).This procedure is repeated five times, generating five different MTPs (as recommended in the MLIP manual [45]), mimicking a cross-validation procedure.From these PRAPs finds the 'best' MTP, defined as the The AP-relaxed structures may be sent for relaxation via AFLOW-DFT, followed by subsequent analysis (green box).The inset schematically illustrates the energy distribution of the structures that can be relaxed with the RP and the AP.Optionally, structures generated using such a procedure may not be exact since the intermediate structures in a DFT relaxation trajectory may have erroneous energies or forces resulting from Pulay stress that arises when the lattice vectors are varied.another procedure (here RANDSPG [56]), labelled as 'Non-DFT Data', may be relaxed via the RP and/or the AP, and this auxillary set of structures also comprises the Relaxation Set.Rectangles with sharp (curved) corners represent structural data (potentials).
one that identifies the most high-and-low-energy structures (ten of each are compared), or if multiple MTPs fulfill this criteria, the MTP with the lowest root-mean-squared (RMS) training error.This step is intended to mitigate the initial random parameters that MLIP assigns at the beginning of training, but users may bypass it if they so desire.To ensure that training is performed on sensible structures, PRAPs can filter out configurations with undesirable interatomic distances (default), or cell volumes.If relaxation trajectories comprise the data set, dozens of ionic steps originating from the same system may be present.Though using many similar configurations for training may be beneficial in some cases, in others it may be desirable to exclude the intermediate steps and only include the final relaxed configuration, and PRAPs provides an option for users to select the desired behaviour.
The RP training (Figure 1, blue box) begins using the PRAPs-determined 'best' MTP.The set of configurations may be optionally augmented with structures lacking any EFS data (Figure 1, left-hand-side, no box), which are combined with the initial DFT data set to form the relaxation set.The relaxation set is used to train the RP by active learning.The active learning begins with the 'best' training set and MTP from the pre-training step (if performed, otherwise PRAPs trains from scratch).During the active learning process, the relaxation set is relaxed using MLIP's built-in functionality and new structures are added to the training set based on the D-optimality criterion (a detailed discussion of this step is available in Section S1).The final output is a potential capable of performing a reasonable relaxation of most any configuration in the relaxation set.
After relaxation with the RP, PRAPs filters out structures with energies above a specified cutoff (default 50 meV/atom) for each composition present.This creates the 'Robust Relaxed' set, which is employed to train the AP via active learning (Figure 1, orange box).The procedure for training the AP is similar to the RP except that PRAPs begins with an empty training set instead of the 'best' pre-trained MTP's training set.We find that this results in better predictions as the AP should not see non-ground-state structures during training.The resulting AP should improve the EFS predictions for these lowenergy structures.For additional efficiency, users can alter the convergence criteria of the AP training at the cost of a few meV/atom in the training error (as described more thoroughly below).
PRAPs performs a number of analyses throughout and at the end of the training procedure.In addition to calculating the training and prediction RMSE and mean absolute error (MAE), and comparing the ten high-and-low-energy structures, it can generate a set of convex hull candidates, invoke AFLOW [32,57] to relax them, and use this data to generate composition-energy convex hull plots (Figure 1, green box).Finally, PRAPs contains a checkpoint system to allow users to start or re-start a job from a certain step.PRAPs is primarily a project management software that performs and automates many menial tasks, and generates plots that may be desired, thereby reducing the amount of human time required to generate MTPs for a particular system and analyze their performance.
In what follows we seek to answer two questions using PRAPs: (i) Can already existing quantum-mechanical data be used to train an MTP that is able to reproduce the DFTcalculated convex hull for a particular ternary alloy?(ii) And, can we use these MTPs to discover thermodynamically stable structures comprising the convex hull that are not present in the DFT-training set?In what follows, we show the answers to both questions is 'yes'.

B. Machine Learned Interatomic Potentials from AFLOW Data
To illustrate how already-existing DFT results can be scraped from large databases and repurposed towards the generation of system specific MLIPs, we chose four ternary alloy systems: CHfTa, CHfZr, CTaTi, and CMoW.TaC, HfC, ZrC, and TiC all adopt the rocksalt structure and are refractory ceramic materials with desirable mechanical properties [58].Though MoC can adopt the same rocksalt structure under pressure [59], at ambient conditions MoC and WC prefer a hexagonal arrangement instead [60].For nearly a century the propensity for transition metal carbides to form high-melting point solid solutions with compositions such as Hf 1−x Ta x C has been known [61,62].The metal formulation can be engineered to contain multiple atoms, and when there are five or more metals, the configurational entropy stabilizes single-phase high entropy carbides (HECs) [47,63,64] and their thin films [48].Recently, ML-IAPs have been developed for various HECs including a deep learning potential for (ZrHfTiNbTa)C 5 [65] and a low rank potential for (TiZrNbHfTa)C 5 [66].Herein, we illustrate that data found within AFLOW, coupled with active learning, can be used to train MTPs for ternary carbides, with the future goal of HECs in mind.
For training we employed DFT relaxation trajectories of ∼210 structures obtained from the AFLOW database, which have been generated through a combination of structure prototyping of naturally occurring compounds [67][68][69] and structure enumeration algorithms [70], followed by relaxation using density functional theory.This resulted in ∼5500-6000 in- dividual configurations on which the MTP was trained (Table I).From this training set, ∼800 configurations, with a minimum interatomic distance greater than 1.1 Å, were chosen randomly and MTPs of levels 10, 16 and 22 were trained in the "basic" mode.For the CMoW system, no level 22 data is reported due to the excessive computational cost required to obtain well-trained potentials.Five trainings were performed, and the best potential was chosen as the pre-Robust Potential (pre-RP).As expected, the training errors for the pre-RP (Table S1) decreased with increasing MTP level, with the average MAE (RMSE) on the energies being 27 (44), 16 (25) and 8 (13) meV/atom, and for the forces being 75 (217), 51 (147) and 26 (78) meV/ Å for levels 10, 16 and 22, respectively.Comparison of the ten highest and ten lowest energy structures as predicted by DFT and the pre-RP revealed that the MTP rarely miscategorized the structures, but had difficulty correctly categorizing more than 7 in the correct highest or lowest energy set (Table S2).
Previous studies have suggested that when used for CSP, MTPs should be trained on a very diverse set of structures [14,25], and as Table I shows, the original AFLOW data was somewhat limited.To obtain this diversity, we used the RANDSPG [56] program to create crystal lattices with up to eight atoms in the unit cell, for all possible ternary compositions.Lattice vectors were constrained to fall between 3-10 Å, and unit cell volumes between 200-600 Å3 , with a minimum interatomic distance of 1.1 Å.For each ternary system, ∼6300 structures were generated (Table I) and combined with the initial AFLOW data to form the relaxation set (Figure 1) used in training the RP.
For a given stoichiometry RANDSPG determines the compatible spacegroups, based upon the Wyckoff positions, and randomly chooses one prior to decorating its sites with atoms, thereby enabling the creation of random symmetric crystal lattices.Employing symmetric structures in the first generation of an evolutionary or particle-swarm directed CSP search greatly decreases the number of configurations that need to be optimized to locate the global minimum in the PES.This can be understood by noting that symmetric structures tend to be either very stable or unstable, spanning a greater amount of the potential energy hypersurface than those generated without symmetry constraints [5].Indeed, tests have shown that the average energy of random structures that are symmetric is higher (less negative) than those that are purely random [56].1).TABLE II.Prediction errors for energies (meV/atom) and forces (meV/ Å) by system, level, and potential.Each box contains the meanabsolute-error and, in parentheses, the root-mean-squared errors (RMSE).Only the Robust Potential (RP) and the Accurate Potential (AP) developed during the active learning scheme (ALS) as illustrated in Figure 1  Table II provides the prediction errors for the energies and forces for both the RP and the AP for different MTP levels.Determining the prediction error when such active learning schemes are used is complicated by the fact that each level of theory (DFT or MTP level) will encounter different structures during the relaxation process.While it is possible to predict the energy of every structure comprising the AFLOW relaxation trajectories via the RP or the AP, such a procedure is fraught with difficulties because of the noise in the DFT relaxation trajectory data.Specifically, the relaxation trajectory may have erroneous energies, forces and stresses resulting from a change in the planewave basis as the lattice vectors are varied during the course of the optimization.Moreover, the AFLOW data contained a few configurations whose magnitude of per-atom-energies were substantially larger than others, though they survived the distance-filtering-criteria (Supplementary Figure 1.Therefore, in Table II we have opted to compare the DFT energies and forces on the final AFLOW relaxed structures, which are accurately calculated, against their RP and AP predicted values. A number of trends can be observed from examination of Table II.First, the MAE are notably smaller than the RMSE, suggesting the presence of poorly-predicted outliers.The prediction errors with the RP generally decrease with increasing level of MTP, but the MAE do not fall below 40 meV/atom for energies and 60 meV/ Å for the forces.While these errors might seem large compared to the < 4 meV/atom and < 160 meV/ Å RMSEs computed for single-component systems wherein the testing and training dataset contained structures that could be derived via perturbations of the ground state crystal [26], they are in-line with some of the errors presented in Ref. [27] that trained MTPs for Li-Al alloys and applied them to a broad range of compositions and crystal lattices.In Ref. [25] the only errors reported for the Al-Ni-Ti system were for the training set with MAEs (RMSEs) of 18 (27) meV/atom for the RP and 7 (9) meV/atom for the AP.In Ref. [14] a RP yielded a training RMSE of 170 meV/atom for allotropes of boron, and an AP yielded errors of 11 meV/atom for the 100 lowest-energy structures that were found.
A key difference between the workflow presented here as compared to previous work by Gubaev and co-workers [25] is the pre-training step on the AFLOW data.To test what effect this may have, we performed the PRAPs procedure starting from an empty MTP, that is one that had not been trained on AFLOW data for a level 16 MTP.Therefore, the basic pretraining denoted schematically in the grey box in Figure 1 was not performed.Table II clearly illustrates that the MAEs for the energies were significantly improved by the AFLOW pre-training, and the RMSEs were slightly improved (by 94 and 16 meV/atom, on average).
A common practice in CSP is to use less accurate, but quicker methods to estimate the energies of a large number of structures, and then to optimize those with low energies with progressively more accurate, but costlier, methods.This workflow saves computational expense by filtering out structures that are unlikely to be stable prior to performing calculations that yield a more predictive rank order.A similar strategy was employed in Refs.[25] and [14] where both robust and accurate MTPs were trained for CSP.In addition to generating 375,000 binary and ternary bcc, fcc and hcp-type unit cells, 1463 Al-Ni-Ti ternary structures were created via decorating prototypes [25].It was noted that the prototypederived structures might possess geometries that were not chemically sensible because of short metal-metal distances and too-small volumes, which could lead to large MTP prediction errors.To circumvent this problem, structures with formation energies (as determined with a Robust MTP) that were within 100 meV/atom of the convex hull were chosen for rerelaxation using active learning starting from an empty MTP to generate an AP.In another related work, the large training error of 170 meV/atom obtained in a CSP search performed on boron, which has a very complex PES, prompted the authors to retrain the MTP on low-energy structures so it could accurately predict their energies [14].
Similar to these two studies [14,25], our RANDSPG generated structures span a wide energy range.Although the RPs obtained from including them in the training set have large errors, they are able to correctly identify most low-energy configurations, and not misidentify them as having high energies.Therefore, by choosing the RP-relaxed structures whose energies were within 50 meV/atom of the minimal energy structure for each composition, we can curate a new set of training data that can then be used to create an AP.Because the low-energy structures are expected to be better behaved than the full relaxation set, possessing nearly optimized geometries that are chemically sensible, the AP should yield lower training and prediction errors than the RP.
In many cases both the training and prediction RMS errors obtained for the AP decreased with increasing MTP level, however they did not fall below 50 meV/atom.Considering that most high-throughput and CSP searches are performed in a stepwise procedure, at each level filtering undesirable structures prior to increasing the accuracy by which the remaining subset is treated, this error is acceptable.Indeed, in both Refs.[25] and [14] all structures with sufficiently low energy, as predicted by the AP, were re-optimized with DFT to obtain an improved energy ordering.Similarly, in CSP searches on Lu-N-H [23] and Sc x H y [24] the lowest energy structures, as predicted by an ML-IAP, were relaxed with DFT to determine a more accurate energy rank.
For an MTP level of 22, the training in some cases took substantially longer than at lower levels.The reason for this is that the ALS, as originally designed [55], terminates when MLIP does not find any configurations that need to be added to the training set.Here a different protocol was used, where the active learning for the AP was stopped when the number of structures to be added to the training set was less than 1% of its size.This choice was motivated by the observation that at times fewer than ten configurations were being added to a training set of thousands in a single iteration, which took a full day to process, for a gain in RMS training error that was < 1 meV/atom.Tests showed that the choice of a variety of early termination criteria (which can be chosen as options in PRAPs) comes at the cost of 2-5 meV/atom in training error.
A key question to be answered is: "Does the ML-aided procedure reduce the total computational expense?"For each ternary the AFLOW and RANDSPG datasets contained ∼210 and ∼6300 individual structures, respectively (Table I).Therefore, ∼6500 geometry optimizations would be needed to relax all of these configurations.Using MTPs, a maximum number of single point DFT energy evaluations performed during the course of the training was ∼3000 (Figure 2) for the CTaTi system at level 22. Dividing the number of configurations comprising the AFLOW data by the number of structures gives an average estimate of the number of steps required per geometry optimization (∼28).Thus, our protocol reduces the number of DFT evaluations that would need to be performed to relax all of the considered structures by a factor of ∼30 or more (as high as 180 for the CHfTa system at an MTP level of 10).We note that as MTP level increases, the number of single-point calculations increases (as does the total training time), as expected (e.g.see Table 3 in Ref. [45]).The reason for this is that higher level MTPs contain more parameters, and therfore require more training data to fit these parameters.
Finally, we consider a few miscellaneous points about the setup and training procedures.Since the training data was predominantly for cubic structures, the HEAs which inspired our choice of elements typically crystallize in a cubic manner [47,48], and the ternary carbides considered herein assume either cubic or hexagonal lattices, it is fair to ask if the results could be improved if the RANDSPG structures were generated using only spacegroups with these symmetries.A test on all four systems, at an MTP level of 10, revealed that both training and prediction errors were very similar when using a reduced set of ∼200 RANDSPG structures across only HCP, FCC, and BCC-type unit cells in addition to the AFLOW data (Table S3).A reader may ask why the pre-training step is carried out here since active learning can be performed on an empty training set [14,25,45].Another test on CHfTa at level 10, showed that omitting the pre-training step resulted in a savings in time at the cost of training and prediction error (Table S4).The user can keep these trade-offs in mind when designing the computational protocol to be employed for their specific situation.

C. Predicted Convex Hulls
Let us now examine the convex hulls and investigate the structures that PRAPs relaxed with the robust and accurate potentials.Since optimization of random symmetric configurations is the first step in CSP, we examined if the aforementioned workflow could discover lattices not found within AFLOW whose energies lie on, or close to the convex hull.We Robust Accurate FIG. 2. Number of density functional theory (DFT) single point energy evaluations performed during the generation of the robust (top) and accurate (bottom) potentials according to the procedure illustrated in Figure 1 for each ternary carbide system at a particular MTP level (see legend).filtered out the AP predicted configurations that were within 50 meV/atom of the lowest energy structure for each composition and DFT-relaxed them via AFLOW.The resulting geometries were then concatenated with the fully relaxed AFLOW data to produce convex hulls, for different MTP levels (Figure 3 and Figures S9 and S10).
The data found within AFLOW for the CHfTa, CHfZr, and CTaTi systems all contained multiple compounds on or near the convex hull with the M 1−x N x C stoichiometry (x = 0.25, 0.33, 0.5, 0.67, 0.75) (Figures S3-S13).Examination of these structures suggested that they could all be obtained via relaxation of metal-carbide rocksalt structures whose metal sites were decorated with two different types of atoms, as would be expected for a solid-solution.On the other hand, for the CMoW system, AFLOW did not contain any ternary carbide phases with a Mo 1−x W x C composition that were on the hull (Figure 3).Unlike the cubic binary carbides comprised of metal atoms from group 5 or 6, isostructural MoC and WC (Figure 4(a)) adopt the hexagonal P 6m2 (#187) spacegroup, suggesting that Mo 1−x W x C stoichiometry structures would be hexagonal as well.Examination of a Imm2 symmetry Mo 0.5 W 0.5 C phase that was 17 meV/atom above the convex hull (teal dot in Figure 3) showed that it could not be derived from a decoration of a hexagonal metal-carbide lattice.Instead, it was related to a high-pressure phase of GaAs (AFLOW prototype AB oI4 44 a b) where the Ga atoms were replaced by C, half of the As atoms were substituted by Mo, and the other half by W.
In addition to the elemental endpoints, as well as MoC and WC, cubic MoW and rhombohedral C 2 Mo 4 (Figure 3) comprise the AFLOW convex hull for the CMoW system.Tetragonal Mo 14 W 2 (labelled by a purple dot) lies 1 meV/atom above the hull -a value that is within the error of the k-mesh and kinetic energy cutoffs employed in our planewave calculations.A different choice of DFT functional, inclusion of zero point energy or finite temperature contributions may place this structure on the hull.We then examined if MTPs trained by PRAPs on AFLOW data could relax structures created with RANDSPG and identify other thermodynamically stable compounds, not found within AFLOW for the CMoW system.
We compare the AFLOW hull with hulls predicted using the PRAPs procedure (Figure 3 and Figures S9 and S10).At an MTP level of 10 no additional structures emerged, and the PRAPsderived hull was identical to the one found within AFLOW.However, for an MTP level of 16, additional thermodynamically stable structures were found, including Imm2 Mo 0.75 W 0.25 C, in addition to the on-and-near-hull compounds present in AFLOW.The Imm2 phase, containing two formula units per primitive cell, resembles hexagonal MoC except that in every second layer half of the metal atoms are replaced by W, and the substituted metal-containing triangular nets are arranged in an • • •ABAB• • • stacking sequence with respect to each other (Figure 4(b)).This phase likely originated from the RANDSPG set, which was subsequently optimized, in an active learning sense, via the robust and accurate potentials.When an MTP of level 10 was used instead, Imm2 Mo 0.75 W 0.25 C was not on the convex hull likely because the relaxation process with the robust potential pushed this particular configuration too high in energy.To test this hypothesis the level 16 data for Imm2 Mo 0.75 W 0.25 C was concatenated with the structures that are present on the level 10 hull, and further analysis revealed that this phase was predicted to be thermodynamically stable.
In addition, four more structures, within 1 meV/atom of the convex hull, lay on the level 16 hull: P 6m2 Mo 0.5 W 0.5 C, P 6m2 Mo 0.333 W 0.666 C, P 6m2 Mo 0.666 W 0.333 C and Cm Mo 0.25 W 0.75 C (Figure 4(c-f)).Though the first has the same composition as the structure present within AFLOW, it is 17 meV/atom lower in energy.In fact, if we do not distinguish between the identities of the metal atoms, the AFLOW structure can be transformed into P 6m2 Mo 0.5 W 0.5 C by doubling it along the b-axis followed by three sets of translations of various subsets of atoms.In both phases the C atoms fall within trigonal-prismatic holes, but in this particular structure the triangular (and square) faces all point along the same crystallographic direction, while in the AFLOW structure half of the prisms are rotated, thereby swapping the axes along which the two sets of faces lie.Importantly, PRAPsfound P 6m2 Mo 0.5 W 0.5 C corresponds to a coloring of the hexagonal CMo/CW prototype structure with an •••ABAB••• arrangement for the metal-containing hexagonal nets (Figure 4(c)).Similarly, the remaining three PRAPsfound structures can be derived from colorings of the hexagonal parent phase, with P 6m2 Mo 0.333 W 0.666 C and P 6m2 Mo 0.666 W 0.333 C being inverses of each other, while Cm Mo 0.25 W 0.75 C can be described as a W-rich same hexagonal prototype.
The identified near-and-on-hull phases lie on a straight line joining the two end-members comprising this CMo/CW series.They represent examples of an ensemble of phases with highly-variable concentrations, suggesting the existence of a solid solution with a very low critical temperature of the miscibility gap.To investigate this, we optimized ∼366 Mo 1−x W x C structures (x = 0.08 3, 0. 1, 0.125, 0.1 6, 0. 2, 0.25, 0. 3, 0.375, 0.41 6, 0. 4, 0.5, 0. 5, 0.58 3, 0.625, 0. 6, 0.75, 0. 7, 0.8 3) with 4-24 atoms in the unit cell, and between 2 and 86 unique structures were optimized per composition.The previously generated level 16 robust and accurate potentials were used to predict their energies and to relax them.Figure 4(g) plots the resulting enthalpies of formation, ∆H, from the monocarbide endpoints: relaxed with the robust potential (RR), subsequently predicted by the accurate potential (AP-RR), and finally relaxed with the accurate potential (AR-RR).
All of the DFT-optimized compounds fell on or within 5.7 meV/atom of the line joining the CMo and CW end points, suggesting that their ∆H is close to 0 meV/atom.For a given composition, various decorations were computed to be nearly isoenthalpic, suggesting that configurational entropy will play a role in the stability of this family of structures.Turning to the results obtained with the generated MTPs, the computed ∆H, as predicted on structures relaxed by the RP, was largely positive (blue dots) with the deviation from the zero-energy line steadily increasing for larger W concentrations.Whereas the distance from the CMo-CW tie-line, averaged over all structures, was calculated as being 0.83 meV/atom (σ = 1.10) via DFT, the robust relaxed protocol resulted in an average tie-line distance of 78.39 meV/atom (σ = 28.95).Prediction of the energies of the robust relaxed structures with the AP (green dots) yielded an average ∆H of 18.82 meV/atom (σ = 14.96).It is only via relaxation with the AP (purple dots) that we obtain an average tie-line distance of 4.20 meV/atom (σ = 3.58).This example illustrates that structural relaxation with the AP is key for obtaining energetics that are in good agreement with those derived from DFT calculations.
The convex hulls discussed and presented above (Figure 3) were optimized with DFT, and the conclusions regarding thermodynamic stability of particular phases were made based upon these hulls.This procedure is common-place [71] but it might make one wonder about the limits of the utility of ML-IAPs in CSP.Part of the answer lies above where we show that ML can significantly reduce the number of required DFT calculations.But, the other part of this answer is in the convex hull candidate structures: the output of PRAPs relaxations and predictions before the final DFT step.The analysis of the CMoW system suggested that relaxation with the AP is key for obtaining energetics that are in-line with DFT results.To further study this aspect, in Figure 5 we plot the convex hulls for the CHfTa system calculated at an MTP level of 22.Comparison of the AFLOW derived hull with one that is obtained after relaxation with the robust potential (RR) shows that the latter predicts a structure that is not found within AFLOW, with Hf 0.5 TaC 0.5 composition, to lie on the hull (after DFT relaxation, it falls 123 meV/atom above the hull) whereas Hf 1−x Ta x C stoichiometries lie around 15 meV/atom above the hull.The rogue Hf 0.5 TaC 0.5 structure disappears after AP prediction, and the energies of the Hf 1−x Ta x C species fall onto-and-just-above the hull.Relaxation with the AP yields a hull that is virtually indistinguishable from the one derived from AFLOW, similar to the results obtained for the CMoW system.In the Supplementary Information we provide these same four convex hulls for each carbide system considered and each MTP level, before and after subsequent relaxation with DFT.Generally speaking, the structures that fall within 60 meV/atom of the convex hull, after being relaxed with the RP followed by the AP, also fall within 60 meV/atom of to finish These examples suggest that relaxations performed with an AP will be useful as a further screening step that can be undertaken prior to performing expensive DFT calculations in the materials prediction workflow.

D. Discussion
The density functional theory (DFT) computed energies, forces and stresses found within the AFLOW database of four ternary carbide systems (HfTaC, HfZrC, MoWC and TaTiC) were employed to train system specific machine learning interatomic potentials of the moment tensor potential (MTP) flavor.A utility package that can be used to generate both robust potentials (RP), capable of roughly relaxing any structure, and accurate potentials (AP), tailored towards the relaxation of low-energy structures, which was employed to automate this training, is described.The AFLOW data was augmented with ∼6300 random symmetric structures resembling those that would be created in the first step of a crystal structure prediction (CSP) search, and these were relaxed with MTPs updated via active learning.For the HfTaC system, relaxation with the AP yielded a convex hull that agreed perfectly with the one found within AFLOW.Moreover, this procedure identified five Mo 1−x W x C stoichiometry compounds, not found within AFLOW, that lay on the convex hull and corresponded to colorings of the hexagonal CMo/CW prototypes, illustrating how the described protocol can accelerate CSP.Subsequently, the RP and AP were used to relax hundreds of Mo 1−x W x C lattices spanning a broad composition range, and it was shown that relaxation with the AP yielded formation enthalpies that were in excellent agreement with those computed via DFT, and the resulting convex hull exhibits regions in which entropy plays a considerable role on the phase stability.The ideas and tools described here may aid in the generation of ML-IAPs from already existing DFT data, to be used for materials prediction.

METHODS E. Computational Details
The density functional theory (DFT) calculations were performed using the Vienna ab initio Simulation Package version 5.4.12 [72] coupled with the Perdew, Burke, Ernzerhof (PBE) gradient-corrected exchange and correlation functional [73] and the projector augmented wave method [74].During the active learning procedure, the VASP calculations were performed using Γ-centered Monkhorst-Pack k-meshes where the number of divisions along each reciprocal lattice vector was chosen such that the product of this number with the real lattice constant was 30 Å.The carbon 2s 2 2p 2 , Hf 6s 2 5d 2 , Ta 6s 2 5d 3 , Zr 5s 2 4d 2 , Ti 4s 1 3d 3 , Mo 5s 2 4d 4 and W 6s 2 5d 4 electrons were treated as valence, and an energy cutoff of 400 eV was employed.After training was complete, the convex hull analysis included a DFT relaxation accomplished by calling AFLOW's management protocol, using the standard settings described in Ref. [57]; the AFLOW hull data was also re-relaxed using this procedure.Structures from the AFLOW database and those generated by RANDSPG [56], as described in the main text, comprised the full relaxation set employed for the development of the MTPs.The crystals whose geometries were relaxed to construct Figure 4 (g) were generated from P 6m2 Mo 0.5 W 0.5 C using the Supercell program, em-

AFLOW RR AR-RR AP-RR
FIG. 5. Structures within 60 meV/atom of the AFLOW derived convex hull (top left), as well as the PRAPs procedure at an MTP level of 22 for the CHfTa system.Structures are colored (see color bar) according to the distances from the hull obtained after relaxing with the robust potential (RR), prediction of the enthalpies of the RR structures with the accurate potential (AP-RR) and relaxing the RR structures with the accurate potential (AR-RR).Black dots are on the hull, and purple dots are within 1 meV/atom of the hull.In contrast to the plots shown in Figure 3, the data illustrated here has not undergone further DFT relaxations.
ploying the merge option to remove duplicate structures [75].
PRAPs was run on each ternary carbide using MTP levels 10, 16, and 22 with a MLIP-relaxation-iteration limit of 100, and an extrapolation grade of 2 < γ ≤ 10.The cutoff distances for generating the MTP were 1.1 Å < x < 5 Å.Active learning was, in most cases, declared to be converged when no additional structures were considered for addition to the training set.In the case of the level 22 trainings, the ALS procedure was stopped when the number of structures to be added to the training set was less than 1% of the number already in the training set.PRAPs filtered out configurations with interatomic distances 1.1 Å < d < 3.1 Å before beginning the pre-training, and when beginning the AP training removed all structures that were higher than 50 meV/atom of the most stable configuration for each composition.
Code Availability.The PRAPs code will be released in a subsequent publication, and in the meanwhile, may be obtained from the corresponding authors upon reasonable request.

Data Availability
The datasets generated during and/or analyzed during the current study are summarized in the supplementary information, and are available from the corresponding author on reasonable request.

S3 Convex Hulls
This section contains a subset of the convex hulls calculated by PRAPs.If desired by the user, PRAPs will calculate eleven convex hulls: two for the DFT data before-and-after relaxation, and nine using the Robust and Accurate Potentials.Six of the hulls are presented here in a table-like format, organized into rows and columns from Figures S3 to S13

FIG. 1 .
FIG.1.Workflow in the Plan for Robust and Accurate Potentials (PRAPs) package, which automates the generation of a moment tensor potential (MTP), given any quantum mechanical training set (here the online AFLOW database[32]).Five MTPs are trained on a set of configurations (cfgs) and the best (the pre-Robust Potential, pre-RP) is chosen (gray box).The pre-RP is employed to initialise the training of a RP via an active learning scheme (ALS, blue box).The RP is generated by relaxing the configurations comprising the Relaxation Set.The RP-relaxed lowest energy configurations are chosen (filtration step) to train an Accurate Potential (AP) via active learning (orange box).The AP-relaxed structures may be sent for relaxation via AFLOW-DFT, followed by subsequent analysis (green box).The inset schematically illustrates the energy distribution of the structures that can be relaxed with the RP and the AP.Optionally, structures generated using such a procedure may not be exact since the intermediate structures in a DFT relaxation trajectory may have erroneous energies or forces resulting from Pulay stress that arises when the lattice vectors are varied.another procedure (here RANDSPG[56]), labelled as 'Non-DFT Data', may be relaxed via the RP and/or the AP, and this auxillary set of structures also comprises the Relaxation Set.Rectangles with sharp (curved) corners represent structural data (potentials).
are shown here; data for the pre-Robust Potential can be found in the SI.The training errors for each potential reflect the RMS values obtained by comparing the energies of that potential's training set as obtained by single-point DFT calculations against those predicted by the respective potential.The prediction errors for the RP indicate comparison of the energies of the original AFLOW data with the energies predicted by the respective potential.The prediction errors for the AP indicate comparison of only the low-energy AFLOW structures with those predicted by the respective potential.

FIG. 3 .
FIG. 3. Convex hulls obtained by DFT-optimizing structures predicted to be within 50 meV/atom of the convex hull obtained from AFLOW (left) and the PRAPs procedure at an MTP level of 16 (right) for the CMoW system.The black dots are structures lying on the convex hull, and phases within 60 meV/atom of the hull are represented by color dots (see color bar).Purple dots are within 1 meV/atom of the hull.

TABLE I .
Number of unique structures for the studied ternary carbides found within the AFLOW database, and the number of individual ionic steps (configurations) from their structural relaxations.The number of individual structures generated by RANDSPG that were employed to develop the robust and accurate potentials using the MLIP program package are also provided.
Calculations performed without pre-training the MTP on the AFLOW data (grey box in Figure *