Efficient prediction of elastic properties of Ti 0.5 Al 0.5 N at elevated temperature using machine learning interatomic potential

.


Introduction
Hard ceramic coatings are utilized to extend the wear resistance and lifetime of machining tools in the cutting zone subjected to excessive heat and mechanical load.In general, linear elasticity at operational temperature is essential to determine forces acting on dislocations in materials and understand plastic deformations [1].It means that the knowledge of elastic constants at elevated temperature supports the design of high-strength, wear resistant industrial materials.On the other hand, phase field simulations have shown the impact of elastic energy on microstructure evolution during spinodal decomposition of alloys [2,3] and Zener's elastic anisotropy has contributed to the understanding of the age hardening in Ti (1-x) Al x N [4][5][6].
Accessing properties of hard coatings at operational temperature, which is up to around 900-1000 • C [7,8] for TiN based materials, is a challenging task for computer simulations based on the laws of quantum mechanics.Besides the freedoms of the electronic motions described through density functional theory calculations the atomic vibrations have to be simulated and the corresponding vibrational entropy should be calculated.Sophisticated approaches utilizing ab initio molecular dynamics (AIMD) simulations have recently been introduced to include the contribution of atomic vibrations to the material's free energy at arbitrary temperatures.Self-consistent phonon theory [9,10], stochastic self-consistent harmonic approximation [11] or the temperature dependent effective potential method (TDEP) [12] are among the approaches used to calculate effective harmonic force constants to predict, for example, lattice thermal conductivity in cubic SrTiO 3 [13], superconducting critical temperature of hydrogen sulfide [14] or phase stability of CrN [15].Recently, the phase diagram of Ti (1-x) Al x N random alloys has been calculated using TDEP [16].The approach has also been extended to predict the elastic tensor of Ti (1-x) Al x N random alloys at elevated temperatures [17].
A bottleneck of all these calculations is the high computational demand of the underlying AIMD simulations, especially in case of metallic alloys, where a large supercell model is needed in combination with dense k-point sampling in the Brillouin zone.This issue gets more severe in case of predicting elastic constants because of the need for accurate stresses.Density functional theory calculation with approximate exchange-correlation functional typically results in 10% error in comparison with experiment on absolute values, but it is much more accurate in predicting trends e.g. on the temperature dependences.In accelerating computational schemes one does not want to add other errors on top of the density functional theory and the highly accurate (within 1 GPa) predictions of the temperature dependences of elastic constants are realistic.One needs ∼1 meV/atom accuracy in energy differences to predict elastic constants with ∼1 GPa accuracy.To moderate the computational demand, approximate AIMD approaches have been introduced; the quasi-static or volume scaled simulations (VSS) [18] and the symmetry imposed force constants (SIFC) approximation combined with the temperature dependent effective potential (SIFC-T-DEP) approach [17].
An alternative to AIMD is to use classical molecular dynamics simulations in which forces and stresses are computed from an empirical interatomic potential (IP).Since classical molecular dynamics does not account for the electronic degrees of freedom the calculation accuracy depends on the applied IP.The classical embedded atom model [19] and modified embedded atom model interatomic potentials [20] can be fitted to carefully performed ab initio calculations to provide quantitative accuracy, however, their applicability for different properties or various materials, phases or compositions is limited by the prescribed functional form of the IP and by the implemented parameter values.
Learning-on-the-fly strategy [21,22] has been suggested to incorporate static first-principles calculations into classical IP models and achieve molecular dynamics trajectories on a quantum mechanical level of accuracy.However, to succeed one needs not only a flexible and complete many-body functional form of the IP but also an effective machine learning strategy.The recently developed moment tensor potential (MTP) approach [23] uses contracted moment tensors as structural descriptors in generating IPs, called MTPs.Furthermore, an active learning method based on the maxvol algorithm [24] has been suggested to train MTPs on-the-fly [25].Machine learned MTPs have successfully been utilized, for example, in crystal structure search [26], predicting thermal conductivity of complex compounds [27] or modeling self-diffusion in aluminum [28].
In this study, as proof of concept, we show that MTPs can be trained on AIMD simulation results to obtain highly accurate elastic constants of a prototypical hard coating alloy, Ti 0.5 Al 0.5 N at 900 K.We note that an active learning approach of MLIP fitting [26] has been applied to accelerate AIMD simulations for elemental metals [29], and for multicomponent alloys [30] reduces the computational effort by two orders of magnitude.A comparison of calculated elastic properties using machine learned MTP with values calculated by highly accurate AIMD, as well as by approximate AIMD-based methods [17,31] shows that with the machine learned MTPs one achieves both, high efficiency and excellent accuracy of the simulated results.

Methodology
Ti 0.5 Al 0.5 N solid solution has been modeled by a (4 ×4 ×4) supercell of the B1 unit cell with totally random arrangement of Ti and Al atoms on the metal sublattice.High accuracy AIMD canonical (NVT) simulations have been performed to characterize the material's properties at four different temperatures from 300 K up to 1500 K.In this study we utilize the results of the simulations performed at 900 K.The calculations of the elastic moduli of Ti 0.5 Al 0.5 N have been carried out using finite distortions with the distortion matrix with η = ±0.01,±0.02 strains applied at the calculated equilibrium lattice parameters.The thermal expansion of Ti 0.5 Al 0.5 N was extracted from Ref. [17].At each η around 15,000 AIMD steps with 0.5 fs time step have been performed.The non-equivalent averaged cubic elastic constants C 11 , C 12 and C 44 have been calculated by performing simulations with three orientation of the coordinate system, (xyz), (yzx) and (zxy).A comprehensive discussion of the AIMD simulations and the derived results is out of the scope of the present study and will be presented elsewhere [32].
The calculation scheme based on the use of the machine learned MTPs, in general, has three major blocks, as shown in Fig. 1.In the first block an arbitrary method should be utilized to generate a set of atomic arrangements (dataset) at the given thermodynamic conditions with corresponding total energies, stresses and forces (labels) calculated with density functional theory.These labels are used with different weights to fit the MTP.We call this labeled dataset as full set of selected labeled data (FS).In our case, only a selected set of configurations from the above discussed AIMD simulations have been used to create FS.In general, three types of labeled datasets have been created from the total ≈60000 labeled data of the full AIMD simulations: FS, a subset of FS,  which is sampled to create training sets (TSs), and a validation set (VS) to validate the trained MTP.FS has been created as follows.Our aim has been to find a reasonably small but still representative TS, therefore, we have utilized the AIMD simulations performed with ±2% strains only, which have served summing up the three lattice orientations altogether 30000 labeled data.This choice is supported by the argument that TS should provide a representation with sufficient volume variation.The selected six AIMD runs have always been treated independently.We have cut away the first 1000 from each simulations corresponding to the initial equilibration processes.The remaining 24,000 configurations with the corresponding density functional theory calculated total energies, stresses and forces have been combined to define FS.
In the second block of Fig. 1 called MLIP calcs.we have created VS and the training sets (TSs).First a sampling set (SS) has been created by combining the first 2500 configurations from each simulations.To create VS, which is uncorrelated from SS we have skipped 1000 configurations after the selected 2500 configurations and combined the remaining around 500 configurations from the six simulations.Thus, SS has contained around 15,000 configurations while VS around 3000.Training set TS 1 with 300 configurations has been created by randomly selecting 2% of each SS, such that equaly 2% has been selected from each simulations.We have made three independent random sampling resulting in TS 1 /A, TS 1 /B and TS1 1 /C.TS 2 training sets have been formed by randomly selecting 4% of each SS.The TS 2 sets have contained 600 at.configurations.
Various levels of approximations have been used in defining the functional form of IP.Depending on the highest degree of polynomiallike basis functions used in the analytic description of the MTP [23], we have used the so called 10g, 16g, 18g, 20g and 24g MTPs.IPs have been derived by fitting the MTPs to TSs using the Broyden-Fletcher-Goldfarb-Shanno method [24] with 2000 iterations and 1.0, 0.01 and 0.01 weights for total energy, stresses and forces in the loss functional.
Elastic constants (C ij ) of Ti 0.5 Al 0.5 N have been calculated in the third block of Fig. 1 called LAMMPS cals.using averaged stresses and strains extracted from the distortions given by the distortion matrix in Eq. ( 1) with ±1, 2% strain.The elastic tensor has been calculated by employing the symmetrization technique [33] to account for the disorder present in the supercell.This means that we have averaged C 11 , C 12 and C 44 elastic constants calculated with three different orientiations of the same supercell labeled with ϕ [111] = {0 ∘ , 120 ∘ , 240 ∘ }.The classical molecular dynamics simulations with the trained IPs were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [34] with 1 fs time step and performing 100.000 steps.Note that Ti 0.5 Al 0.5 N has been modeled with exactly the same supercell as in the AIMD NVT-simulations.

Results
Table 1 gives an overview of the accuracy of the MTP fitting.By comparing the value of root mean square (RMS) error of the fitting with the validation RMS error one can state that in each cases the choice of TS is representative.Furthermore, the small variation of RMS with changing TS 1 /A to B and C rules out overfitting.To predict elastic constants  2 for all values.

Table 2
Elastic constanst of Ti 0.5 Al 0.5 N at 900 K using TS 1 and TS 2 training sets created from AIMD simualations and fitted to various levels of MTP potentials.The number of parameters in MTPs are given in parenthesis.The supercell size in Lammps calculation is 128 atoms.Sym.stands for symmetry averaged elastic constants.* Ref [17] In the following we compare elastic constants calculated by different methods.The values are listed in Table 2. C 11 =362 GPa, C 12 =145 GPa and C 44 =187 GPa are the cubic averaged elastic constants obtained by direct highly accurate AIMD simulations [32].These values are used as reference hereafter.The table also lists the values obtained by applying the SIFC-TDEP [17] and VSS [31] approach.
To calculate elastic constants we have selected the model from Table 1, which results in the smallest RMS error in validation.For example, in the choice of 10g MTP and TS 1 the random selection C, that is the 10gTS 1 /C model results in the smallest RMS.However, in some cases the LAMMPS simulations used to calculate Cij have led to unphysical results for the selected model, e.g.(i) to unphysically large pressure or (ii) to unrealistically large value of the bulk modulus.In such cases we have selected the model with next smallest RMS error in validation, followed by the nex one if the unphysical features remain.For example, in case of 18g and TS 2 the model with second lowest RMS error (18gTS 2 /A) has resulted in elastic constants.The obtained elastic constants are listed in Table 2.
A comparison of the relative errors of C ij s obtained with different methods is shown in Fig. 2. One sees that from 16g MTPs one obtaines all C ij s with less than 2% error wrt. the highest accurate AIMD values.Zener's type elastic anisotropy values have been calculated and compared in Fig. 3.The range of 10% error from the AIMD calculated value is marked with two dashed lines.One sees that all the MTP potentials result in less than 10% error.However, SIFC-TDEP results in 25% error while VSS in 15%.One sees that from and above the order of 16g MTP the error in predicted Zener's anisotropy is below 4%.The 20g MTP results in less than 1.5% error.In summary, using 20g MTP one achieves an accurate elastic chracterization of Ti 0.5 Al 0.5 N at 900 K with less than 2% relative error with respect to highly accurate AIMD simulations.

Conclusions
Elastic constats of Ti 0.5 Al 0.5 N at 900 K have been calculated using classical molecular dynamics simulations with machine learned moment tensor potentials.Using 16 as the highest degree of polynomial-like basis functions in the analytic description of the moment tensor potential (16g) the accuracy of the predicted elastic constants is within the ± 2% of the values derived from the highly accurate direct ab initio molecular dynamics simulations.In comparison with recently proposed approximate ab initio molecular dynamics based acceleration schemes, which result in 24% error for the Zener elastic anisotropy, 16g moment tensor potentials result in less than 4% difference.Though, in this study the moment tensor potentials have been created using data of ab initio molecular dynamics simulations, we foresee that the potentials can be efficiently trained using an active learning leading to accurate values of elastic constants and greatly improving the efficiency of simulations of mechanical properties of hard coatings at operational conditions of cutting tools.

Fig. 1 .
Fig. 1.Flow diagram of the calculation process.See text for the detailed explanation.

Fig. 3 .
Fig. 3. Zener elastic anisotropy values of Ti 0.5 Al 0.5 N calculated by ab initio molecular dynamics (AIMD) simulations without additional approximations, symmetry imposed force constants approximation combined with the temperature dependent effective potential (SIFC-TDEP),the quasi-static or volume scaled simulations (VSS) approximation and different MTPs.

Table 1
Total energy errors in eV unit for the fitting Err.(TS) and validation Err.(VS) using a supercell with 128 atoms.The number of parameters in MTP is given in parenthesis, which are optimized using 2000 iterations with energy weight 1.0, force weight 0.01 and stress weight 0.01.TS 1 training sets contain 300 configurations, while TS 2 600 ones.Validation error is checked on 3362 configurations.MAD means maximum absolute difference and RMS stands for root mean square.