Semi-Empirical Force-Field Model for the Ti1−xAlxN (0 ≤ x ≤ 1) System

We present a modified embedded atom method (MEAM) semi-empirical force-field model for the Ti1−xAlxN (0 ≤ x ≤ 1) alloy system. The MEAM parameters, determined via an adaptive simulated-annealing (ASA) minimization scheme, optimize the model’s predictions with respect to 0 K equilibrium volumes, elastic constants, cohesive energies, enthalpies of mixing, and point-defect formation energies, for a set of ≈40 elemental, binary, and ternary Ti-Al-N structures and configurations. Subsequently, the reliability of the model is thoroughly verified against known finite-temperature thermodynamic and kinetic properties of key binary Ti-N and Al-N phases, as well as properties of Ti1−xAlxN (0 < x < 1) alloys. The successful outcome of the validation underscores the transferability of our model, opening the way for large-scale molecular dynamics simulations of, e.g., phase evolution, interfacial processes, and mechanical response in Ti-Al-N-based alloys, superlattices, and nanostructures.


Introduction
Ti 1−x Al x N (0 ≤ x ≤ 1) alloys are widely used for wear and oxidation protection of metal cutting and forming tools, aerospace components, and automotive parts [1,2]. The parent binary phases B1-TiN (space group Fm-3m) and B4-AlN (space group P6 3 mc) are immiscible at room temperature and at thermodynamic equilibrium [3]. However, far-from-equilibrium conditions, which prevail during vapor-based thin-film deposition [4], enable synthesis of metastable B1-Ti 1−x Al x N solid solutions for x ≤ 0.7 [5]. These metastable alloys undergo, in the temperature range ≈1000 to ≈1200 K, spinodal decomposition into strained isostructural Al-rich and Ti-rich B1-Ti 1−x Al x N domains, which in turn, results in an age-hardened material with superior high-temperature mechanical and oxidation performance [6].
Over the past two decades, numerous theoretical and experimental studies [6][7][8][9][10][11][12][13] have focused on the mechanical properties and chemical stability, kinetics and energetics of structure formation, and phase stability of Ti 1−x Al x N (0 ≤ x ≤ 1) alloys, during vapor-based thin-film synthesis and high-temperature operation. Despite the knowledge generated from these investigations, atomic-scale processes that drive phase transformations (including spinodal decomposition) and elastic/plastic response in B1-Ti 1−x Al x N alloys (0 ≤ x ≤ 1) are poorly understood. This is because the time and length scales at which these processes occur often lie beyond the temporal and spatial resolution An effective and complete description of all pair and triplet interactions within the Ti-Al-N system requires the selection of values for ≈10 2 MEAM parameters (see Table S1 in Supplementary Materials). We address this formidable challenge by using a stochastic algorithm based on an adaptive simulated annealing (ASA) scheme [37], which is explained in detail in Section S2 and illustrated in Figure S3 in Supplementary Materials. This algorithm allows us to determine parameter values that optimize the potential's predictions relative to experimentally-and theoretically-determined physical properties of phases and configurations of Al, Ti, and N and their binary and ternary combinations, as explained later in the present section. This is a necessary, yet not sufficient, condition for developing a MEAM potential which is fully transferable among different compositions with varying x values (0 ≤ x ≤ 1) and metal-sublattice atomic arrangements within the Ti 1−x Al x N alloy system. To reduce the size of the multi-dimensional parameter-space that needs to be scanned during the ASA optimization procedure, we constrained the range of key MEAM parameters, i.e., interatomic distance r e , cohesive energy E c , bulk modulus B (E c , B, and r e determine the MEAM parameter α), and attractive d + and repulsive dinteraction, to be close to values that best fit density functional theory (DFT) total energies versus equilibrium volumes (i.e., equation of states), for all reference elemental and binary phases. More details and results on the above-mentioned DFT calculations can be found in Section S1 and in Figures S1 and S2 in Supplementary Materials.
For the phases listed above, we used ASA to determine sets of MEAM parameters that optimize the potential's predictions for the following 0 K physical properties: (i) cohesive energies (E c ); (ii) lattice parameters (a, c), and elastic constants (bulk moduli B and C ij ) for all phases; (iii) mixing enthalpies (∆H mix ) for B1-Ti 1−x Al x N random solid solutions calculated with respect to the energies of B1-TiN and B1-AlN parent binary phases (∆H mix = E c,B1-Ti1−xAlxN − (1-x)·E c,B1-TiN − x·E c,B1-AlN ); and (iv) point-defect (vacancies V and interstitials I) formation energies E f V/I for B1-TiN and B1-AlN. Bulk moduli and lattice parameters were determined using least-square fitting of the total energy versus equilibrium volume curves, while the elastic constants C ij were calculated using least-square fitting of energy versus strain curves for structures at their equilibrium volumes. E f V/I values were computed with respect to the chemical potential µ of Ti, Al, and N in hcp-Ti, fcc-Al, and N 2 molecules, respectively. Specifically, E f V/I = E D ± µ -E 0 , where E D is the energy of a crystal containing one point-defect, E 0 is the energy of the defect-free crystal, while the sign preceding µ is positive for vacancies and negative for interstitials. All physical properties mentioned above were determined via MEAM conjugate-gradient energy minimization calculations at 0 K using lattices consisting of 250 to 450 atoms. Exception to the latter are the B1-Ti 1−x Al x N systems, where all properties were evaluated by averaging over the results obtained for 10 different random cation-sublattice configurations using supercells formed by 1000 atoms. For reference, we also calculated N and Ti vacancy formation energies in B1-TiN with respect to the chemical potential of isolated Ti and N atoms by means of both MEAM and DFT. The latter was done using the VASP code [38], the Perdew-Burke-Ernzerhof (PBE) approximation [39] for the exchange and correlation energy, the projector augmented wave (PAW) method to describe electron-core interactions [40], and supercells formed of 215 atoms + 1 vacancy. The total energy of the relaxed TiN system was evaluated to an accuracy of 10 -5 eV/supercell, using 3 × 3 × 3 k-point grids and 400 eV cutoff-energy for planewaves. The energy of isolated N and Ti atoms was computed accounting for spin relaxation.
The quality of each MEAM parameter set was further assessed (during the ASA optimization) by evaluating the potential's performance with respect to structural stability, melting points (T m ), and phase transitions induced in AlN upon changing temperature T and/or pressure P. This was accomplished by means of 5 to 30 ps long CMD simulations at finite T and P values for cell sizes of ≈1000 atoms. Structural stability was tested by verifying that the simulated systems exhibit no unphysical behavior during 30 ps and that total energies are conserved during microcanonical NVE sampling at temperatures of ≈0.5 T m . Melting points were rapidly estimated as the temperature values for which the derivative of equilibrium volumes versus T becomes sharply discontinuous by performing 5 ps CMD isobaric-isothermal NPT simulations at different temperatures near experimental T m values. The relative stability of competing phases was assessed by calculating the free energies of such phases on a grid of T and/or P values via thermodynamic integration (TI) [41]. The reference energy for TI is the free energy of an Einstein crystal, i.e., a system of non-interacting harmonic oscillators with frequencies determined from the mean square displacements (MSD) of the material system under consideration [42]. Using the assessment process described above, along with the ASA optimization scheme, we defined the best set of MEAM parameters, which is validated as explained in Section 2.2.

Potential Validation Methodology
In order to explore the transferability of the parameters determined using the procedure described in Section 2.1, we studied the following properties and dynamic processes that are not explicitly included in the ASA optimization scheme: (i) lattice thermal expansion of B4-AlN and B1-Ti 1−x Al x N (0 ≤ x ≤ 1); (ii) finite-temperature phonon spectra of B1-TiN, B1-AlN, and B4-AlN; (iii) nitrogen-and metal-vacancy migration energies for B1-TiN, B1-AlN, and B4-AlN, as well as diffusion energetics of Al interstitials in B1-TiN; (iv) lattice and elastic constants of sub-stoichiometric B1-TiN y (0.7 ≤ y ≤ 1), B2-TiN (space group Pm-3m), body-centered tetragonal (bct) Ti 2 N (space group I4 1 /amd); (v) free energies of B1-AlN and B4-AlN and dynamics of B4-to B1-AlN phase transformation; and (vi) B1-Ti 1−x Al x N alloy mixing free energies with respect to B1-TiN and B1-AlN.
The first step for studying the effect of temperature on the lattice parameter is to determine the temperature-dependent equilibrium volumes of different phases. This was done by performing CMD simulations employing NPT sampling of the phase space, using the Langevin thermostat coupled to the Parrinello-Rahman barostat [43], with settings suggested in Ref. [41]. The supercells used for these simulations are typically formed of ≈500 atoms. Equilibrium volumes and lattice parameters versus T (0 < T ≤ 3300 K) curves were then fitted to second-order polynomials to obtain the lattice thermal expansion coefficient.
Finite-temperature phonon spectra and free energies (including vibrational entropies) were obtained using post-simulation processing of CMD forces and atomic displacements via the temperature-dependent effective-potential (TDEP) method [44][45][46], using both second-order and third-order force constants. The TDEP method, which allows for the evaluation of the dynamic matrix directly at a temperature of interest, has been successfully applied for understanding the vibrational effects on the dynamic stabilization of crystal lattices [47] and solid-solid phase transitions [48]. For TDEP simulations in this work, we used canonical NVT sampling at equilibrium volumes previously determined via NPT sampling at a given temperature and pressure, and we controlled the temperature via the Nose-Hoover thermostat. Moreover, we compared CMD phonon-spectra to those obtained via TDEP applied to density-functional ab initio molecular dynamics (AIMD) [49] simulation results. AIMD/TDEP data are available in the literature for B1-AlN and B4-AlN [50]. For B1-TiN, we carried out additional AIMD simulations (VASP, PAW, PBE, 300 eV cutoff energy, and Γ-point sampling) of 10-ps duration on supercells formed of 512 atoms at temperatures of 300 and 1200 K. For calculating the free energy of mixing for B1-Ti 1−x Al x N random solid solutions, we imposed a symmetric force constant matrix in the framework of TDEP. This was implemented by using reference ideal B1 primitive cells with the metal sublattice occupied using a single elemental species Q of mass M Q = (1-x)·M Ti + x·M Al , as proposed in Ref. [51]. The configurational entropy S = -k B [x·ln(x) + (1-x)·ln(1-x)] per formula unit of B1-Ti 1−x Al x N alloys was evaluated via the mean-field theory approximation. Binodal and spinodal compositions x were determined as a function of temperature. Binodal points, which separate stable from metastable alloy regions in the composition space, were obtained from the tangents to the B1-Ti 1−x Al x N mixing free energy curves [52]. Spinodal points, which separate metastable from unstable B1-Ti 1−x Al x N regions, were determined from the change in sign in the second derivative of alloy mixing free energies versus x [52].
Vacancy and interstitial migration energies, as well as minimum-energy paths, were obtained via static MEAM nudged-elastic band (NEB) [53,54] calculations. The elastic constants of B1-TiN y (0.7 ≤ y ≤ 1) were determined by averaging over the values calculated for ≈100 different configurations (supercells containing 512 lattice sites) using the scheme described in Ref. [55]. The dynamics of B4-to B1-AlN transformation was studied via NPT CMD simulations at 300 K using a timestep of 0.1 fs. First, a supercell containing ≈20000 B4 AlN atoms was thermally equilibrated. Then, a pressure ramp with an average rate ≈7 GPa·ps -1 was used during CMD/NPT simulations for a total duration of ≈0.2 ns. A video of the AlN transformation, produced with the Visual Molecular Dynamics [56] software, can be found in the Supplemental Material.
All 0 K calculations and finite-temperature CMD simulations, described in Sections 2.1 and 2.2, were performed using the LAMMPS package [57] (software updated to version of August 10, 2015). Timesteps of 1 fs (or less) were used for the integration of the equations of motion during CMD and AIMD runs. All CMD data presented in Sections 3 and 4 were obtained using our MEAM Ti 1−x Al x N (0 ≤ x ≤ 1) potential.

Potential Parametrization Results
The MEAM parameters, as determined by the procedure described in Section 2.1, are listed in Table S1 in Supplemental Material and are also available online [58,59]. In addition, we provide scripts that allow a prompt start of CMD simulations for interested readers (see Supplemental Material). The predictions of the potential with respect to the properties of the metallic reference phases (see Section 2.1) are provided in the Supplemental Material file (Section S3, Table S3, Figure S4-S7), where the cohesive energy and the interatomic distances for N 2 and N 3 molecules can also be found (Table S2) (see Supplemental Materials).
The predicted cohesive energies, lattice parameters, and elastic constants for B1-TiN, ε-Ti 2 N, B1-AlN, B3-AlN, and B4-AlN along with their respective experimental (in brackets) and DFT (in parentheses) reference literature values [36, are listed in Table 1. MEAM results for lattice parameters are in excellent agreement with DFT and experimental data. In addition, elastic properties calculated using our MEAM potential are, in general, in very good agreement with the reference values. Furthermore, our potential reproduces the established trend qualitatively with respect to the energetics (i.e., E c ) of AlN polymorph structures; B4-AlN is more stable than B3-AlN, which is, in turn, more stable relative to B1-AlN. From a quantitative point of view, the cohesive energy difference E c,B3-AlN -E c,B4-AlN of 30 meV/atom calculated from MEAM is within the range of DFT values (≈20-40 meV/atom), while the MEAM E c,B1-AlN -E c,B4-AlN value of 68 meV/atom is lower than DFT predictions (≈150-200 meV/atom). Moreover, we find that the MEAM parameters yield an asymmetric B1-Ti 1−x Al x N ∆H mix versus x curve that exhibits a downward concavity and a maximum of ≈230 meV/f.u at x = 0.6. This is within the range of DFT predictions by Alling et al. [3,[90][91][92], Wang et al. [93], and Mayrhofer et al. [94], which found ∆H mix maxima in the interval ≈210-280 meV/f.u. for x ≈ 0. 62 [98,99] are in the range between -3 and 4 eV for N V and -2 and 10 eV for Al V , i.e., in the same order of magnitude with our MEAM values. It is important to note that the formation energy of point defects in semiconductors is largely determined by their respective charge state. Since the MEAM formalism does not account for charge states, a strict quantitative comparison between classical and ab initio data for vacancy formation energies in AlN is not meaningful.
The lattice parameters of B1-Ti 1−x Al x N solid solutions versus x calculated from MEAM at 0 K (see Figure 1) are in excellent agreement with DFT (see figure 2 in Ref. [100] and figure 6 in Ref. [92]) and experimental (see figure 6 in Ref. [92]) results. The elastic constants and Zener's elastic anisotropy factor A = 2·C 44 /(C 11 -C 12 ) obtained from MEAM calculations for B1-Ti 1−x Al x N are plotted as a function of x in Figure 2a. The bulk moduli B and C 11 elastic constantly decrease, while C 44 , C 12 , and the factor A increase monotonically with increasing x. This evolution is in excellent agreement with DFT data [101]. Moreover, our potential yields the same as in DFT calculations [101], stoichiometry (x = 0.28) at which B1-Ti 1−x Al x N solid solutions are elastically isotropic (A = 1). The latter is a strong indication that our MEAM parameters reproduce, in a reliable fashion, interatomic forces and energetics of Ti-Al-N systems. The lattice parameters of B1-Ti1−xAlxN solid solutions versus x calculated from MEAM at 0 K (see Figure 1) are in excellent agreement with DFT (see figure 2 in Ref. [100] and figure 6 in Ref. [92]) and experimental (see figure 6 in Ref. [92]) results. The elastic constants and Zener's elastic anisotropy factor A = 2·C44/(C11-C12) obtained from MEAM calculations for B1-Ti1−xAlxN are plotted as a function of x in Figure 2(a). The bulk moduli B and C11 elastic constantly decrease, while C44, C12, and the factor A increase monotonically with increasing x. This evolution is in excellent agreement with DFT data [101]. Moreover, our potential yields the same as in DFT calculations [101], stoichiometry (x = 0.28) at which B1-Ti1−xAlxN solid solutions are elastically isotropic (A = 1). The latter is a strong indication that our MEAM parameters reproduce, in a reliable fashion, interatomic forces and energetics of Ti-Al-N systems.  This representation is known as Blackman's diagram [102] and entails information on the elastic anisotropy and bonding characteristics of the crystal. The black solid line in Figure 2b represents the condition C12 = C44. The regions above (below) the C12 = C44 line indicate positive (negative) Cauchy's pressures C12 − C44. Based on phenomenological observations, Pettifor suggested that the Cauchy pressure can be used to assess the bonding character in cubic systems: positive values indicate metallic bonds, while negative values indicate more directional and covalent-like bonds [103]. We see that, for B1-Ti1−xAlxN alloys, C12 − C44 become progressively more negative for increasing AlN contents (blue curve in Figure 2b). This is due to the increasing directionality of the bonds caused by the presence of Al, and it is in agreement with DFT results [101].
The melting points Tm of all investigated elementary systems, estimated from CMD simulations, are in reasonable agreement with their reference literature values (see SM). In the case of AlN, we estimated a Tm value of ≈3200 K. To our knowledge, the melting point of AlN has not been experimentally determined, and previous ab initio predictions indicate that AlN melts at approximately 3000 K [60]. For B1-TiN, our MEAM parameters yield a Tm of approximately 6000 K, which is considerably higher than the experimental value (≈3250 K) [104]. It should be emphasized that CMD simulations using NPT sampling in defect-free crystals may lead to overestimations of melting points relative to values predicted using CMD modelling of solid-liquid interfaces at equilibrium [105]. It is reasonable to expect that such overestimations are greater for TiN versus AlN since the absence of polymorphs competing with the B1-TiN at atmospheric pressure reduces the possibility of nucleating heterogeneous sites that can facilitate melting. Figure 2b presents a plot of C 12 /C 11 versus C 44 /C 11 ratios calculated from our MEAM potential. This representation is known as Blackman's diagram [102] and entails information on the elastic anisotropy and bonding characteristics of the crystal. The black solid line in Figure 2b represents the condition C 12 = C 44 . The regions above (below) the C 12 = C 44 line indicate positive (negative) Cauchy's pressures C 12 − C 44 . Based on phenomenological observations, Pettifor suggested that the Cauchy pressure can be used to assess the bonding character in cubic systems: positive values indicate metallic bonds, while negative values indicate more directional and covalent-like bonds [103]. We see that, for B1-Ti 1−x Al x N alloys, C 12 − C 44 become progressively more negative for increasing AlN contents (blue curve in Figure 2b). This is due to the increasing directionality of the bonds caused by the presence of Al, and it is in agreement with DFT results [101].
The melting points T m of all investigated elementary systems, estimated from CMD simulations, are in reasonable agreement with their reference literature values (see Supplemental Material). In the case of AlN, we estimated a T m value of ≈3200 K. To our knowledge, the melting point of AlN has not been experimentally determined, and previous ab initio predictions indicate that AlN melts at approximately 3000 K [60]. For B1-TiN, our MEAM parameters yield a T m of approximately 6000 K, which is considerably higher than the experimental value (≈3250 K) [104]. It should be emphasized that CMD simulations using NPT sampling in defect-free crystals may lead to overestimations of melting points relative to values predicted using CMD modelling of solid-liquid interfaces at equilibrium [105]. It is reasonable to expect that such overestimations are greater for TiN versus AlN since the absence of polymorphs competing with the B1-TiN at atmospheric pressure reduces the possibility of nucleating heterogeneous sites that can facilitate melting.

Potential Validation Results
Classical potentials often lack transferability, i.e., they fail to accurately calculate interatomic forces in chemical environments different than those used in parameter fitting schemes. This means that the ability of the MEAM parameters listed in the supplemental material to reproduce 0 K physical properties of reference phases is a necessary, but not sufficient, condition for describing the dynamic evolution and thermodynamic properties of Ti1−xAlxN (0 ≤ x ≤ 1) alloys. Hence, in this section we present a thorough validation of the proposed set of Ti-Al-N MEAM parameters focusing on finitetemperature properties, which are difficult to implement within the ASA scheme. We study lattice thermal expansion and dynamics, defect migration energetics, and structural stability and transformations in key binary Ti-N and Al-N phases. These properties, in combination with studies of phase energetics in a multitude of Ti1−xAlxN (0 < x < 1) alloy compositions, provide a solid foundation for assessing the robustness and the reliability of the potential beyond the phases and configurations used in the ASA procedure.

Lattice Thermal Expansion
The lattice thermal expansion is an important feature in crystalline solids, as it reflects changes in interatomic forces as a function of temperature. Determination of thermal expansion using CMD simulations is a simple and rapid test for unravelling unphysical lattice instabilities and structural transformations, which may indicate poor potential parameter quality and/or insufficient transferability of the potential formalism.
The temperature-induced variation in B1-Ti1−xAlxN (x = 0, 0.25, 0.5, 0.75, and 1) lattice parameters obtained via CMD simulations is represented in Figure 3. The CMD results (red solid lines in Figure  3a-e) are in good agreement with previous experimental (blue dots in Figure 3a-c) and AIMD simulation data (green solid lines in Figure 3a-e) and black solid line in Figure 3e) [50,71].
Closer inspection of Figure 3a shows that our MEAM parameters yield a room-temperature B1-TiN lattice constant of 4.26 Å, which is slightly larger (<1%) than the corresponding experimental value of 4.24 Å [106]. In addition, static (i.e., 0K) calculations using our MEAM potential show that aB1-TiN(0 K) = 4.252 Å (see Table 1), which is within the range of 0 K DFT calculations (4.188-4.256 Å) based on different electronic exchange/correlation approximations [36]. The relative increase in lattice parameter values calculated for T = 2000 K via CMD simulations is equal to 1.2%, which is comparable to AIMD predictions (2.0%) from Ref. [59] (see solid green line in Figure 3a). The mean linear thermal

Potential Validation Results
Classical potentials often lack transferability, i.e., they fail to accurately calculate interatomic forces in chemical environments different than those used in parameter fitting schemes. This means that the ability of the MEAM parameters listed in the Supplemental Material to reproduce 0 K physical properties of reference phases is a necessary, but not sufficient, condition for describing the dynamic evolution and thermodynamic properties of Ti 1−x Al x N (0 ≤ x ≤ 1) alloys. Hence, in this section we present a thorough validation of the proposed set of Ti-Al-N MEAM parameters focusing on finite-temperature properties, which are difficult to implement within the ASA scheme. We study lattice thermal expansion and dynamics, defect migration energetics, and structural stability and transformations in key binary Ti-N and Al-N phases. These properties, in combination with studies of phase energetics in a multitude of Ti 1−x Al x N (0 < x < 1) alloy compositions, provide a solid foundation for assessing the robustness and the reliability of the potential beyond the phases and configurations used in the ASA procedure.

Lattice Thermal Expansion
The lattice thermal expansion is an important feature in crystalline solids, as it reflects changes in interatomic forces as a function of temperature. Determination of thermal expansion using CMD simulations is a simple and rapid test for unravelling unphysical lattice instabilities and structural transformations, which may indicate poor potential parameter quality and/or insufficient transferability of the potential formalism.
The temperature-induced variation in B1-Ti 1−x Al x N (x = 0, 0.25, 0.5, 0.75, and 1) lattice parameters obtained via CMD simulations is represented in Figure 3. The CMD results (red solid lines in Figure 3a-e) are in good agreement with previous experimental (blue dots in Figure 3a-c) and AIMD simulation data (green solid lines in Figure 3a-e) and black solid line in Figure 3e) [50,71].
Closer inspection of Figure 3a shows that our MEAM parameters yield a room-temperature B1-TiN lattice constant of 4.26 Å, which is slightly larger (<1%) than the corresponding experimental value of 4.24 Å [106]. In addition, static (i.e., 0 K) calculations using our MEAM potential show that a B1-TiN (0 K) = 4.252 Å (see Table 1), which is within the range of 0 K DFT calculations (4.188-4.256 Å) based on different electronic exchange/correlation approximations [36]. The relative increase in lattice parameter values calculated for T = 2000 K via CMD simulations is equal to 1.2%, which is comparable to AIMD predictions (2.0%) from Ref. [59] (see solid green line in Figure 3a). The mean linear thermal expansion coefficient of B1-TiN is thus calculated to monotonically increase from 6.7 × 10 -6 K -1 at 300 K to 8.1 × 10 -6 K -1 at 3000 K. The MEAM estimates are slightly lower than the corresponding experimental values, which lay in the range (7-10) × 10 -6 K -1 [104,[107][108][109]. Even though our potential underestimates the thermal expansion of B1-TiN, the discrepancy between MEAM and experimental a B1-TiN is within ≈1% from room temperature up to 3000 K.
Materials 2018, 11, x FOR PEER REVIEW 9 of 26 expansion coefficient of B1-TiN is thus calculated to monotonically increase from 6.7 × 10 -6 K -1 at 300 K to 8.1 × 10 -6 K -1 at 3000 K. The MEAM estimates are slightly lower than the corresponding experimental values, which lay in the range (7-10) × 10 -6 K -1 [104,[107][108][109]. Even though our potential underestimates the thermal expansion of B1-TiN, the discrepancy between MEAM and experimental aB1-TiN is within ≈1% from room temperature up to 3000 K. Figure 3 also shows that, according to CMD simulations, the variation ∆aB1-TiAlN becomes more pronounced with increasing AlN content x. This is consistent with the results of DFT calculations based on the quasi-harmonic approximation, AIMD simulations [59], and synchrotron X-ray diffraction [109]. For instance, at T = 2000 K, ∆aB1-TiAlN(T) values determined using CMD are found to increase from 1.2 to 3.3% for x increasing from 0 (B1-TiN) to 1 (B1-AlN). For reference, AIMD simulation results [71] showed that ∆aB1-TiAlN(2000 K) increases from 2.0 to 2.5% for the same x range.
Our CMD estimation of the B1-AlN thermal expansion closely matches recent AIMD results (figure 2a in Ref. [50]) obtained for temperatures between 0 and 3000 K (compare black vs. red curves Figure 3e). It should be noted that accurate experimental determination of the B1-AlN thermal expansion is difficult due to the metastable nature of this phase. To our knowledge, the only experimental data available [110], obtained for temperatures between 300 and 450 K, are scattered between 5 and 10 × 10 -6 K -1 . For the same temperature range, MEAM results, Figure 3e, yield a thermal expansion coefficient of ≈13 × 10 -6 K -1 .
The temperature variation of B4-AlN equilibrium volumes between 0 and 3000 K is presented in Figure 4a. There, the B1-AlN data from Figure 3e are plotted again for the sake of clarity. Reproducing the temperature-dependence in equilibrium volumes and lattice parameters of hexagonal wurtzite crystals is particularly challenging due to anisotropic (in-plane vs. out-of-plane) structural properties and interatomic forces [86]. This notwithstanding, our CMD predictions (red solid lines) are in excellent agreement with both experimental (blue solid lines) [86] and AIMD (green solid lines) [50] results. Moreover, both CMD and AIMD curves presented in Figure 4a indicate that B4-AlN has a smaller thermal expansion than B1-AlN, which is consistent with synchrotron X-ray diffraction results from Ref. [109].
The B4-AlN a and c lattice parameters, as well as the c/a ratio are plotted versus T in Figure 4b. The results from our CMD simulations (red solid lines) are within 1% deviation from experimental values (blue solid lines) obtained via x-ray powder diffractometry measurements between 300 and 1400 K [86]. In addition, c/a versus T CMD data are also in agreement with experimental data in Figure 4b.  Figure 3 also shows that, according to CMD simulations, the variation ∆a B1-TiAlN becomes more pronounced with increasing AlN content x. This is consistent with the results of DFT calculations based on the quasi-harmonic approximation, AIMD simulations [59], and synchrotron X-ray diffraction [109]. For instance, at T = 2000 K, ∆a B1-TiAlN (T) values determined using CMD are found to increase from 1.2 to 3.3% for x increasing from 0 (B1-TiN) to 1 (B1-AlN). For reference, AIMD simulation results [71] showed that ∆a B1-TiAlN (2000 K) increases from 2.0 to 2.5% for the same x range.
Our CMD estimation of the B1-AlN thermal expansion closely matches recent AIMD results (figure 2a in Ref. [50]) obtained for temperatures between 0 and 3000 K (compare black vs. red curves Figure 3e). It should be noted that accurate experimental determination of the B1-AlN thermal expansion is difficult due to the metastable nature of this phase. To our knowledge, the only experimental data available [110], obtained for temperatures between 300 and 450 K, are scattered between 5 and 10 × 10 -6 K -1 . For the same temperature range, MEAM results, Figure 3e, yield a thermal expansion coefficient of ≈13 × 10 -6 K -1 .
The temperature variation of B4-AlN equilibrium volumes between 0 and 3000 K is presented in Figure 4a. There, the B1-AlN data from Figure 3e are plotted again for the sake of clarity. Reproducing the temperature-dependence in equilibrium volumes and lattice parameters of hexagonal wurtzite crystals is particularly challenging due to anisotropic (in-plane vs. out-of-plane) structural properties and interatomic forces [86]. This notwithstanding, our CMD predictions (red solid lines) are in excellent agreement with both experimental (blue solid lines) [86] and AIMD (green solid lines) [50] results. Moreover, both CMD and AIMD curves presented in Figure 4a indicate that B4-AlN has a smaller thermal expansion than B1-AlN, which is consistent with synchrotron X-ray diffraction results from Ref. [109].
The B4-AlN a and c lattice parameters, as well as the c/a ratio are plotted versus T in Figure 4b. The results from our CMD simulations (red solid lines) are within 1% deviation from experimental values (blue solid lines) obtained via x-ray powder diffractometry measurements between 300 and 1400 K [86]. In addition, c/a versus T CMD data are also in agreement with experimental data in Figure 4b.  [50] and experiments [86]. The experimental curve (blue) was obtained from the a and c/a ratio variation versus T determined in Ref. [86]. (b) Temperature dependence of a, c, and c/a ratio of B4-AlN lattice parameters as predicted using CMD (red solid line) and measured using X-ray powder diffractometry (blue solid line) in the temperature range 300 to 1400 K from Ref. [86].

Lattice Dynamics
In order to ensure reliability of our MEAM model for large-scale CMD simulations, it is also necessary to test the way in which it reproduces vibrational properties of Ti-Al-N crystals. The study of lattice dynamics using the TDEP enables us to calculate phonon dispersions and phonon densities of states of Ti-Al-N systems directly at finite temperatures. This is useful to, e.g., verify the dynamic stability of crystal structures of interest and assess the effect of vibrational entropy on finitetemperature material properties.
The B1-TiN phonon dispersion curves are presented in Figure 5 for two different temperatures (300 and 1200 K). CMD results are shown in Figure 5a,b, while phonon spectra obtained using AIMD simulations are plotted in Figure 5c,d. For comparison, Figure 5a includes experimental data obtained at room temperature for B1-TiN0.98 compounds. Vibrational frequencies measured via neutron scattering [111] are shown as green circles (transversal modes) and red squares (longitudinal modes). The 300 K CMD dispersion curves (Figure 5a) are in good qualitative agreement with both experimental data ( Figure 5a) and AIMD results of Figure 5c. Moreover, our CMD simulations show that an increase of temperature from 300 (Figure 5a) to 1200 K (Figure 5b) does not significantly affect the vibrational frequencies in B1-TiN, in qualitative agreement with the corresponding AIMD data in Figure 5c,d. On a quantitative level, we observe that CMD predicted optical phonon bandwidths of ≈6 THz (Figure 5a,b), which is larger than the corresponding values from experiments (≈4 THz, Figure 5a) [100] and AIMD (≈2 THz, Figure 5c,d). We also found that CMD simulations cannot not reproduce the phonon softening observed in experimental and AIMD data (Kohn anomalies on acoustic modes at L zone boundaries and on the Γ→X and Γ→K paths), since the screening function employed in the MEAM formalism [21] removes long-range interactions.
Phonon density of states (PDOS) data in Figure 6 demonstrate that our MEAM potential correctly yield dynamically stable (all phonon frequency values are positive, i.e., real) B1-and B4-AlN structures at both room and elevated temperatures (2400 K). The same held for B3-AlN at 300 K (data not shown). Moreover, total PDOS curves and site-projected populations of acoustic and optical states obtained from our simulations are quantitatively close to those obtained via AIMD at both 300 and 2400 K [50]. The dynamic stability evidenced in Figure 6, combined with the accurate reproduction of 0 K cohesive energies in Section 3, shows that our MEAM potential provides a complete description of all AlN polymorphs with a single parameter set.  [50] and experiments [86]. The experimental curve (blue) was obtained from the a and c/a ratio variation versus T determined in Ref. [86]. (b) Temperature dependence of a, c, and c/a ratio of B4-AlN lattice parameters as predicted using CMD (red solid line) and measured using X-ray powder diffractometry (blue solid line) in the temperature range 300 to 1400 K from Ref. [86].

Lattice Dynamics
In order to ensure reliability of our MEAM model for large-scale CMD simulations, it is also necessary to test the way in which it reproduces vibrational properties of Ti-Al-N crystals. The study of lattice dynamics using the TDEP enables us to calculate phonon dispersions and phonon densities of states of Ti-Al-N systems directly at finite temperatures. This is useful to, e.g., verify the dynamic stability of crystal structures of interest and assess the effect of vibrational entropy on finite-temperature material properties.
The B1-TiN phonon dispersion curves are presented in Figure 5 for two different temperatures (300 and 1200 K). CMD results are shown in Figure 5a,b, while phonon spectra obtained using AIMD simulations are plotted in Figure 5c,d. For comparison, Figure 5a includes experimental data obtained at room temperature for B1-TiN 0.98 compounds. Vibrational frequencies measured via neutron scattering [111] are shown as green circles (transversal modes) and red squares (longitudinal modes). The 300 K CMD dispersion curves (Figure 5a) are in good qualitative agreement with both experimental data ( Figure 5a) and AIMD results of Figure 5c. Moreover, our CMD simulations show that an increase of temperature from 300 (Figure 5a) to 1200 K (Figure 5b) does not significantly affect the vibrational frequencies in B1-TiN, in qualitative agreement with the corresponding AIMD data in Figure 5c,d. On a quantitative level, we observe that CMD predicted optical phonon bandwidths of ≈6 THz (Figure 5a,b), which is larger than the corresponding values from experiments (≈4 THz, Figure 5a) [100] and AIMD (≈2 THz, Figure 5c,d). We also found that CMD simulations cannot not reproduce the phonon softening observed in experimental and AIMD data (Kohn anomalies on acoustic modes at L zone boundaries and on the Γ→X and Γ→K paths), since the screening function employed in the MEAM formalism [21] removes long-range interactions.
Phonon density of states (PDOS) data in Figure 6 demonstrate that our MEAM potential correctly yield dynamically stable (all phonon frequency values are positive, i.e., real) B1-and B4-AlN structures at both room and elevated temperatures (2400 K). The same held for B3-AlN at 300 K (data not shown). Moreover, total PDOS curves and site-projected populations of acoustic and optical states obtained from our simulations are quantitatively close to those obtained via AIMD at both 300 and 2400 K [50]. The dynamic stability evidenced in Figure 6, combined with the accurate reproduction of 0 K cohesive energies in Section 3, shows that our MEAM potential provides a complete description of all AlN polymorphs with a single parameter set.  CMD acoustic phonon frequencies calculated for B1-and B4-AlN at 300 and 2400 K are, overall, in good agreement with AIMD results (compare Figure 6 with figure 4 in Ref. [50]). It is important to note that MEAM cannot reproduce the splitting of degenerate longitudinal and transverse optical phonons [112], which is inherent to polarizable materials such as AlN (see, for example, optical frequencies near the Γ point in Figure 6 and those obtained ab initio in Ref. [50]). This is due to the fact that the MEAM formalism does not explicitly include Coulomb interactions.

Point-Defect Migration Energies
The set of MEAM parameters selected via the ASA procedure was also tested with respect to the energetics of point-defect migration in Ti-N and Al-N systems. This is motivated by the fact that a number of structural transformations in solid solutions, including spinodal decomposition of B1-Ti1−xAlxN, is believed to be kinetically controlled by diffusion of lattice vacancies [6,113].   CMD acoustic phonon frequencies calculated for B1-and B4-AlN at 300 and 2400 K are, overall, in good agreement with AIMD results (compare Figure 6 with figure 4 in Ref. [50]). It is important to note that MEAM cannot reproduce the splitting of degenerate longitudinal and transverse optical phonons [112], which is inherent to polarizable materials such as AlN (see, for example, optical frequencies near the Γ point in Figure 6 and those obtained ab initio in Ref. [50]). This is due to the fact that the MEAM formalism does not explicitly include Coulomb interactions.

Point-Defect Migration Energies
The set of MEAM parameters selected via the ASA procedure was also tested with respect to the energetics of point-defect migration in Ti-N and Al-N systems. This is motivated by the fact that a number of structural transformations in solid solutions, including spinodal decomposition of B1-Ti1−xAlxN, is believed to be kinetically controlled by diffusion of lattice vacancies [6,113]. CMD acoustic phonon frequencies calculated for B1-and B4-AlN at 300 and 2400 K are, overall, in good agreement with AIMD results (compare Figure 6 with figure 4 in Ref. [50]). It is important to note that MEAM cannot reproduce the splitting of degenerate longitudinal and transverse optical phonons [112], which is inherent to polarizable materials such as AlN (see, for example, optical frequencies near the Γ point in Figure 6 and those obtained ab initio in Ref. [50]). This is due to the fact that the MEAM formalism does not explicitly include Coulomb interactions.

Point-Defect Migration Energies
The set of MEAM parameters selected via the ASA procedure was also tested with respect to the energetics of point-defect migration in Ti-N and Al-N systems. This is motivated by the fact that a number of structural transformations in solid solutions, including spinodal decomposition of B1-Ti 1−x Al x N, is believed to be kinetically controlled by diffusion of lattice vacancies [6,113]. Figure 7 presents the results of the 0 K MEAM-NEB calculations for metal (Ti V and Al V ) and N (N V ) vacancy migration in B1-TiN (Figure 7a), B1-AlN (Figure 7b), and B4-AlN (Figure 7c,d). We found that the migration energies for Ti V and N V in B1-TiN along <110> directions are 4.13 and 4.24 eV, respectively, with saddle-point transition states located halfway between initial and final atomic positions (Figure 7a). The MEAM predictions are in good agreement with DFT-NEB results, which yield energy barriers of ≈3.8 eV for N V [114] and 4.26 eV for Ti V [97]. Experimental values for vacancy migration energies in B1-TiN are scarce in the literature. Kodambaka et al. [115] determined a global (i.e., both Ti V and N V ) value of 4.5 ± 0.2 eV for vacancy diffusion energies in bulk B1-TiN. Hultman et al. [116] estimated activation energies for metal interdiffusion at TiN/ZrN superlattice interfaces that range between 2.6 and 4.5 eV, while Wood and Paasche [117] and Anglezio-Abautret et al. [118] reported that activation barriers for N diffusion in B1-TiN lie in the range 1.8 to 5.5 eV.  (Figure 7b), and B4-AlN (Figure 7c,d). We found that the migration energies for TiV and NV in B1-TiN along <110> directions are 4.13 and 4.24 eV, respectively, with saddle-point transition states located halfway between initial and final atomic positions (Figure 7a). The MEAM predictions are in good agreement with DFT-NEB results, which yield energy barriers of ≈3.8 eV for NV [114] and 4.26 eV for TiV [97]. Experimental values for vacancy migration energies in B1-TiN are scarce in the literature. Kodambaka et al. [115] determined a global (i.e., both TiV and NV) value of 4.5 ± 0.2 eV for vacancy diffusion energies in bulk B1-TiN. Hultman et al. [116] estimated activation energies for metal interdiffusion at TiN/ZrN superlattice interfaces that range between 2.6 and 4.5 eV, while Wood and Paasche [117] and Anglezio-Abautret et al. [118] reported that activation barriers for N diffusion in B1-TiN lie in the range 1.8 to 5.5 eV. Besides TiV and NV, we also studied formation and diffusion of Al interstitials (AlI) in B1-TiN. Using the energy of an Al atom in fcc-Al as a chemical potential, we calculated that the energy required for forming an AlI at tetrahedral sites is 4.92 eV, which is in reasonable agreement with the DFT value of 3.81 by Mei et al. [96]. Then, the corresponding AlI migration energy across tetrahedral B1-TiN sites was found to be equal to 2.42 eV, which matches perfectly the DFT value from Ref. [96].
Our calculations for NV diffusion within and across the B4-AlN (0001) plane yield activation energies of 1.97 and 2.19 eV, respectively (see Figure 7c). These values are within the experimental uncertainty range of O and N interdiffusion activation energies at Al2O3/AlN interfaces (2.49 ± 0.42 eV) [119]. For AlV, we found that migration across the (0001) B4-AlN lattice planes requires an activation energy of 2.29 eV, which is lower than the value of 2.72 eV for in-plane diffusion ( Figure  7d). Moreover, our potential predicts that AlV transport in B1-AlN occurs with an activation energy of 2.47 eV, which is similar to that required in the B4-AlN polymorph. Concurrently, NV migration in B1-AlN requires significantly higher activation energy (4.00 eV) than that in the B4 structure ( Figure  7b). We also note that no experimental and/or theoretical data on diffusion of point defects in B1-and B4-AlN are available in the literature to compare with our MEAM results. Besides Ti V and N V , we also studied formation and diffusion of Al interstitials (Al I ) in B1-TiN. Using the energy of an Al atom in fcc-Al as a chemical potential, we calculated that the energy required for forming an Al I at tetrahedral sites is 4.92 eV, which is in reasonable agreement with the DFT value of 3.81 by Mei et al. [96]. Then, the corresponding Al I migration energy across tetrahedral B1-TiN sites was found to be equal to 2.42 eV, which matches perfectly the DFT value from Ref. [96].
Our calculations for N V diffusion within and across the B4-AlN (0001) plane yield activation energies of 1.97 and 2.19 eV, respectively (see Figure 7c). These values are within the experimental uncertainty range of O and N interdiffusion activation energies at Al 2 O 3 /AlN interfaces (2.49 ± 0.42 eV) [119]. For Al V , we found that migration across the (0001) B4-AlN lattice planes requires an activation energy of 2.29 eV, which is lower than the value of 2.72 eV for in-plane diffusion (Figure 7d). Moreover, our potential predicts that Al V transport in B1-AlN occurs with an activation energy of 2.47 eV, which is similar to that required in the B4-AlN polymorph. Concurrently, N V migration in B1-AlN requires significantly higher activation energy (4.00 eV) than that in the B4 structure (Figure 7b). We also note that no experimental and/or theoretical data on diffusion of point defects in B1-and B4-AlN are available in the literature to compare with our MEAM results.
Ab initio data for point-defect diffusion in Ti-Al-N solid solutions are not available in the literature. Nevertheless, experimental studies of structural evolution in annealed B1-Ti 1−x Al x N samples by Mayrhofer et al. [6] and Norrby et al. [120] attributed activation energies of ≈3.3 eV (in Ti 0.36 Al 0.64 N) and ≈3.6 eV (in Ti 0.55 Al 0.45 N) to spinodal decomposition and B1-to-B4 transformations within AlN-rich domains, which in turn, is primarily attributed to diffusion of metal and nitrogen vacancies. The experimental estimates of atomic migration energies during spinodal decomposition in Ti-Al-N (3.3-3.6 eV) are consistent with the range of values that we obtain for cation and anion diffusion in binary Ti-N and Al-N.
Although determination of point-defect formation and migration energies in Ti-Al-N solid solutions is strongly dependent on the local chemical environment [121] (and therefore lies outside the scope of this work), we used CMD simulations to ensure stability of B1-Ti 1−x Al x N alloys containing defects, i.e., monovacancies, divacancies, interstitials, and interstitial pairs at different temperatures. We found that point defect migration and point-defect/point-defect interactions do not cause unphysical structural transformations within the alloys during the investigated time scales, which are of the order of one nanosecond. This, together with the results presented above, lends confidence that our model potential is suitable to investigate phase transformation phenomena in Ti-Al-N solid solutions.

Equilibrium Volumes and Elastic Properties of B1-TiN y , B2-TiN, and bct-Ti 2 N
The proposed set of MEAM parameters (see Supplemental Material) was carefully fitted to the properties of stoichiometric B1-TiN and ε-Ti 2 N. To verify transferability of our model to a variety of Ti-N lattice configurations and bonding geometries, which may be encountered during simulations of dynamic processes in Ti-Al-N alloys, we present here the elastic and structural properties calculated for nitrogen-deficient B1-TiN y (0.69 ≤ y < 1) as well as high-pressure B2-TiN and bct-Ti 2 N phases.
Ti-N compounds maintain the cubic B1 structure over a wide compositional range [122]. In B1-TiN y , understoichiometry (y < 1) is primarily accommodated by anion vacancies [123]. In addition, control of the N content during synthesis can be used to tune the TiN y optical [124], electrical [106], and mechanical [84,125] properties. In Figure 8, we plot the lattice parameter a TiNy (Figure 8a), Young's modulus E (Figure 8b), and elastic constants C 11 and C 44 (Figure 8c,d, respectively) as predicted using 0 K calculations with our MEAM potential along with experimental and DFT data for comparison. We observe (Figure 8a) that a TiNy increase monotonically with increasing the N content y, in excellent agreement (maximum discrepancy ≈0.5%) with experimental measurements [84]. Remarkably, our MEAM potential reproduces the experimental a TiNy versus y trend better than DFT [126]. Other results (not included in Figure 8a) show that, depending on the choice of the electronic exchange-correlation approximation, DFT overestimates or underestimates the lattice parameter of stoichiometric TiN by up to 1.3% [36]. DFT calculations [126,127] and experiments (acoustic wave velocities and nanoindentation tests [84,85,128]) have demonstrated that E, C 11 , and C 44 , in B1-TiN y increase monotonically as function of y, as shown in Figure 8b

Phase Stability and Transitions
The cohesive energies of B1-TiN and ε-Ti2N phases, which are the ground-state configurations for the TiN and Ti2N systems, were reproduced via the ASA parametrization (see Section 3 and Table  1). Post-parametrization calculations were carried out to verify that our model predicted correct trends in the energetics for high-pressure metastable B2-TiN and bct-Ti2N polymorph structures. The MEAM potential predicts that the B1-TiN structure was ≈1.6 eV/atom more stable than the B2-TiN phase, in fair agreement with present DFT calculations which yield an energy difference of ≈0.9 eV/atom between the two phases. Moreover, we find that ε-Ti2N is 98 meV/atom more stable than the bct-Ti2N polymorph, which is qualitatively consistent with DFT results that have shown an energy difference of 16 meV/atom in favor of the ε-Ti2N structure [129].

Al-N.
At atmospheric pressure and for T = 300 K, B4 is the thermodynamically stable AlN structure. However, subject to compression, polycrystalline B4-AlN samples transform into the metastable B1-AlN polymorph [67,74,130]. Experiments indicate that the B4-to-B1 transition pressure decreases from ≈14 GPa (at 300 K) to ≈11-12 GPa at temperatures between 1000 and 2000 K [67,74,130]. For comparison, DFT-based studies report a transition pressure that decreases monotonically with increasing T from ≈13-17 GPa at 0 K to ≈6 GPa at 3000 K [60,77,131]. Other first-principles investigations based on the quasi-harmonic approximation [132] demonstrate that the B4-to B1-AlN transition pressures predicted via DFT strongly depend on the approximation used for the electronic exchange-correlation energy: the B4/B1 phase boundary calculated ab initio with two different exchange-correlation functionals decreases, for temperatures increasing from 0 to ≈2200 K, from 7 to 5 GPa and from 12 to 10 GPa (figure 2 in Ref. [132]). Moreover, the AlN phase diagram obtained via AIMD simulations indicates that the transition pressure changes from 13 GPa (at 0 K) to 0 GPa (at 3200 K) (see figure 1 in Ref. [50]). The free energies of B1-and B4-AlN calculated as a function of temperature and pressure via CMD (see Section 2.2 for details) show that, at 300 K, the B1-AlN polymorph structure becomes the most thermodynamically stable at a pressure of ≈5 GPa, which is quite below the experimentally-determined value (≈14 GPa). CMD results yield a B4/B1 phase

Ti-N
The cohesive energies of B1-TiN and ε-Ti 2 N phases, which are the ground-state configurations for the TiN and Ti 2 N systems, were reproduced via the ASA parametrization (see Section 3 and Table 1). Post-parametrization calculations were carried out to verify that our model predicted correct trends in the energetics for high-pressure metastable B2-TiN and bct-Ti 2 N polymorph structures. The MEAM potential predicts that the B1-TiN structure was ≈1.6 eV/atom more stable than the B2-TiN phase, in fair agreement with present DFT calculations which yield an energy difference of ≈0.9 eV/atom between the two phases. Moreover, we find that ε-Ti 2 N is 98 meV/atom more stable than the bct-Ti 2 N polymorph, which is qualitatively consistent with DFT results that have shown an energy difference of 16 meV/atom in favor of the ε-Ti 2 N structure [129].

Al-N
At atmospheric pressure and for T = 300 K, B4 is the thermodynamically stable AlN structure. However, subject to compression, polycrystalline B4-AlN samples transform into the metastable B1-AlN polymorph [67,74,130]. Experiments indicate that the B4-to-B1 transition pressure decreases from ≈14 GPa (at 300 K) to ≈11-12 GPa at temperatures between 1000 and 2000 K [67,74,130]. For comparison, DFT-based studies report a transition pressure that decreases monotonically with increasing T from ≈13-17 GPa at 0 K to ≈6 GPa at 3000 K [60,77,131]. Other first-principles investigations based on the quasi-harmonic approximation [132] demonstrate that the B4-to B1-AlN transition pressures predicted via DFT strongly depend on the approximation used for the electronic exchange-correlation energy: the B4/B1 phase boundary calculated ab initio with two different exchange-correlation functionals decreases, for temperatures increasing from 0 to ≈2200 K, from 7 to 5 GPa and from 12 to 10 GPa (figure 2 in Ref. [132]). Moreover, the AlN phase diagram obtained via AIMD simulations indicates that the transition pressure changes from 13 GPa (at 0 K) to 0 GPa (at 3200 K) (see figure 1 in Ref. [50]). The free energies of B1-and B4-AlN calculated as a function of temperature and pressure via CMD (see Section 2.2 for details) show that, at 300 K, the B1-AlN polymorph structure becomes the most thermodynamically stable at a pressure of ≈5 GPa, which is quite below the experimentally-determined value (≈14 GPa). CMD results yield a B4/B1 phase boundary that remains near ≈4 GPa up to the AlN melting point estimated by MEAM (T m ≈ 3200 K). This is consistent with the result that MEAM underestimates the energy difference E c,B4-AlN -E c,B1-AlN (see Table 1).
Despite the fact that the AlN phase diagram based upon CMD free-energy calculations is only in reasonable agreement with experimental and ab initio findings, our model reproduces the correct dynamics of the B4-to-B1 AlN transition. This is seen in Figure 9, which shows the evolution of stress in a B4-AlN crystal along the three orthogonal directions upon applying external pressure that increases at an average rate of ≈7 GPa·ps -1 , at 300 K (see more details in Section 2.2). After a monotonic increase up to a value of 110 ± 10 GPa, the stress accumulated in B4-AlN is partially relieved. As we discuss in detail below, this stress relief is due to the formation of small B1-AlN grains.
The pressure value at which the B4-to B1-AlN phase transition was initiated in our CMD simulations is in excellent agreement with the results of AIMD modeling (100-120 GPa, see figure 1 in Ref. [133]). The B4-to-B1 AlN transition pressure values obtained during present CMD and previous AIMD [133] simulations at 300 K are considerably above the interval of pressures (5-14 GPa) predicted in our AlN phase diagram, as well as in AIMD and experimental AlN phase diagrams [50,74]. This is an expected discrepancy; while CMD simulations are carried out for defect-free crystals during extremely short times (≈10 -10 s), compression experiments are conducted at close-to-equilibrium conditions (≥10 2 s) and on polycrystalline B4-AlN samples. Hence, pressures much larger than those predicted by thermodynamics (phase diagrams) are necessary to quickly activate the B4-to-B1 AlN transition during the timescales accessible via molecular dynamics simulations. The fact that our CMD results (Figure 9) match AIMD predictions [133] suggests that the kinetic free-energy barrier that separates the B4-and B1-AlN phases is accurately reproduced by MEAM. . This is consistent with the result that MEAM underestimates the energy difference Ec,B4-AlN -Ec,B1-AlN (see Table 1). Despite the fact that the AlN phase diagram based upon CMD free-energy calculations is only in reasonable agreement with experimental and ab initio findings, our model reproduces the correct dynamics of the B4-to-B1 AlN transition. This is seen in Figure 9, which shows the evolution of stress in a B4-AlN crystal along the three orthogonal directions upon applying external pressure that increases at an average rate of ≈7 GPa·ps -1 , at 300 K (see more details in Section 2.2). After a monotonic increase up to a value of 110 ± 10 GPa, the stress accumulated in B4-AlN is partially relieved. As we discuss in detail below, this stress relief is due to the formation of small B1-AlN grains.
The pressure value at which the B4-to B1-AlN phase transition was initiated in our CMD simulations is in excellent agreement with the results of AIMD modeling (100-120 GPa, see figure 1 in Ref. [133]). The B4-to-B1 AlN transition pressure values obtained during present CMD and previous AIMD [133] simulations at 300 K are considerably above the interval of pressures (5)(6)(7)(8)(9)(10)(11)(12)(13)(14) predicted in our AlN phase diagram, as well as in AIMD and experimental AlN phase diagrams [50,74]. This is an expected discrepancy; while CMD simulations are carried out for defect-free crystals during extremely short times (≈10 -10 s), compression experiments are conducted at close-toequilibrium conditions (≥10 2 s) and on polycrystalline B4-AlN samples. Hence, pressures much larger than those predicted by thermodynamics (phase diagrams) are necessary to quickly activate the B4to-B1 AlN transition during the timescales accessible via molecular dynamics simulations. The fact that our CMD results (Figure 9) match AIMD predictions [133] suggests that the kinetic free-energy barrier that separates the B4-and B1-AlN phases is accurately reproduced by MEAM.  Figure 10 presents the detailed structural evolution observed in AlN under compression during CMD simulations (full movie available in the SM, Video S1). The defect-free B4-AlN structure is retained up to a pressure of ≈100 GPa (simulation time of 16 ps). The occurrence of local lattice slip at a pressure of ≈110 GPa (shown in the magnified B4-AlN regions at 16.6 ps in Figure 10) leads to the formation of small B1-AlN grains, as seen in the insets of the snapshot for 16.7 ps. This causes a pressure drop from ≈110 to ≈70 GPa. From that point, B1-AlN grains continue growing at the expense of the B4-AlN crystal. At approximately 20 ps (≈85 GPa), a single phase polycrystalline B1-AlN structure is formed. Thereupon, smaller B1-AlN grains coalesce with the larger ones, see e.g, the areas enclosed in the white brackets in the simulation snapshots corresponding to times between 24 and 36 ps (≈140-330 GPa). At a simulation time of ≈37 ps (≈350 GPa), the original defect-free B4-AlN has transformed into defective B1-AlN crystal. The B1-AlN structure is retained up to pressures of ≈1500 GPa (not shown in Figure 10), which is the largest pressure used in our simulations. This result is The abrupt drop in stress at ≈17 s (≈110 GPa) corresponding to the formation of small B1-AlN grains, which is the first step for the B4-to B1-AlN phase transformation. Figure 10 presents the detailed structural evolution observed in AlN under compression during CMD simulations (full movie available in the Supplementray Materials, Video S1). The defect-free B4-AlN structure is retained up to a pressure of ≈100 GPa (simulation time of 16 ps). The occurrence of local lattice slip at a pressure of ≈110 GPa (shown in the magnified B4-AlN regions at 16.6 ps in Figure 10) leads to the formation of small B1-AlN grains, as seen in the insets of the snapshot for 16.7 ps. This causes a pressure drop from ≈110 to ≈70 GPa. From that point, B1-AlN grains continue growing at the expense of the B4-AlN crystal. At approximately 20 ps (≈85 GPa), a single phase polycrystalline B1-AlN structure is formed. Thereupon, smaller B1-AlN grains coalesce with the larger ones, see e.g., the areas enclosed in the white brackets in the simulation snapshots corresponding to times between 24 and 36 ps (≈140-330 GPa). At a simulation time of ≈37 ps (≈350 GPa), the original defect-free B4-AlN has transformed into defective B1-AlN crystal. The B1-AlN structure is retained up to pressures of ≈1500 GPa (not shown in Figure 10), which is the largest pressure used in our simulations. This result is consistent with X-ray diffraction observations, which excluded the occurrence (at equilibrium conditions) of a second solid-to-solid AlN phase transition at pressures up to 132 GPa [79].

Ti-Al-N
B1-Ti 1−x Al x N solid solutions with high (x > 0.5) AlN contents are known to decompose via the spinodal route, at temperatures near 1200 K [7,8], leading to formation of coherent Ti-and Al-rich B1-Ti 1−x Al x N domains. This hinders dislocation motion across strained domain interfaces and results in the age-hardening effect, which is of extreme importance for metal cutting at elevated temperatures [134,135]. A further increase in temperature causes B1-to-B4 structural transformation within Al-rich domains, which is detrimental for hardness and coating performance [7,8,134,135]. To date, the atomistic mechanisms that control spinodal decomposition and subsequent formation of B4-AlN in Ti-Al-N coatings are unclear. These processes can be elucidated by means of large-scale CMD simulations using our MEAM potential. The suitability of the potential for such simulations is explored in this section by discussing the B1-Ti 1−x Al x N phase diagram predicted using MEAM, in comparison to first-principle and experimental results.
The B1-Ti 1−x Al x N binodal and spinodal curves calculated on the basis of CMD free energies of mixing are presented in Figure 11. It is important to note that our phase diagram ends at the AlN melting point (T m ≈ 3200 K); the temperature beyond which the free-energy of pure B1-AlN (necessary to compute the alloy free energy of mixing) is not defined. This, however, does not imply that B1-AlN-rich regions may not form within B1-Ti 1−x Al x N solid solutions during CMD simulations at T > 3200 K. Our estimated binodal and spinodal regions, Figure 11, are in very good agreement with recent thermodynamic assessments (see figure 2 in Ref. [136]), as well as ab initio calculations accounting for vibrational effects on alloy free energies (see figure 3 in Ref. [93]). Our CMD results also match reasonably well with the Ti-Al-N phase diagram determined via AIMD simulations (see figure 2 in Ref. [137]).

Ti-Al-N.
B1-Ti1−xAlxN solid solutions with high (x > 0.5) AlN contents are known to decompose via the spinodal route, at temperatures near 1200 K [7,8], leading to formation of coherent Ti-and Al-rich B1-Ti1−xAlxN domains. This hinders dislocation motion across strained domain interfaces and results in the age-hardening effect, which is of extreme importance for metal cutting at elevated temperatures [134,135]. A further increase in temperature causes B1-to-B4 structural transformation within Al-rich domains, which is detrimental for hardness and coating performance [7,8,134,135]. To date, the atomistic mechanisms that control spinodal decomposition and subsequent formation of B4-AlN in Ti-Al-N coatings are unclear. These processes can be elucidated by means of large-scale CMD simulations using our MEAM potential. The suitability of the potential for such simulations is explored in this section by discussing the B1-Ti1−xAlxN phase diagram predicted using MEAM, in comparison to first-principle and experimental results.
The B1-Ti1−xAlxN binodal and spinodal curves calculated on the basis of CMD free energies of mixing are presented in Figure 11. It is important to note that our phase diagram ends at the AlN melting point (Tm ≈ 3200 K); the temperature beyond which the free-energy of pure B1-AlN (necessary to compute the alloy free energy of mixing) is not defined. This, however, does not imply that B1-AlN-rich regions may not form within B1-Ti1−xAlxN solid solutions during CMD simulations at T > 3200 K. Our estimated binodal and spinodal regions, Figure 11, are in very good agreement with recent thermodynamic assessments (see figure 2 in Ref. [136]), as well as ab initio calculations accounting for vibrational effects on alloy free energies (see figure 3 in Ref. [93]). Our CMD results also match reasonably well with the Ti-Al-N phase diagram determined via AIMD simulations (see figure 2 in Ref. [137]). Figure 11. Spinodal and binodal curves for B1-Ti1−xAlxN solid solutions calculated using CMD free energies of mixing with respect to B1-TiN and B1-AlN. The inset compares our results to experimental data (black and red stars mark B1-Ti1−xAlxN and spinodally-decomposed solid solutions, respectively) collected in Ref. [138].
Previous DFT calculations, in which temperature effects on alloy mixing free energies were neglected or added as mean-field theory corrections to mixing enthalpies determined at 0 K [3,91,92], Figure 11. Spinodal and binodal curves for B1-Ti 1−x Al x N solid solutions calculated using CMD free energies of mixing with respect to B1-TiN and B1-AlN. The inset compares our results to experimental data (black and red stars mark B1-Ti 1−x Al x N and spinodally-decomposed solid solutions, respectively) collected in Ref. [138].
Previous DFT calculations, in which temperature effects on alloy mixing free energies were neglected or added as mean-field theory corrections to mixing enthalpies determined at 0 K [3,91,92], estimated a consolute temperature T c in the range ≈7500-9000 K. Subsequent first-principles investigations by Wang et al. [93] showed that accounting for vibrational entropies S vib in Ti-Al-N free energies of mixing via the Debye-Grüneisen approximation significantly lowered the calculated T c (≈3800 K). More recently, AIMD results by Shulumba et al. demonstrated that the explicit inclusion of anharmonic effects in S vib further reduces the predicted T c to 2900 K [137]. Consistent with the observations of Refs. [93,137], our CMD-based evaluations also show that the implicit inclusion of S vib yields narrower miscibility gaps at temperatures above ≈1500 K, as compared to those close to room temperature.
At elevated temperatures, the CMD phase diagram presented in Figure 11 exhibits discrepancies with previous ab initio [93,137] and Calculation of Phase Diagrams (CALPHAD) [136] results. For example, at a temperature of 2500 K, the spinodal region predicted by the MEAM potential spans in the AlN content x range 0.35 to 1, while Refs. [93,136,137], indicate spinodal regions for 0.4 < x < 0.9, 0.4 < x < 0.9, and 0.6 < x < 0.9, respectively. Nonetheless, consistent with several first-principles results [3,[91][92][93]136,137], our Ti-Al-N phase diagram is skewed toward the AlN end, indicating that TiN-rich alloys are more stable than their AlN-rich counterparts. In stark contrast with the predictions of our model and Refs. [3,[91][92][93]136,137], another recent thermodynamic evaluation [138] reports a B1-Ti 1-x Al x N phase diagram which is skewed toward the TiN end and with the lowest T c (≈2100-2200 K) calculated so far (see figure 1 in Ref. [138]).
In their paper [138], Zhou et al. provide a detailed overview of experimental results collected during the past two decades for annealed B1-Ti 1-x Al x N solid solutions. These are included for comparison with our results in the inset of Figure 11, where black stars denote stable or metastable cubic alloy compositions, and red stars mark spinodally-decomposed Ti-Al-N samples. A first apparent inconsistency between theoretical and experimental assessments of miscibility gaps can be seen for T ≤ 1000 K; the experimental observations (black stars) indicate stability (or metastability) of B1-Ti 1−x Al x N alloys for AlN concentrations well within the spinodal compositional range (0.25 < x < 1) predicted by the CMD phase diagram. Refractory metal-nitride compounds exhibit relatively low atomic mobilities [139], which may shift the onset of structural transformation to temperatures larger than those predicted in theoretical models. In order to explore possible effects of low point-defect mobilities on the outcome of annealing experiments conducted on Ti-Al-N, we consider the case of the parent compound TiN, for which nitrogen diffusivities D(T) are available in the literature [36] (in TiN, Ti migration is even slower than N migration [97]). We consider the model system B1-TiN 1-y with nitrogen vacancy concentrations at the dilute limit (y ≈ 10 -5 ) and diffusivities D(1100 K) ≈ 10 -20±2 cm 2 ·s -1 and D(800 K) ≈ 10 -33±2 cm 2 ·s -1 , as shown in figure 11 of Ref. [36]. Annealing B1-TiN 1−y at a temperature T = 1100 K during one day (t ≈ 10 5 s) would induce nitrogen transport over distances d ≈ [6·t·D(1000 K)] 1 2 of a few nm. A modestly lower annealing temperature (800 K) would require considerably longer times (years) to observe diffusion over comparable length scales. Given the similarity in chemical bonds and crystal structures, it is reasonable to assume that nitrogen and metal atoms in Ti-Al-N migrate at rates comparable to those discussed above for TiN, thus providing a possible explanation for the retained stability of B1-Ti 1−x Al x N solid solutions at T ≤ 1000 K (inset in Figure 11).
The experimental results are in consistent with theory for temperatures between 1000 K and ≈1500-1800 K: (i) no decomposition was observed for x < ≈0.2 (black stars); (ii) spinodal regions (red stars) are within the range ≈0.25 < x < ≈0.8. Nonetheless, for T > 1700 K, the experimental observations [138] indicate that Ti-Al-N solid solutions with AlN contents in the range ≈0.25 < x < ≈0.6 maintain the cubic structure. At such high temperatures, the scenario of spinodal decomposition being kinetically limited is not valid. This may indicate that the miscibility gap closes at temperatures as low as ≈1700 K, which is well below the consolute temperature predicted in all theoretical TiAlN phase diagrams (including the one in the present study) reported in the literature.

Summary and Outlook
In this work, we presented a MEAM interatomic potential for the Ti 1−x Al x N (0 ≤ x ≤ 1) alloy system. This was achieved by developing an adaptive simulated annealing (ASA) algorithm for scanning the MEAM multi-dimensional parameter space and identifying a set of parameters that optimize the potential's predictions relative to 0 K experimentally-and ab-initio-determined equilibrium volumes, cohesive energies, elastic constants, defect formation energies, and mixing enthalpies of elemental, binary, and ternary phases and configurations within the Ti-Al-N system. We then validated the transferability of the parameters selected via the ASA procedure by establishing the ability of the potential to reproduce known finite-temperature kinetic and thermodynamic properties, not explicitly included in the ASA optimization, for B1-Ti 1−x Al x N (0 ≤ x ≤ 1), B4-AlN, and various Ti-N phases. We found that, overall, the potential reproduces the following well: (i) temperature-dependence of equilibrium volumes in B1-Ti 1−x Al x N and B4-AlN; (ii) phonon dispersion curves of B1-TiN, B1-AlN, and B4-AlN; (iii) point-defect migration energies in B1-TiN, B1-AlN, and B4-AlN; (iv) lattice and elastic constants of sub-stoichiometric B1-TiN y (0.7 ≤ y ≤ 1), B2-TiN, and bct-Ti 2 N; (v) free energies of B1-AlN and B4-AlN, and dynamics of B4-to B1-AlN pressure-induced phase transformation; and (vi) B1-Ti 1−x Al x N alloy mixing free energies and binodal/spinodal phase boundaries with respect to B1-TiN and B1-AlN.
The outcome of the present study is an interatomic potential that can enable large-scale classical molecular dynamics simulations in one of the most important engineering materials systems, i.e., Ti-Al-N. These simulations may include thermal stability of B1-Ti 1−x Al x N solid solutions, multilayers, and self-organized nanostructures and can provide critical insights onto the atomistic pathways that drive segregation via the spinodal route and phase transformation via nucleation and growth. Our MEAM potential can also be used to simulate elastic and plastic responses of B1-Ti 1−x Al x N alloys, which together with thermal stability simulations can guide the design of protective coatings with superior mechanical and oxidation performance. Other areas where our MEAM set of parameters can be used include diffusion of Al in B1-Ti 1−x Al x N used as barrier layer and contact material in microelectronic devices [140][141][142]. Moreover, the present set of MEAM parameters can be used as starting point for developing potentials that can simulate vapor-based growth of Ti-Al-N films. Beyond Ti-Al-N, the methodologies presented in this work can serve as the foundation for developing potentials for quaternary systems, e.g., Ti-Al-Si-N, Ti-Al-Ta-N, and Ti-Al-Nb-N, where complex metal-sublattice configurations are used to enhance surface wear and oxidation resistance, and inhibit wurtzite-phase formation subsequent to spinodal decomposition, which allows extending the age hardening effect to temperatures above 1200 K [143][144][145][146][147][148][149][150].

Conflicts of Interest:
The authors declare no conflict of interest.