Using Machine-Learned Force Fields for Describing Heat-Transport Related Quantities in AlGaN and Derived Materials

1. Introduction

The rapid development of microelectronic devices, coupled with the ongoing demand for higher integration density and smaller feature sizes, presents significant design challenges. Effective management of waste heat has become increasingly critical, making a detailed quantitative understanding of thermal transport at the micro- and nanoscales essential. To achieve this, several computational approaches are available, such as molecular dynamics [1], anharmonic lattice dynamics [2], and atomistic Green’s functions (AGF)[3]. All of these approaches require an accurate description of interatomic forces for comparably large supercells, which often need to be calculated millions of times. This makes first-principles methods like density functional theory (DFT) computationally too expensive. As a result, empirical force-field potentials (FFPs) are commonly used instead. These are parameterized analytical functions that approximate interatomic interactions. The parameters are typically obtained by fitting the empirical FFPs to experimental or first-principles data for material properties relevant to the intended application, and often, empirical FFPs are designed to work for a whole class of materials. Unfortunately, this approach does not guarantee a highly accurate description of the underlying forces, stresses, and energies. By fitting material properties, it is merely guaranteed that the configurations associated with these properties are described well. The parameterized empirical FFPs are frequently incapable of extrapolating to unfamiliar configurations. Recently, machine-learned FFPs have gained popularity as a promising alternative [4,5]. They aim to achieve close to first-principles accuracy while retaining the low computational cost of empirical FFPs. This is usually achieved by expressing the potential energy in terms of highly flexible analytical basis functions, or descriptors, thereby enabling systematic improvements in accuracy by including higher-order terms. Unlike empirical FFPs, the parameters in machine-learned FFPs are determined from ab initio datasets that include forces, stresses, and energies for relevant atomic configurations. This approach ensures that both the derived material properties and the underlying physics are accurately represented. In this work, we generate machine-learned potentials using MLIP (version 2) [6], a software package that uses moment tensors as basis functions to expand the potential energy as a sum of atomic energies. In particular, we evaluate how well the trained moment tensor potentials (MTPs) describe key material properties, including lattice and elastic constants, harmonic and anharmonic phonon properties, thermal expansion, and thermal transport. This is done for materials relevant to modern microelectronics: wurtzite AlN, GaN, and ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N. To do so, we compare MTP predictions with experimental data and first-principles results, and analyse how the choice of training-dataset hyperparameters affects model performance.

2. Methodology

2.1. MLIP

In this work, we use MLIP (Machine-Learned Interatomic Potentials) [6] to train machine-learned force field potentials of the moment tensor potential type (MTP), which we will refer to as MTPs throughout this paper. Additional details on moment-tensor potentials can be found in the Supporting Information. Two approaches for training machine-learned force field potentials are commonly used throughout the literature, passive training [7,8], which uses a fixed training dataset, and active training, where training configurations are dynamically selected, in an iterative procedure, from molecular dynamics (MD) trajectories [9,10,11]. While active learning enables a more targeted selection of relevant configurations, its implementation requires a well-structured workflow to manage tasks such as model updates, uncertainty estimation, and data selection. Recent developments are focusing on software solutions that tackle this complexity [12]. In passive training, configurations can be obtained either from molecular dynamics (MD) trajectories—e.g., via ab initio molecular dynamics [13]—or from supercells with randomly displaced atoms. In contrast to active learning, the MD-based approach is non-iterative and typically relies on single simulations, which can be computationally expensive. The randomly displaced approach, by contrast, is simpler to implement and offers a more uniform sampling of configuration space, thereby reducing the likelihood that rare or strongly correlated outlier configurations dominate or destabilise the model’s performance. Especially in the initial stages, when only a limited number of training configurations is available, this results in more stable force fields [14]. The downside is that the training configurations are not necessarily representative of actual physical conditions, making it more difficult to control certain aspects of the training data, such as the anharmonicity of the displacements, which, in MD-based methods, can be directly controlled via the temperature of the MD trajectory.

In this work, we choose the passive training approach. Our training data consist of DFT-calculated forces, stresses, and energies of random supercell configurations (Section 2.2). The parameters of our MTPs are determined by minimizing the cost function

$$\sum _{n=1}^N{W}_E{\left(E_n^{\mathrm{MTP}}\left(\overrightarrow{\xi }\right)-E_n^{DFT}\right)}^2+W_F\sum _{n=1}^N\sum _{m=1}^{M_n}{\left|{\overrightarrow{\mathrm{F}}}_{n,m}^{\mathrm{MTP}}-{\overrightarrow{\mathrm{F}}}_{nm}^{\mathrm{DFT}}\right|}^2+W_S\sum _{n=1}^N{\left|{\tilde{S}}_n^{\mathrm{MTP}}-{\tilde{S}}_n^{\mathrm{DFT}}\right|}^2$$

where $n=1,\dots ,N$ labels the configurations, $m=1,\dots ,M_n$ the number of atoms in each configuration. $E_n$ are the potential energies, ${\overrightarrow{F}}_{nm}$ the forces acting on individual atoms and ${\tilde{S}}_n$ the components of the stress tensors. For the forces, the regular vector norm is used and for the stress tensors, the Frobenius norm. The weights for forces, stresses and energies are $W_F=0.1$, $W_S=1.0$ and $W_E=1.0$. The optimization is run for 4000 steps of the underlying BFGS optimization routine [15], which is well converged for the MTP levels and dataset sizes considered here. We use a fixed number of steps because different MTP levels and training sets required noticeably different BFGS convergence thresholds to obtain reliable results, likely due to the fact that the BFGS stopping criterion is based on absolute cost values. In practice, for the cases considered here, the change in cost typically becomes negligible after about 1000–1500 steps.

2.2. Datasets

Our MTPs are trained on datasets consisting of forces, stresses, and energies obatined from DFT calculations for a variety of (distorted) atomic configurations relevant to the target system. For solid bulk materials, as in the present study, these configurations are typically supercells. In this study such configurations are obtained from supercells with randomly displaced atoms and deformed cell vectors (passive training). Specifically, we consider $2\times 2\times 2$ supercells consisting of 32 atoms, created from the wurtzite unit cells of AlN, GaN and AlGaN (determined with Vegard’s law). In these supercells all atoms are simultaneously displaced in random directions by random distances ranging from zero to a maximum displacement, ${\Delta }_d$. The equilibrium supercell vectors ${\overrightarrow{a}}_i$ are varied according to ${\overrightarrow{a}}_i^{\prime }=\overrightarrow{a}i+\alpha \left|{\overrightarrow{a}}_i\right|\overrightarrow{n}$, where $\alpha$ is a random factor between zero and a maximum deformation ${\Delta }_c$ given as a fraction of the original cell vector length, and $\overrightarrow{n}$ is a normalized random direction vector. After deforming the cell vectors, the atomic positions are shifted so that the scaled atomic positions remain unchanged. This formulation allows for fully triclinic cell deformations. Therefore, configurations are defined by two parameters: the maximum atomic displacement, ${\Delta }_d$ (in Å), and the maximum fractional change in cell vectors, ${\Delta }_c$. In practice, two datasets were considered, each consisting of $N=300$ randomly generated configurations. The primary training dataset, referred to as "set 1", was generated with a maximum atomic displacement of ${\Delta }_d=0.2$ Å and a maximum cell deformation of ${\Delta }_c=0.02$ (2%). The secondary dataset "set 2", employs smaller atomic displacements (${\Delta }_d=0.05$ Å) and smaller cell deformations ${\Delta }_c=0.005$ (0.5%). Considering these two sets allows investigating the effects of sampling different anharmonic regimes. In both cases, an additional $N=50$ random configurations were generated for validation purposes using the same parameters. Since "set 1" serves as the primary dataset for training the material-specific MTPs, separate "set 1" datasets were constructed for each material (AlN,GaN, ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N). For alloys (${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N), the configurations were evenly distributed across all alloy compositions that can be realized with $2\times 2\times 2$ supercells. Since each supercell contains 32 atoms, 16 of which are Al/Ga sites, the possible Al fractions are given by $n/16$ with $n\in \{0,\dots ,16\}$. Details of the datasets are summarized in Table 1.

2.3. Statistical Variability and Overall Performance Cost

At the start of the fitting process, MLIP initializes the MTPs’ parameters randomly, introducing considerable variability among the resulting MTPs. This variability arises because the cost function that is minimized during the fitting is highly complex, containing numerous local extrema due to the large number of fitting parameters and the substantial size of the training data set. Since the number of fitting parameters in MTPs is typically too large for global optimization methods, local optimization algorithms suitable for high-dimensional problems are used instead. In MLIP this is the BFGS routine [15]. Depending on the initial position in parameter space, the optimization converges to one of the many local minima. To account for these inherent fluctuations, we first focus on average trends and train ensembles of 10 independently parameterized MTPs. The ensembles allow us to estimate the variability of predicted results and investigate the impact of different training-related hyperparameters. Here we focus on three key hyperparameters: the MTP level, which measures how many basis functions are included in the MTP expansion (see Supporting Information for details); the number of configurations included in the training data set; and the maximum atomic displacements used in the training data. For a straightforward assessment of the general quality of an MTP, without having to analyze individual properties, we define a performance cost in the form of a second cost function, separate from the one being optimized by MLIP, that combines the deviations between MTP and DFT results for several material properties. The performance cost function is constructed as a weighted sum of the normalized root mean square errors (NRMSE) between DFT and MTP simulated material properties studied in this work. These include lattice constants (1.0), elastic constants (1.0), forces (2.0), stresses (2.0), and energies (2.0) of configurations in the validation set, phonon band structure (2.0), mode-resolved thermal conductivities (3.0), and the thermal conductivity tensor (3.0), where the values in parentheses denote dimensionless weights that control the relative importance of each property. The motivation behind these weights is given below. To combine errors associated with different material properties, which may be characterized by distinct magnitude scales and different numbers of target data points, an appropriate and scale-consistent metric is required. In this work, we employ the normalized root-mean-square error (NRMSE) [16]. As an average error measure, the NRMSE is independent of the number of data points associated with a given property. Furthermore, normalization by the full range of the corresponding target values ensures a consistent scale across different properties, enabling a balanced and meaningful comparison of errors. It is given by

$$S_y={\displaystyle \frac{\sqrt{\frac{1}{N}{\sum }_{i=1}^N{(y_i-y_{i,\mathrm{DFT}})}^2}}{\max \left(y_{\mathrm{DFT}}\right)-\min \left(y_{\mathrm{DFT}}\right)}},$$

where y denotes a material property predicted by the MTP and represented as an array, i indexes it’s elements, $y_{i,\mathrm{DFT}}$ are the corresponding DFT reference values, N is the number of data points, and min and max refer to the minimum and maximum values of the DFT results. Another aspect that needs to be taken into account is that the supposedly more accurate MTPs with higher levels (and, thus, more fitting parameters) are significantly slower than those with lower levels. To what extent this is relevant depends strongly on the type of simulations to be performed with the MTPs: for example, for calculating harmonic phonon frequencies the number of simulations to be performed as well as the computational effort per simulation are much smaller than when calculating molecular dynamics trajectories for large systems over many nanoseconds. Therefore, for the overall performance it is useful to also evaluate the impact of the MTP level on computational efficiency. To this end, we determined the average time per MD step for NPT simulations of 50 ps duration (using a 0.5 fs timestep) perfomed on $10\times 10\times 10$ supercells. The weighting factors are not uniquely defined and therefore introduce a degree of arbitrariness. In the present work, they were selected to reflect the authors’ assessment of the relative importance of accurately reproducing the respective physical properties. For instance, lattice constants are well described by all MTPs and are, therefore, assigned a relatively low weight of 1.0. Similarly, elastic constants are given a weight of 1.0, as they are less critical for the intended application of thermal transport simulations. Harmonic phonon frequencies, forces, stresses, and energies, while important, are comparatively easier to reproduce and are assigned a slightly higher weight of 2.0. In contrast, mode-resolved and bulk thermal conductivities, which are central to the application, receive the highest weight of 3.0.

2.4. MTP-Based Material Property Calculations

To assess the MTPs capability to accurately capture the phyiscs in wurtzite AlN, GaN and their alloys, structural, harmonic and quasi-harmonic vibrational, and anharmonic thermal transport properties are benchmarked. These include lattice and elastic constants, the volume thermal expansion, the phonon band structure, the Grüneisen parameters and the thermal conductivities. All MTP-based calculations are performed using LAMMPS (version 20220623) [17], an open-source molecular dynamics simulation package. Due to the complexity of evaluating material properties of the alloy ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N, we restrict the benchmarking to just ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N and model it as a periodic system whose primitive unit cell contains one aluminum and one gallium atom.

The static equilibrium lattice constants are obtained with LAMMPs, using the minimize command together with fix box/relax. This performs a structural optimization, which simultaneously minimizes the total energy and stress tensor. The optimization is terminated when the norm of the global force vector, a composite vector of the forces acting on all atoms, falls below ${10}^{-10}$ eV/Å. We also report the temperature-dependent volume and volume thermal expansion coefficients predicted by our MTPs. These are determined from molecular dynamics simulations and from lattice dynamics using the quasi-harmonic approximation. For the MD results, we run NPT simulations on $10\times 10\times 10$ supercells relaxed for 500 ps using a Nose-Hoover style barostat. A timestep of ${\mathrm{t}}_{\mathrm{step}}=0.5$ fs was used and the temperature and pressure damping parameters were ${\mathrm{T}}_{\mathrm{damp}}=0.05$ ps and ${\mathrm{P}}_{\mathrm{damp}}=0.5$ ps, respectively. The first 200 ps were discarded as they are associated with the equilibration of the system. The volume thermal expansion coefficient ${\alpha }_{\mathrm{V}}\left(\mathrm{T}\right)={\displaystyle \frac{1}{\mathrm{V}}}{\displaystyle \frac{\mathrm{dV}}{\mathrm{dT}}}$ is obtained by fitting the volumes calculated at different temperatures with a model function and taking the derivative with respect to the temperature. For classical results, such as those produced by MD, we use a cubic polynomial ansatz and for experimental results, which include quantum effects, we use an ansatz based on the Einstein model [18]. They are given by

$$\begin{array}{cc}\hfill \mathrm{V}\left(\mathrm{T}\right)& =c_0+c_1\mathrm{T}+c_2{\mathrm{T}}^2+c_3{\mathrm{T}}^3\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \mathrm{V}\left(\mathrm{T}\right)& =c_0+c_1\left(1+{\displaystyle \frac{2}{exp\left(\frac{c_2}{\mathrm{T}}\right)-1}}\right),\hfill \end{array}$$

(2)

where $c_i$ are the fit parameters. The second approach for calculating the thermal expansion is the quasi-harmonic approximation (QHA), which is implemented in Phonopy [19], an open-source package for phonon calculations. In this approach, the phonon frequencies are allowed to vary with volume. This allows one to determine the Gibbs free energy to extract the temperature dependence of several key material properties. These include the equilibrium volume $V\left(T\right)$, the corresponding volumetric thermal expansion coefficient ${\alpha }_v\left(T\right)$, the bulk Grüneisen parameter $\gamma \left(T\right)$ and the bulk modulus $B\left(T\right)$. The approach requires second-order force constants to be computed at multiple volumes. In this study, we use a set of 11 volumes obtained by isotropically straining the cell vectors (the same factor for all lattice constants) ranging from -2.5% to +2.5% around the static equilibrium volume.

The elastic constants describe the relationship between applied stress and the resulting strain. They are directly related to the propagation of sound waves, and, in the long-wavelength limit, they determine the propagation of acoustic phonons. The MTP-based elastic constants were calculated using LAMMPS, using the example script provided on the official LAMMPS website. This approach is based on the linear stress–strain relationship, $s_{ij}=C_{ijkl}{ϵ}_{kl}$, where $C_{ijkl}$ are the elastic constants, and all indices run over the three Cartesian directions. The elastic constants are obtained by evaluating the derivatives of the stress tensor components $s_{ij}$ with respect to the strain components ${ϵ}_{kl}$. In practice, unit cell vectors are first deformed by ${10}^{-4}$%, after which the atomic positions are relaxed until the norm of the total force vector acting on all atoms falls below ${10}^{-10}$ eV/Å. The resulting stress tensors are then computed with LAMMPs to extract the elastic constants.

Phonon calculations were performed with Phonopy [20], employing a $5\times 5\times 5$ supercells for second-order force constants (FC2) and a $20\times 20\times 20$ Brillouin zone grid for sampling the phonon reciprocal space. The FC2 are computed using the finite-difference method implemented in Phonopy, with a default atomic displacement amplitude of 0.01 Å.

The mode-Grüneisen parameters ${\gamma }_{i\nu }$ provide a quantitative measure of the anharmonicity of the interatomic potential and are related to thermal transport properties. They are defined as the (negative) relative changes in phonon frequencies as a function of relative changes in unit cell volumes. Here we determine the mode-Grüneisen parameters with ShengBTE, the software used for calculating the thermal conductivity (see below). ShengBTE uses a different expression that yields identical results but links them to the third-order force constants (FC3), which more directly relate them to thermal transport. Details can be found in [21]. The results provided by the FC3-based approach are usually more accurate for MTPs, since the FC3 are calculated from equilibrium supercells with small atomic displacements, which are far easier to predict than the far-from equilibrium, isotropically strained configurations needed for the first approach. This is particularly relevant for these highly harmonic systems, where significant deviations from equilibrium configurations would be necessary to properly capture the variation of phonon frequencies with respect to volume.

Lastly, the thermal conductivity is calculated using ShengBTE [21], a program based on anharmonic lattice dynamics (ALD). In this perturbative approach, the phonon Boltzmann transport equation (BTE) is solved using phonon scattering rates obtained from perturbation theory. Anharmonic phonon–phonon interactions are included through higher-order terms in the expansion of the lattice potential energy. In practice, this expansion is typically truncated at third-order (cubic) terms due to the immense complexity and computational cost associated with evaluating higher-order scattering processes [22], and ShengBTE follows this standard approximation. As a consequence of neglecting higher-order anharmonicity, the method becomes less accurate at elevated temperatures and in materials with strongly anharmonic interactions. To assess the influence of these higher-order effects, we also performed calculations with the FourPhonon module [22], an extension to ShengBTE that includes fourth-order terms. Because these calculations are computationally demanding, we restricted them to two temperature points (300 K and 800 K) for GaN, which is sufficient to illustrate the impact of higher-order anharmonicity. The required third (FC3) and fourth (FC4) order force constants were computed with the finite-displacement method using Python scripts provided by ShengBTE and FourPhonon. For FC3 we used $2\times 2\times 2$ supercells and the default atomic displacement of $0.01\text{Å}$. The forces in these supercells were obtained with LAMMPS [17]. The Peierls–Boltzmann equation is solved within the relaxation-time approximation on a $20\times 20\times 20$ reciprocal-space phonon grid in the case of ShengBTE and $12\times 12\times 12$ grid for FourPhonon. We observed, that for the materials studied here, the change between RTA and a full iterative solution is small $<5\%$ at room temperature and above, therefore we restrict ourself to the RTA. The convergence with respect to the grid size is shown in Figure S7a (3-phonon) and Figure S7b (4-phonon). For all materials, isotope scattering was included for a better comparison with experimental data. This is only really relevant for GaN, where the impact is stronger [23]. For thermal transport simulations of the alloy ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N, instead of treating it as a hypothetical ordered compound as done for the other observables (see beginning of this section), we employed the virtual crystal approximation (VCA). In this approach, the lattice parameters, atomic masses, and second- and third-order force constants (FC2 and FC3) of AlN and GaN were linearly interpolated according to the alloy fraction x. The effects of mass disorder were additionally included using the perturbative Tamura model [24].

Finally, to evaluate the contributions of different spectral regions to the thermal conductivity, we calculated the average frequency-resolved thermal conductivity ${\langle \kappa \rangle }_{\omega }$. Since the mode-thermal conductivities are calculated on a q-point grid this quantity is approximated in pratice by Gaussian smearing

$${\langle \kappa \rangle }_{\omega }={\displaystyle \frac{{\sum }_{\nu }{\sum }_i{\kappa }_{i\nu }\delta \left(\omega -{\omega }_{i\nu }\right)}{{\sum }_{\nu }{\sum }_i\delta \left(\omega -{\omega }_{i\nu }\right)}}={\displaystyle \frac{{\sum }_{\nu }{\sum }_i{\kappa }_{i\nu }\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{{\left(\omega -{\omega }_{i\nu }\right)}^2}{2{\sigma }^2}\right)}{{\sum }_{\nu }{\sum }_i\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{{\left(\omega -{\omega }_{i\nu }\right)}^2}{2{\sigma }^2}\right)}},$$

(3)

where i and $\nu$ label the reciprocal grid points and phonon branch indices, ${\kappa }_{i\nu }$ are the mode-thermal conductivities, ${\omega }_{i\nu }$ are the phonon frequencies on the Brillouin zone grid and $\sigma$ is the standard deviation assumed for the smearing. We use a standard deviation of half the frequency resolution of $\omega$ for smearing.

2.5. First-Principles Calculations

All DFT calculations were performed with Quantum ESPRESSO [25], using the PBESol exchange–correlation functional along with pseudopotentials of the precision version from the Standard Solid-State Pseudopotential (SSSP) library [26]. The wavefunction and density cutoffs used are those suggested by the authors of the pseudopotential library, which we independently verified to be sufficiently converged. Convergence tests can be found in the supplementary material (Section S6). The wavefunction cutoffs are 30 Ry (408 eV) for Aluminium, 90 Ry (1225 eV) for Ga and 80 Ry (1088 eV) for Nitrogen. In practice, the largest cutoff among the elements present in the system is used: 80 Ry (1088 eV) for AlN and 90 Ry (1225 eV) for both GaN and ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N. The density cutoffs used are 8 times the wavefunction cutoff. All DFT calculations use a self-consistency threshold for the electronic energy of 1e-12 Ry. The reciprocal-space grid size for electrons is application-specific and reported below.

The first-principles lattice constants were calculated with QE using the vc-relax option, minimizing both energy and stresses by adjusting atomic coordinates and cell vectors. The calculation used a 12×12×12 k-point grid for self-consistently calculating the electron density and convergence thresholds of ${10}^{-6}$ Ry ($1.36·{10}^{-5}$eV) and ${10}^{-6}$ Ry/Bohr ($2.57·{10}^{-5}$ eV/Å) for the total energy (etot_conv_thr) and forces (forc_conv_thr) were chosen.

Elastic constants were obtained from the QE driver package thermo_pw [27]. The calculation was performed on a 12×12×12 reciprocal-space grid.

DFT-based second-order force constants were obtained using Quantum ESPRESSO through density functional perturbation theory (DFPT), where the dynamical matrices are computed on a 5×5×5 phonon reciprocal-space grid and subsequently Fourier transformed into real space. The electronic structure calculations performed in the course of the DFPT simulations employed a 10×10×10 k-point grid and a phonon self-consistency threshold (tr2_ph) of ${10}^{-16}$.

Third-order force constants (FC3) were computed for 2×2×2 supercells, using the finite displacement method implemented in a python script provided by ShengBTE. The required forces are calculated with QE on a 4×4×4 electron k-point grid, which combined with the real space supercell size is sufficient for convergence. The reference supercell configurations for the training and validation of MTPs are computed with the same settings as the FC3 supercells.

3. Results

3.1. Statistics and Trends

In this section, we discuss general observations concerning the influence of the MTP level, the number of configurations in the training datasets, and the magnitude of atomic displacements and cell deformations in those training configurations. The effect of the MTP level is relatively straightforward and is illustrated in the top panel of Figure 1. It shows the performance cost, evaluated according to the procedure described in Section 2.3, for 10 MTPs at levels 10, 12, 14, 16, and 18. As expected, the performance cost decreases with increasing MTP level, demonstrating improved overall accuracy with MTP level. The variance of cost values within each ensemble also genereally tends to decrease with MTP level. Because the MTP basis at higher levels provides a better representation of the true potential energy surface, the fitting error is dominated by approximation bias rather than statistical variance when an abundance of data is available. As a result, at higher MTP level independently trained models converge toward similar solutions and thus, exhibit reduced ensemble variability of predicted material properties. This behavior reflects the well-known bias–variance decomposition [28], whereby increasing model capacity primarily reduces bias when sufficient training data constrain the parameters.

Despite the continous improvement with MTP level, the best performaning MTP at level 18, achieves only a $\sim 20\%$ lower performance cost than the best performing MTP at level 10, presumably due to the comparably simple structure of the inorganic semiconductor materials investigated here. At the same time, the computational cost, illustrated in the center panel of Figure 1, increases by roughly a factor of seven. A practical approach is, therefore, to tailor the MTP to the specific application, as some of us did for conjugated polymers [29,30]. If an application requires very high accuracy, i.e., lattice dynamics or other non-MD calculations, a higher MTP level is justified. In this work, however, we aim for a general-purpose potential that can also be used efficiently in large-scale MD simulations, and therefore mainly focus on MTPs of level 10.

The bottom panel of Figure 1 shows how the number of training configurations affects the accuracy of the results. Increasing the number of training configurations from 100 to 300 clearly improves the overall performance of the MTPs.

To obtain a more comprehensive picture, we next briefly discuss the impact of these hyperparameters, for some of the individual material properties studied in this work. In general, all properties improve with higher MTP level. This can be clearly seen in Figure S1.

The average increase in accuracy is substantial, generally ranging between $\sim 30-60\%$ depending on the property, with phonon frequencies and mode-thermal conductivities on the lower end (clsoer to ∼30%) and forces, stresses, energies, lattice constants and elastic constants on the higher end (closer to $\sim 60\%$). However, due to the huge variability among MTPs of level 10, especially for derived properties such as the lattice constants, elastic constants, phonon frequencies and mode-thermal conductivities, the variations between individual MTPs of the same level clearly exceeds the overall increase in accuracy. The differences in errors between the best performing level-10 MTP and those of the average level-18 MTP is generally much smaller. For example, for the phonon-related properties such as the frequencies or mode-thermal conductivity the difference is $8\%$ and $4\%$ respectively. This justifies our approach, especially when it comes to phonon-related computations.

Notably, as far as the lattice constant are concerned, all trained MTPs were able to reproduce the DFT results remarkably well, with mean absolute errors (MAEs) of approximately ${10}^{-3}\text{Å}$. The average error for DFT drops with the MTP level up to MTP level 14, at which point it remains roughly the same (Figure S1d), indicating that a maximum accuracy for the given training dataset has been achieved. Increasing the number of training configurations from 100 to 300 also does not improve the lattice constants (Figure S2d) and even slightly worsens the results. This suggests that MTPs of level 10, in this particular case, have reached their representational limit with respect to the number of training data, and the remaining error is likely due to limited model flexibility rather than insufficient data, as shown in Figure S2d. The MTPs trained on the dataset with smaller maximum atomic displacements (set 2, Table 1) achieve better results for the lattice constants (see Figure S3d). This can be explained by a denser sampling of near-equilibrium configurations due to the smaller strain, allowing the level-10 MTP to more accurately capture the equilibrium geometry.

The elastic constants, unlike the lattice constants, improve significantly with increasing dataset size (see Figure S2e) because they depend on the second derivative, that is, the curvature of the potential energy surface with respect to strain, whereas the lattice constants are determined by the first derivative. As a result, accurately describing the elastic constants requires a more thorough sampling of configurations around equilibrium. The effect of cell deformations on the elastic constants is twofold. First, a minimum amount of strain must be present in the training data for the MTP to interpolate stresses reliably; however, this threshold is extremely small ($\ll {10}^{-3}$) and was comfortably exceeded in all datasets, becoming noticeable only when no strain at all was included. Second, when the configurations are too strongly strained, the sampling density near equilibrium decreases, which, similar to the case of the lattice constants, can lead to less accurate predictions (Figure S3e). A similar observation has been made by Reicht et al. [29], when training MTPs for describing polymers using an active-learning approach: also there for elastic constants (and phonon-related properties), a better agreement with DFT data was observed when training the MTPs on low-displacement structures (in that case generated via low-temperature MD trajectories). Another important factor to consider for the hexagonal cells of AlN, GaN, and AlGaN is that the the training data contains cells that are sufficiently triclinically strained. Otherwise, components of the elastic constants other than $C_{11}$ and $C_{33}$ are poorly described.

Phonon frequencies benefit slightly from larger datasets sizes (Figure S2f) and become slightly worse with larger displacements (Figure S3f). The reasons are the same already discussed for the elastic constants. Anharmonic (FC3-based) phonon properties, such as the mode-Grüneisen parameters (Figure S2g) and mode-thermal conductivity (Figure S2h) show slight improvements with larger dataset sizes. The Grüneisen parameters (FC3-based) also improve with larger atomic displacements (Figure S3g), while the mode-thermal conductivities (Figure S3h) remain mostly unchanged.

In summary, for our purposes, MTPs of level 10 provide the best compromise between performance and accuracy and will therefore be the focus for the remainder of this work.

So far, we have analyzed the performance of the MTPs in terms of the overall effective performance cost, defined as the sum of the weighted, normalized mean error contributions. We have also examined the trends in the errors of individual properties with respect to the training-related hyperparameters. In both cases, the analysis focused on ensemble-averaged properties of the trained MTPs. Because most practical applications require the use of a single force field throughout a calculation, and because reporting benchmarking results for all MTPs in the ensemble across the three materials and all considered properties would be unnecessarily exhaustive while providing limited additional insight, a representative force field is selected from the main ensemble for each material based on its overall performance (see Figure 2). The performance cost implicitly reflects our prioritization of specific material properties through the chosen weights (see Section 2.3). As the present work focuses on thermal transport, the mode-thermal conductivity is assigned the highest weight (3.0) followed by the phonon band structure (2.0). Consequently, the selected MTPs exhibit the lowest, or near-lowest (where differences are small), errors for these properties. For properties assigned lower weights (e.g., lattice and elastic constants), the errors are generally at least below average in all cases.

In the following section, we present a detailed comparison of the predictions obtained from the selected MTPs with results from DFT calculations and available experimental data for all considered observables. Thus, in the remainder of this work when we refer to our MTPs, we explicitly mean the selected optimal level-10 MTPs for each material.

3.2. Prediction of Forces, Stresses and Energies

First, we evaluate the MTPs’ performance on the validation set. Figure 3 presents the mean absolute errors (MAEs) for forces, stresses, and energies. All MAEs fall below 40 meV/Å for forces, 0.15 GPa for stresses, and 1.4 meV/atom for energies. These values are comparable to those reported for state-of-the-art machine-learned FFPs for semiconductor materials. For example, Minamitani et al. [31] train high-dimensional neural network potentials and achieve a MAE for the forces of $25.5$ meV/Å for Si and $37.8$ meV/Å for GaN and Gao et al. [32] use a Deep Potential Model on phosphorus-doped silicon and report a MAE of 33 meV/Å. In our case, depending on the MTP and material, errors can be even lower. For example, in GaN, the MAEs for forces, stresses, and energies are around 5 meV/Å, 0.03 GPa, and 0.3 meV/atom, clearly outperforming current standards. The larger force errors observed for AlN compared to GaN may reflect it’s stiffer, more ionic Al–N bonds and the resulting steeper potential-energy surface, which tend to produce larger forces for similar atomic displacements. These factors make it more challenging for FFPs to accurately reproduce forces, stresses and energies, at least at the fairly low MTP levels considered here.

As expected, the material-specific MTPs perform better because they are trained exclusively on AlN or on GaN configurations and do not need to interpolate, regarding the material’s composition. In contrast, the alloy MTPs are trained on configurations spread uniformly across all alloy fractions, leaving relatively few data points near the ends of the composition range. Despite this, our alloy MTP still achieves highly competitive results.

In addition to reporting average deviations, it is useful to compile correlation plots that illustrate the resemblance between MTP and DFT results for all calculated quantities. Correlation plots for forces, energies and strains are shown in Figure 4. Such plots would reveal systematic biases or scaling issues. In correlation plots, a nonzero intercept of a hypothetical line through the datapoints would indicate a systematic offset (bias), while a slope of that line different from one would reveal a scaling issue with the predicted values being proportionally too large or too small. To illustrate that such problems do not occur here, all plots in Figure 4 contain straight lines through the origin with a slope of exactly one (red lines). To be more quantitative, we also calculated the Pearson correlation coefficient (R value), which measures the linearity and agreement between the two data sets, with values above 0.99 indicating excellent correspondence. Our MTPs achieve R values of 0.999 and higher, thus providing an essentially DFT-level accurate and reliable description of forces, stresses, and energies.

Figure 3. The figure shows the mean absolute errors in forces, stresses, and energies predicted by the selected MTPs for the validation datasets of AlN, GaN and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N. Results for the force fields fit to the individual material are shown in red and those for the alloy force field in green.

Figure 3. The figure shows the mean absolute errors in forces, stresses, and energies predicted by the selected MTPs for the validation datasets of AlN, GaN and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N. Results for the force fields fit to the individual material are shown in red and those for the alloy force field in green.

Figure 4. Correlation between DFT and MTP-calculated forces, stresses, and energies for structures taken form the validation sets of AlN. The panels on the left show results for the material-specific MTP and those on the right show results for the alloy MTP. The red lines represent the linear regression fit between the DFT and MTP-predicted data.

Figure 4. Correlation between DFT and MTP-calculated forces, stresses, and energies for structures taken form the validation sets of AlN. The panels on the left show results for the material-specific MTP and those on the right show results for the alloy MTP. The red lines represent the linear regression fit between the DFT and MTP-predicted data.

Figure 5. Correlation between DFT and MTP-calculated forces, stresses, and energies for structures taken form the validation sets of GaN. The panels on the left show results for the material-specific MTP and those on the right show results for the alloy MTP. The red lines represent the linear regression fit between the DFT and MTP-predicted data.

Figure 5. Correlation between DFT and MTP-calculated forces, stresses, and energies for structures taken form the validation sets of GaN. The panels on the left show results for the material-specific MTP and those on the right show results for the alloy MTP. The red lines represent the linear regression fit between the DFT and MTP-predicted data.

3.3. Lattice Constants

In our study, we focus on the wurtzite phase of the materials investigated, which has a hexagonal crystal structure with space group ${\mathrm{P6}}_3$mc. The primitive cell is defined by two lattice parameters, a and c. The DFT- and MTP-calculated static values of these lattice constants, as well as the results of low-temperature diffraction experiments are listed in Table 3. For ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N, experimental data are only available for thin films, which are strained and, therefore, deviate from the bulk values [33]. Instead, we apply Vegard’s law [34] to estimate the alloy lattice constants from the average of the AlN and GaN values listed in Table 2. The agreement between the MTP and DFT simulations is excellent, with lattice constants agreeing to within four significant digits. For the material-specific MTPs, the mean relative deviations between MTP and DFT results are 0.01%, 0.01%, and 0.08% for AlN, GaN, and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N, respectively. The alloy MTP also excellently reproduces the AlN and GaN lattice constants, with slightly larger errors of 0.09% and 0.03%, respectively. This shows that, for lattice constants, the MTP perfectly reproduces the results of the parent methodology used to train it.

Both the DFT and MTP results obtained by static geometry optimizations agree extraordinarily well with the low-temperature experimental data. This agreement, however, must be taken with a grain of salt, as in the simulations, phonon effects and, in particular, the zero-point motion of the atoms, have been neglected in the structural optimization. When phonon contributions are included using the quasi-harmonic approximation (QHA), the lattice constants increase by $\sim 0.25\%$ for all materials (see DFT${}^{}\mathrm{anh}$, MTP${}^{}\mathrm{anh}$, MTP${}^{}\mathrm{anh}_{\mathrm{alloy}}$ in Table II). As a consequence, the agreement between theory and simulations becomes slightly worse. Notably, because MTP and DFT results show very similar changes, this is not due to the MTP parametrisation but rather an intrinsic shortcoming of the underlying DFT methodology. The slight overestimation of the lattice constants when zero-point motion is included in the simulations is most likely due to the PBESol exchange-correlation functional used in the calculations. Functionals based on the generalized gradient approximation (GGA), such as PBE, are known to underbind, which generally leads to overestimated lattice constants [35]. The PBESol functional was developed to compensate for the underbinding tendency and performs well for many materials [36]. In our case, when disregarding the zero-point motion, a partial error cancellation between the remaining tendency to underbind and the neglected zero-point motion occurs. This error cancellation is eliminated by explicitly considering the zero-point motion, and the lattice constants are overestimated. While this complication leads to an offset in the “0K” lattice constants, it is not expected to impact the thermal expansion of the materials.

This is indeed the case, as illustrated in Figure 6, which shows the temperature-dependent volumes for AlN (top-left), ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N (top-center), and GaN (top-right). The corresponding volume thermal expansion coefficients are displayed in the bottom panels. The figure presents molecular dynamics (MD) results obtained with the material-specific MTPs (red data points and lines) and the alloy MTP (green data points and lines). Additionally, results obtained using the quasi-harmonic approximation (QHA, see Section 2.4) are shown for the MTPs (red and green dashed lines) and for DFT (blue dashed lines). Most importantly, the QHA results reveal an excellent agreement between the simulations performed with DFT and those using an MTP. This is particularly true for the system-specific force fields but also the alloy MTP, which provides a satisfactory agreement but somewhat overestimates the crystal volume close to 0 K and slightly underestimates the thermal expansion at higher temperatures. At elevated temperatures the MD results calculated with the MTPs also agree well with the DFT-based QHA predictions. This agreement, however, deteriorates at lower temperatures, which is not surprising considering that in the MD simulations, the effective phonon calculations follow the equipartition theorem [37], neglecting the quantum effects manifested in the Bose-Einstein statistics. This is consistent with the results of Ramírez et al. [18], who observed a similar effect in silicon carbide (SiC) when comparing classical molecular dynamics simulations with path-integral molecular dynamics.

Regarding a direct comparison with experimental measurements (unconnected datapoints in Figure 6), one observes the above-discussed overestimation of the low-temperature unit cell volume. Nevertheless (as already suggested above), the change in volume with temperature and, thus, the volumetric thermal expansion coefficients are excellently reproduced by all QHA simulations. More specifically, for AlN and GaN, the QHA results obtained with the material-specific MTPs (red dashed line) and DFT (blue dashed line) accurately reproduce experimental measurements, and the agreement only starts to deteriorate for GaN above approximately 400 K, where the QHA somewhat overestimates the thermal expansion. We attribute this to additional anharmonic effects that are not captured within the QHA framework and that are more pronounced in GaN, which is more anharmonic in nature, as can be seen by the lower thermal conductivity, which will be discussed in section Section 3.7, and the smaller phonon life times [48].

3.4. Elastic Constants

Table 3 presents the elastic constants calculated using DFT and the MTPs, together with the range of experimental values found in the literature. The full collection of experimental results from the literature is given in Table S1. Overall, the elastic constants predicted by our MTPs agree well with the DFT results, with average relative errors of only a few percent across all investigated materials. For the material-specific MTPs, when considering all five independent elastic constants, we obtain mean relative errors of $2.4\%$, $0.5\%$, and $3.6\%$ for AlN, GaN, and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N, respectively. Applying the alloy MTP to AlN and GaN yields mean errors of $6.3\%$ and $2.7\%$. Our calculated values also show good agreement with experimental data, with relative errors of $5.6\%$ and $4.8\%$ for the DFT results, $7.5\%$ and $5.5\%$ for the material-specific MTPs, and $12\%$ and $3.9\%$ for the alloy MTP. It is important to note that the experimental measurements were conducted at T = 300 K, while the calculated values correspond to T = 0 K.This is, however, not expected to introduce significant discrepancies, as temperature-dependent measurements for AlN by Kim et al. [28] and for GaN by Adachi et al. [29] indicate that the elastic constants are only weakly temperature-dependent. A general trend is that compared to experiments, both our DFT and MTP results underestimate the elastic constants ${\mathrm{C}}_{11}$ and ${\mathrm{C}}_{33}$. This can be attributed to the underbinding caused by the chosen exchange–correlation functional (see Section 3.3), which results in an underestimation of the stiffness of the bonds and smaller elastic constants.

Table 3. DFT and MTP calculated static elastic constants of wurtzite AlN, GaN and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N. Also shown is the range of experimental values at room temperature found in the literature (the full table of experimental data can be found in S1). MTP and ${\mathrm{MTP}}_{\mathrm{alloy}}$ dennote the MTPs trained specifically for the material and for the ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N alloy, respectively. Units are given in GPa.

	AlN				GaN				${\mathrm{Al}}_{\mathbf{0}.\mathbf{5}}{\mathrm{Ga}}_{\mathbf{0}.\mathbf{5}}$N
	DFT	MTP	${\mathrm{MTP}}_{\mathrm{alloy}}$	Exp.	DFT	MTP	${\mathrm{MTP}}_{\mathrm{alloy}}$	Exp.	DFT	${\mathrm{MTP}}_{\mathrm{alloy}}$
${\mathbf{C}}_{\mathbf{11}}$	381	382	380	401−413	346	342	346	365−390	363	362
${\mathbf{C}}_{\mathbf{33}}$	356	360	389	368−390	384	385	387	379−398	381	385
${\mathbf{C}}_{\mathbf{44}}$	112	111	103	120−127	92	92	89	90−109	98	100
${\mathbf{C}}_{\mathbf{12}}$	137	144	145	127−149	128	128	133	106−145	133	139
${\mathbf{C}}_{\mathbf{13}}$	108	113	109	96−119	93	96	97	70−114	98	108

Table 3. DFT and MTP calculated static elastic constants of wurtzite AlN, GaN and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N. Also shown is the range of experimental values at room temperature found in the literature (the full table of experimental data can be found in S1). MTP and ${\mathrm{MTP}}_{\mathrm{alloy}}$ dennote the MTPs trained specifically for the material and for the ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N alloy, respectively. Units are given in GPa.

	AlN				GaN				${\mathrm{Al}}_{\mathbf{0}.\mathbf{5}}{\mathrm{Ga}}_{\mathbf{0}.\mathbf{5}}$N
	DFT	MTP	${\mathrm{MTP}}_{\mathrm{alloy}}$	Exp.	DFT	MTP	${\mathrm{MTP}}_{\mathrm{alloy}}$	Exp.	DFT	${\mathrm{MTP}}_{\mathrm{alloy}}$
${\mathbf{C}}_{\mathbf{11}}$	381	382	380	401−413	346	342	346	365−390	363	362
${\mathbf{C}}_{\mathbf{33}}$	356	360	389	368−390	384	385	387	379−398	381	385
${\mathbf{C}}_{\mathbf{44}}$	112	111	103	120−127	92	92	89	90−109	98	100
${\mathbf{C}}_{\mathbf{12}}$	137	144	145	127−149	128	128	133	106−145	133	139
${\mathbf{C}}_{\mathbf{13}}$	108	113	109	96−119	93	96	97	70−114	98	108

3.5. Phonon Band Structure

Figure 7 shows the phonon frequencies along the high symmetry directions in the first Brillouin zone for AlN (left), ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N (center), and GaN (right). The DFT- and MTP-calculated frequencies show very good agreement, with only minor deviations in the higher optical branches. As expected, the alloy MTP performs slightly worse than the material-specific MTPs, however, the difference is small, and its overall performance remains very good. The minor discrepancies between DFT and MTPs at higher frequencies are attributed to the presence of long-range Coulomb interactions in semiconductors with polar bonds, which must be treated separately because their dipole–dipole nature leads to non-analytic contributions to the dynamical matrix in the long-wavelength limit [58]. In DFT phonon calculations, these long-range forces are excluded from the short-range interatomic force constants and incorporated instead through the non-analytic term correction (NAC), ensuring an accurate description of the LO–TO splitting in polar materials. Unfortunately, the same approach cannot be applied to the MTPs, since long-range Coulomb forces are already, at least in part, implicitly included in the training dataset forces and, consequently, in any forces predicted by the MTPs. Thus, the long-range forces would have to be removed from the training data beforehand to avoid double-counting and an oscillatory behaviour near the Gamma point when the NAC is applied. However, this is beyond the scope of this work and is not relevant to thermal transport (the final target of the MTP parametrizations). This is because deviations occur in the upper optical branches that usually don’t contribute much. Therefore, in this work NACs are not considered in any of our calculations. Despite this limitation, the mean relative deviation between DFT-and MTP-calculated phonon frequencies, averaged over the entire, homogeneously sampled Brillouin zone, is in the range of only $1\%$.

Our data also agree very well with the experimental data from Schwoerer et al. (b)[56] for AlN (Figure 7 left, black circles) and Ruf et al. [57] for GaN (Figure 7 center, black circles). Only for the LO-TO splitting near the $\Gamma$-point, larger deviations are observed for the reasons discussed above.

In addition to the comparison between theory and experiment discussed above, Tables S2 and S3 in the Supplementary Material present a comparison with available Raman measurements. Good agreement is also observed in this case, with mean relative errors (MRE) in the frequencies of approximately 2% for all materials. The Raman data further reveal a systematic tendency of both the MTP and DFT results to slightly underestimate the experimental frequencies. This trend is consistent with the underbinding behavior discussed in Section 3.3. In principle, the observed frequency shift could also be attributed to a temperature-dependent phonon renormalization; however, previous studies [59,60] have shown that such effects are small ($<3{\mathrm{cm}}^{-1}$) and, thus, insufficient to explain the discrepancy.

3.6. Mode-Grüneisen Parameters

The mode-Grüneisen parameters characterize the relative shift in phonon frequencies with respect to changes in the unit cell volume (see Section 2.4). They are often used as a measure of a material’s anharmonicity, making them a useful quantity for benchmarking in thermal transport applications. These parameters are among the most challenging material properties to reproduce accurately, as the changes in phonon frequencies with volume are often comparable in magnitude to the errors in the MTP-predicted phonon frequencies. Despite these challenges, our MTPs reproduce the DFT reference data well, as shown in Figure 8. Some deviations persist, especially in the acoustic and lower-optical branches, which correspond to the lower-frequency region of the two. Relative errors of the material-specific MTPs with respect to DFT are 15.9%, 13.0% and 20.9% for AlN, GaN and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N. For alloy MTP applied to AlN and GaN, the errors are 14.3% and 26.0%, respectively. It should be noted that relative errors are dominated by modes with small Grüneisen parameters, where they tend to blow up. The mean absolute errors (MAE) are actually fairly small with 0.048, 0.041 and 0.049 for AlN, GaN and AlGaN in the case of the material-specific MTP, and 0.048 and 0.074 for the alloy MTP applied to AlN and GaN. When compared with the average Grüneisen parameters of the respective materials, these correspond to errors in the range of roughly $4\%$ to $8\%$. Experimental data are available only for the Raman-active modes (Table S4). Table S5 provides an overview of the DFT- and MTP-calculated mode Grüneisen parameters for the Raman-active modes, along with the range of experimental values reported in the literature. The agreement is acceptable, with mean relative errors with respect to experiments of approximately $4\%$ for DFT, $6\%$ for the material-specific MTPs, and $12\%$ for alloy MTP.

3.7. Lattice Thermal Conductivity

The final goal of our MTPs is to accurately simulate thermal transport, which requires them to reproduce both the full thermal conductivity tensor as well as the frequency-dependent thermal conductivity contributions. The focus below will be on calculating the thermal conductivity from harmonic phonon frequencies and anharmonic lifetimes, using the Boltzmann transport equation in the relaxation-time approximation. Molecular-dynamics based studies will be left for follow-up studies. Figure 9 shows the temperature-dependent values of the thermal conductivities of AlN, GaN, and ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N in the c-direction of the crystals, ${\kappa }_{zz}$. These results are calculated using anharmonic force constants obtained with DFT and our MTPs. The simulations are compared to a series of experimental data. The focus here is on heat transport in the c-direction because, in most applications, thin films are grown in the c-orientation. For the sake of completeness, Table 4 also reports the diagonal components of the thermal conductivity tensors at 300K. When comparing the solid red lines obtained from system-specific MTPs to the blue dashed lines calculated using DFT, an impressive agreement is revealed between the ab initio results and those generated by the machine-learned force fields. This testifies to the enormous potential of MTPs of the type discussed here for accurately and efficiently predicting heat transport in inorganic semiconductors. The performance of the alloy MTP is slightly worse but still acceptable.

Comparing the simulations to the experiments on single crystals (adopted from literature), one observes again an excellent agreement for AlN up to 400K especially when considering the data by Inyushkin et al. [62]. Unfortunately, there are no experimental data beyond 400 K available from that source. When comparing the simulations to the data by Slack et al. [63] a similarly favourable picture evolves at low temperatures, but one also sees that at 600K the simulations overshoot the experiments by $\sim 25\%$, with the deviation becoming worse at 1000 K. For GaN the agreement between theory and experiments at low temperatures is again satisfactory, but here already at 300K the simulations overestimate the experiments by about 10%. The deviation gets significantly worse at higher temperatures. The temperature-dependence of the deviation, as well as the fact that they are worse in the more anharmonic GaN (see above) suggest that the deviations are likely a consequence of considering only 3rd-order force constants (i.e., 3-phonon scattering processes) in the determination of phonon lifetimes when using ShengBTE. To verify this hypothesis, the impact of 4th order force constants (i.e., inlcuding 4-phonon scattering processes) was tested at 300 K and at 800 K, using an extension module for ShengBTE that incorporates four-phonon scattering [22] (see Methods section). Indeed, including higher-order anharmonicities distinctly reduces the calculated thermal conductivity, getting it closer to the experimental values (see green crosses at 300K and at 800K in the central panel of Figure 9). At 300 K, the thermal conductivity agrees remarkably well with the largest reported value for GaN by Shibata et al. [66] with a relative error of only 1.4%. At 800 K, the deviation between theory and experiment remains sizable. This is not surprising, considering that at very high temperatures even higher-order scattering processes may play a role. Additionally, factors such as impurity scattering, which strongly affects high-frequency phonons ($\propto {\omega }^4$) and thus has a more significant contribution at 800 K, are not considered here. Shibata et al. [66] report Si impurities of $2.1·{10}^{17}{\mathrm{cm}}^{-3}$, which could explain the discrepancy between their measurements (Figure 9 orange diamonds) and our results (Figure 9 lime crosses) at 800 K, even when fourth-order contributions are included.

For ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N, the situation is more involved. First, there are only a few measurements on ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N in general and even fewer that consider the target alloy fraction $x=0.5$. Moreover, to the best of the author’s knowledge, no bulk measurements are available in the literature, which is particularly unfortunate, since alloys are known to display strong size effects on the thermal conductivity. Therefore, no reliable data on the thermal conductivity of ideal bulk ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N are known to date. The closest match for our case is the data by Daly et al. [61] for a single-crystalline thin-film ${\mathrm{Al}}_{0.44}{\mathrm{Ga}}_{0.56}$N sample. From a theoretical point of view, the standard ALD-based approach has serious limitations when applied to alloys. First, the virtual crystal approximation (VCA) is used to enforce periodicity on the material, which is required for ALD. This is already a fairly big approximation. Second, the missing randomness is approximated using a lowest-order perturbative treatment of mass disorder via the Tamura model [24], while disorder in the force constants (i.e., variations in interatomic bonding strength) is neglected. A perturbative approach cannot be rigorously justified in the case of large alloy fractions. Despite these glaring problems, the method is widely used and accurately describes the thermal conductivity of many alloys [23,69,70,71]. With these caveats aside, our ${\mathrm{Al}}_{0.5}{\mathrm{Ga}}_{0.5}$N results show a reasonable agreement with the available experimental data but overeshoots the experiment, likely due to size effects. It also has to be mentioned, that the axis orientation of the sample by Daly et al. is not known and that ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N has an extremly large anisotropy ($\sim 50\%$) for the thermal conductivity as can be seen in Table 4). The scalar thermal conductivity (average over trace elements) agrees much better with the data reported by Daly et al.

To analyze the contributions of various phonon regions to the thermal conductivity, Figure 10 illustrates the average frequency-resolved contributions to the thermal conductivity ${\langle \kappa \rangle }_f$, as determined by eq:kappaDOS. As usual, the results for the material-specific MTPs, for the alloy MTP and for DFT are compared. The positions, shapes and widths of main peaks of the DFT curve are excellently reproduced by the MTPs. The only minor deviations are shifts in peak positions, which can be traced back to the minor differences in phonon bands. This clearly illustrates that the MTPs presented in this work can accurately capture even the microscopic details of thermal transport.

4. Conclusion

This work demonstrates that machine-learned force fields can be a highly effective tool for simulating the static and dynamic structural properties of ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N and related materials, including thermal transport. They combine the computational efficiency of empirical force fields with the reliability and accuracy of first-principles calculations. Specifically, moment tensor potentials (MTPs) were trained for three materials commonly used in microelectronics, AlN, GaN, and ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N, using datasets of supercells with varying atomic displacements and cell deformations. The performance of the MTPs was assessed by benchmarking thermal transport–related properties against experimental and DFT reference data.

We investigated the impact of dataset parameters, such as the number of configurations and the magnitude of atomic displacements and cell deformations. MTP prediction quality generally improves with increasing MTP levels and with larger dataset sizes. Increases beyond a certain point (∼ 300 in our case) yield little further improvement, at least at low MTP levels. Larger atomic displacements and cell deformations were found to negatively affect certain properties, such as the lattice and elastic constants. For the materials studied here, level-10 MTPs are found to be sufficient, but careful selection from an ensemble is necessary due to the large variability at lower MTP levels. When chosen this way, optimal MTPs can achieve accuracy comparable to much higher-level models while requiring only a fraction of the computational cost.

The MTPs show excellent overall agreement with DFT, with relative errors below 5% for most benchmarked properties. In some cases, such as lattice constants and phonon band structures, the MTPs reproduce the DFT results nearly exactly. Notably, the alloy MTP performs only slightly worse than MTPs trained specifically for AlN and GaN, highlighting its potential for studying ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$N alloys and heterostructures, including interfaces and superlattices. Their combination of high computational efficiency and near-DFT accuracy makes them a powerful tool for MD-based simulations.

In general, our MTP- and DFT-based calculations show good agreement with experiment. When minor discrepancies arise, they are systematic and occur in the same direction for both MTP and DFT results, suggesting that their origin lies in the underlying DFT description. We attribute these deviations to the choice of exchange–correlation functional (PBEsol) employed in our simulations, which is known to exhibit slight underbinding.