A deep neural network interatomic potential for studying thermal conductivity of β-Ga₂O₃

As a wide-bandgap transparent oxide semiconductor, β-Ga₂O₃ has attracted much scientific and technological attention.¹ Due to its large critical electric field and relatively inexpensive growth method,² β-Ga₂O₃ is a promising candidate for applications in power electronic devices,³ solar-blind UV photodetectors,⁴ and gas sensors.⁵ In particular, for high-power and high-voltage switching devices (e.g., transistors), the performance and reliability can be significantly affected by the high temperature induced by localized Joule heating in the channel of oxide semiconductors.^6,7 An understanding of the thermal properties of β-Ga₂O₃ is, therefore, critical for better thermal management and device packaging design.

A few experimental studies have reported the thermal conductivity values of β-Ga₂O₃.^8–11 It is shown that the measured thermal conductivity (κ) of β-Ga₂O₃ is anisotropic due to the monoclinic lattice structure.¹² Recent computational studies using the first-principles Boltzmann transport equation (BTE) method,^13,14 however, failed to predict κ values in good agreement with the experimental values, which is likely due to factors like the difficulties in obtaining converged κ with respect to the grid size and the broadening coefficient for screening the scattering phase space. In addition, given the large conventional cell of β-Ga₂O₃ containing 20 atoms, the first-principles calculations are highly time-consuming especially when extensive convergence tasks need to be performed. Ab initio molecular dynamics (AIMD) can also accurately predict thermal conductivity,¹⁵ but simulations of β-Ga₂O₃ will be severely limited by the system size (a few hundred atoms) and timescales (tens of picoseconds). Molecular dynamics (MD) simulation using empirical potentials enables us to efficiently predict properties of both bulk materials and nanostructures, over a wider range of system sizes (up to hundreds of nanometers)¹⁶ and timescales (up to tens of nanoseconds),¹⁷ even though the prediction accuracy highly depends on the fidelity of the potentials. However, such accurate potentials do not exist for β-Ga₂O₃, and developing them is usually a challenging task, which commonly involves fitting pre-selected functional forms with the ab initio energy surface.¹⁸ Depending on the specific problem of interest, the fixed functional form of empirical potentials can limit the applicability of the fitted potential, especially when higher order terms (e.g., force and anharmonicity) are important. Studying thermal transport is one such problem.

Thanks to the recent advances in machine learning (ML)-based techniques,¹⁹ one can reconstruct the ab initio potential energy surface using certain ML models, such as artificial neural networks²⁰ and Gaussian process regression²¹ with very high accuracy. Such ML models can be used as reliable fitting tools to develop interatomic potentials that precisely map the atomic configurations to the corresponding ab initio energies and forces. Several different ML potentials have been proposed to predict thermal properties of crystalline solids such as Si,²² Zr,²³ graphene,²⁴ and MoS₂.²⁵ Recent studies proved that even materials with more complex structures can be well simulated using ML potentials,^26,27 and a unified potential has been developed to model different phases of Si and, at the same time, predict thermal conductivity values accurately.²⁸ It has also been proved that long-range interactions are important for the correct description of thermal vibrations in covalent materials, but two-body empirical potentials using radial descriptors cannot properly capture the forces and interatomic force constants due to the lack of flexibility.²⁹ Since ML potentials can provide such flexibility and no accurate potential exists for β-Ga₂O₃, a ML potential for β-Ga₂O₃ will be highly impactful to the community.

In this work, we develop a neural network potential (NNP) to predict the thermal conductivity of β-Ga₂O₃ with MD simulations. The potential is obtained by training a deep neural network (DNN) model with the ab initio data of a β-Ga₂O₃ supercell. Accurate κ values can be achieved using the NNP, which agree well with the experimental results at different temperatures. Normal mode analysis has been performed to quantify the contributions to thermal transport from different phonon modes.

In order to build the NNP for β-Ga₂O₃, we first need to generate a reference dataset for the training process. AIMD simulations³⁰ are conducted to collect the needed data (potential energy, atomic forces, atomic coordinates, and supercell lattice vectors) of β-Ga₂O₃. We perform AIMD simulations at different finite temperatures, along with static self-consistent density functional theory (DFT) calculations of randomly displaced structures at 0 K to diversify the sampled phase space. As for the AIMD simulations at finite temperatures, snapshots are collected from different molecular trajectories in the isobaric–isothermal ensemble with temperature ranging from 50 to 600 K at zero pressure. The system is first equilibrated for 1.0 ps with a time step of 1.0 fs, and then the snapshots are collected every four steps during production runs for 5.0 ps at different temperatures. The ab initio calculations are carried out using the QUICKSTEP algorithm implemented in the CP2K package, which is based on the Gaussian and plane wave (GPW) approaches.^31,32 We use the triple-ζ valence plus polarization (TZVP) basis sets and the Goedecker–Teter–Hutter (GTH) pseudopotentials³³ to describe core–valence interactions. The plane wave cutoff is set to be 800 Ry, and the Brillouin zone is sampled at the Γ point. All the calculations are performed for a relaxed monoclinic 2 × 4 × 2 β-Ga₂O₃ supercell, with the lattice parameters a = 12.376 Å, b = 3.084 Å, c = 5.893 Å, and β = 103.83° [Fig. 1(a)], in good agreement with the experimental results.³⁴ A total of 9200 ab initio snapshots are generated for the training process.

An end-to-end DNN-based potential scheme is used to obtain the NNP for β-Ga₂O₃. This MD potential scheme, namely, the Smooth Edition of Deep Potential (DeePot-SE) developed by Zhang et al.,³⁵ represents the lattice potential energy as a sum of atomic energies, which can be determined by the local environment of each atom. Within DeePot-SE, all the natural symmetries are preserved by mapping the local environments to an embedded feature space. More details about this scheme and related models can be found in the literature.^36,37 In this study, our DNN model is composed of an embedding network and a fitting network, and the only inputs are the chemical species and atomic coordinates. We have used a two-layer embedding network to map the original inputs to symmetry-preserving components. The fitting network is a three-layer feedforward neural network containing 160 nodes in each layer, which maps the embedded features to the atomic energies. The DNN model is trained with the ab initio dataset to minimize the loss function, which is a weighted sum of mean square errors in energy and force. The optimization is solved by TensorFlow's implementation of the Adam stochastic gradient descent method.³⁸ The DNN structures and training hyper-parameters are also optimized to minimize the prediction errors. All training data and potential files are included in Ref. 39. Then, we use the LAMMPS package⁴⁰ with this established NNP to perform MD simulations of β-Ga₂O₃.

The fitting accuracy can be verified by the agreement between NNP predictions and the benchmark AIMD results. We compare the NNP-predicted energies and forces for a separate testing dataset containing 600 AIMD snapshots that are randomly selected and not included in the training dataset. It is found that our NNP can provide very similar results compared to the testing data at different temperatures, with the root mean square errors of 0.0005 eV per atom for the potential energy and 0.078 eV/Å for the atomic forces. The mean percent errors in energy and force are 0.0003% and 9.1%, respectively. For comparison, it has been reported that popular empirical potentials such as Tersoff⁴¹ and Stillinger–Weber⁴² potential for crystalline Si result in 35% and 210% error in forces, respectively.⁴³ Figure 1(b) shows the comparison of AIMD energies and forces and NNP-predicted quantities, indicating that the NNP accurately reproduces the ab initio energies and forces.

To evaluate the accuracy of harmonic force constants, we calculate the phonon dispersion relation along selected high symmetry paths. The harmonic force constants are computed by fitting displacement–force data using the Alamode program,⁴⁴ where symmetry-irreducible displacements of 0.01 Å are imposed and atomic forces are calculated with this NNP. As shown in Fig. 2(a), the phonon dispersion using NNP agrees well with that from the first-principles calculation by Santia et al.,¹³ while existing empirical potentials fail to capture that.^45,46 Considering the low symmetry and high anisotropy of β-Ga₂O₃, it is reasonable to see that the phonon dispersion profile is much more complex than that of many other wide-bandgap materials, such as GaN and ZnO.⁴⁷ At the same time, it is understandable why it has been very difficult to achieve an excellent agreement in phonon dispersion using a fitted empirical potential composed of simple functions. However, because of the versatility arising from high complexity of the DNN, harmonic force constants can be reproduced accurately with NNP, which are important for the prediction of thermal conductivity.

Equilibrium MD (EMD) simulations are preformed to evaluate the thermal conductivity using the Green–Kubo approach.^48–50 The Green–Kubo formula relates the instantaneous fluctuations of heat current, in terms of autocorrelation function, to thermal conductivity as

Here, κ_α denotes the thermal conductivity in the α direction, k_B is the Boltzmann constant, V is the volume of the simulation cell, T is the temperature, τ₀ is the integral upper limit, which should be infinite theoretically, t represents the correlation time, $$ denotes the ensemble average, and J_α is the component of heat current J in the α direction, which is expressed as⁵¹ 

where E_i, v_i, U_j, and r_ij are the site energy, atomic velocity, potential energy, and distance between atoms i and j, respectively. Conventionally, one defines $∑iviEi$ as the kinetic part, which is found to barely contribute to the thermal conductivity for solids and can be simply neglected.⁵² We, thus, only focus on the potential part $∑i∑j≠irij∂Uj∂rij·vi$ for the heat current in this work.

We perform EMD simulations on a bulk β-Ga₂O₃ system consisting of 4 × 13 × 7 conventional unit cells consisting of 7280 atoms (at least 4 nm long in each crystallographic direction), which is usually not feasible to simulate using first-principles methods. Periodic boundary conditions are applied in all three directions, and the isobaric-isothermal ensemble is first used to equilibrate the system at zero pressure for 200 ps with a time step of 1.0 fs. Then, the simulations are switched to the microcanonical ensemble and run for another 2 ns to collect the heat current data. To suppress the statistical uncertainty, 15 independent MD runs with different initial velocity distributions are performed to obtain the average κ values at 200, 300, 400, and 500 K. As mentioned in Ref. 28, MD simulations using NNP are usually slower than those using empirical potentials. With the potential in this work, it takes around 6 days to simulate a 7280-atom system for 2 ns using a single Intel Xeon Silver 4110 CPU core and a NVIDIA TITAN Xp GPU, which suggests that our calculations of thermal conductivity are in a reasonable computing time range.

Figure 3(a) shows the average integral of heat current autocorrelation functions along three crystallographic directions at 300 K. The inset shows the thermal conductivity values for various simulation cell sizes, with the 7280-atom system giving results within the error bar of the case for a 5040-atom system, which are considered converged with respect to the domain size. We have observed the anisotropy in thermal conductivity, with the calculated values at 300 K to be κ_[100] = 10.68, κ_[010] = 20.78, and κ_[001] = 12.61 W m⁻¹ K⁻¹. The temperature dependence of κ is plotted in Fig. 3(b) along with the results from time-domain thermoreflectance (TDTR) measurements.^8,11 Predictions from first-principles BTE calculations are also included for comparison.^13,14 The calculated κ values using NNP show good agreement with experimentally observed trends, with the highest κ along the [010] direction. However, for the [100] and [001] directions, both NNP predicted results and experimental results are lower than the first-principles results. One of the reasons could be that MD simulations with NNP naturally include higher order anharmonicities, whereas the BTE method truncates the anharmonic force constants to the third order.⁵³ Moreover, another observation is that the BTE-calculated κ values in the [010] and [001] directions are essentially the same, which contradicts with both experimental measurements^8,11 and our MD calculations. We, thus, cast our doubts on the BTE results,¹³ which are challenging to obtain accurately given the complexity of the β-Ga₂O₃ unit cell and uncertainties in convergence with respect to different parameters in such calculations as previously mentioned. At high temperatures (>400 K), our calculated κ values in all directions agree quantitatively with experimental measurements. At lower temperatures (e.g., 200 K), the MD results are lower than experiments, which is likely due to the size effect of the simulation domain and the overestimation of the scattering rates that cutoff the contributions from long wavelength phonons, which are more important at lower temperatures due to their increased mean free path. Here, we have not accounted for any quantum effect in our calculations, and there can be two competing quantum effects: classical statistics gives larger modal heat capacities at high frequencies but smaller phonon scattering times at low frequencies compared to quantum statistics.⁵⁴ Since the quantum correction is still controversial in accounting for quantum effects, we do not include quantum correction in this work. However, for practical applications, β-Ga₂O₃ thermal conductivity above room temperature is more important, where our NNP-EMD is able to yield good results.

To further study the phonon transport physics, we calculate the contribution of different phonon modes to the anisotropic thermal conductivity using Green–Kubo modal analysis (GKMA).^55,56 The GKMA method combines the lattice dynamics with the Green–Kubo formula such that the thermal conductivity becomes a direct summation of modal contributions. For a system with N atoms vibrating around their equilibrium positions, one can obtain 3N collective modes of vibration by lattice dynamics calculations on the entire supercell. Then, the MD trajectory can be projected onto the modes by relating the velocity of each atom j to the normal mode velocity $Ẋn$ for mode n,

where m_j is the mass of atom j and e_j,n is the eigenvector that gives the magnitude and direction of motion for atom j in mode n. The normal mode velocity can be obtained by reverse transformation,

where * represents the complex conjugate. Neglecting the kinetic part, we can write the heat current in Eq. (2) as a sum of individual modal contributions,

By substituting Eq. (5) into Eq. (1), we have the contribution of mode n in thermal conductivity

Using the modes obtained from the finite difference method, we perform a modal analysis based on the GKMA method for β-Ga₂O₃ at 300 K. The accumulation of thermal conductivity with respect to mode frequency is shown in Fig. 4. Our GKMA results using NNP show that around 80% of κ in three crystal directions are contributed by the phonon modes below 10 THz. It is also evident that the optical phonons with frequency above 5 THz play a critical role in the thermal transport, especially in the [010] direction. At 300 K, optical phonons with frequency above 5 THz contribute over 45% of the total κ_[010]. The high optical phonon contribution may result from the hybridization of low energy optical modes with the longitudinal acoustic mode.^13,14 Furthermore, compared to [100] and [001] directions with relatively flat optical phonon branches [Fig. 2(b)], the [010] direction shows larger phonon group velocities, which can lead to a higher κ contributed by optical phonons.

In summary, we have developed an interatomic potential for β-Ga₂O₃ using a deep neural network model. This NNP is trained with ab initio energies and forces such that it can reproduce the phonon dispersion in very good agreement with the first-principles results. Using MD simulations with this NNP, we predict the lattice thermal conductivity in three crystallographic directions over a range of temperatures from 200 to 500 K. The thermal conductivity of β-Ga₂O₃ is found to exhibit large anisotropy, and our predictions agree well with experimental results in the literature especially above room temperature. Furthermore, the modal analysis using the GKMA method reveals that optical phonon modes contribute significantly to the overall thermal conductivity at 300 K, especially in the [010] direction. These results demonstrate the accuracy of our NNP in predicting thermal properties of β-Ga₂O₃. Since the accurate prediction of thermal conductivity implies the accuracy of our NNP in both second and third order force constants, this potential can be a foundation for MD simulations of other β-Ga₂O₃ properties, such as mechanical and thermo-mechanical properties. In addition to β-Ga₂O₃, more NNPs can be developed in the same vein to study thermal transport in other complex semiconductor materials.

The authors would like to thank ONR MURI (No. N00014-18-1-2429) for the financial support. The simulations are supported by the Notre Dame Center for Research Computing and NSF through the eXtreme Science and Engineering Discovery Environment (XSEDE) computing resources provided by Texas Advanced Computing Center (TACC) Stampede II under Grant No. TG-CTS100078.