Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Neural‐Network Potential Simulation of Defect Formation Induced by Knock‐On Irradiation Damage in GaN

Neural‐Network Potential Simulation of Defect Formation Induced by Knock‐On Irradiation Damage in... IntroductionGallium nitride is widely used in optoelectronic and electronic devices, such as high electron mobility transistors , light emitting diodes (LEDs), laser diodes, and photodetectors, mainly because it has a wide direct bandgap (3.4 eV) and high electron mobility (2000 cm2 V−1 s−1).[1–4] It also exhibits excellent resistance against irradiations, making it widely used in aerospace and military fields for irradiation‐resistant devices. Its irradiation resistance comes from the strong ionic‐covalent bonding, great thermal stability, and high displacement energy.[5] For both radiation‐resistant and electronic/optoelectronic applications, irradiation damage not only happens in the radiative environment, but also commonly exits in the ion implantation process during device fabrication.[6] The defects caused by irradiation damage greatly impact the device performance, degradation, and reliability. Therefore, understanding the microscopic mechanism of the irradiation damage of GaN is of great importance.Due to the complexity and the picosecond‐scale of the irradiation damage process, it is rather difficult to directly locate the defects and monitor the evolution of the defects using experimental techniques.[7] As a result, the threshold displacement energy (denoted as Ed), which means the least kinetic energy required for an atom to make displacement happen, can only be measured indirectly. For example, by observing the luminescence of GaN LEDs under irradiation, Ionascut‐Nedelcescu et al.[8] determined the Ed of the Ga atoms was determined to be 19±2 eV based on the assumption that the change of luminescence is mainly caused by the displacement of Ga atoms. On the contrary, the investigation of the microscopic mechanism of the irradiation damage mostly relies on the theoretical simulations. For example, Look et al.[9] analyzed the transport characteristics of GaN thin films irradiated by electron beam and obtained the Ed of GaN through fitting the relativistic cross section of atomic displacement as the function of electron energy (Ed (Ga) = 20.5 eV, Ed (N) = 10.8 eV). A more common theoretical approach to utilize molecular dynamics (MD), including the classical MD and the ab initio MD (AIMD), to simulate the atomic displacement process and determine the threshold displacement energy. Nord et al.[10] proposed an empirical Ziegler–Biersack–Litmark (ZBL) shielding potential of GaN and study the irradiation damage, however, the minimal and the average Ed of Ga atoms are determined to be lower than N atoms. He et al.[11] improved the potential by including the short‐range interaction, and the average Ed of Ga is determined to be 45 eV. Interestingly, using ab initio MD, Xiao et al.[12] also evaluated the average Ed of Ga and N as 73.2 and 32.4 eV, indicating the discrepancies between classical MD and AIMD.Conventionally, AIMD simulations are regarded to be more accurate and reliable than classical MD, because the empirical potential cannot accurately describe the saddle states. However, first‐principles calculations are slower and typically limited to hundreds of atoms, preventing long or large‐scale simulations as in classical MD simulations. Therefore, it is of great significance if one can construct a classical MD potential that is able to provides similar accuracy as AIMD. In recent years, significant progress has been achieved in the construction of the potentials with the help from artificial intelligence and machine learning, such as the potentials for Al, MoS2, SiO2, ZnO, etc. Wang et al.[13] constructed a potential for Al based on a hybrid scheme that interpolates the ZBL screened nuclear repulsion potential with a deep learning potential (DP) energy model. Using this potential, they not only obtained the accurate threshold displacement energy of Al, but also simulated the high‐energy radiation damage processes using a supercell containing 600 000 atoms. They also used this scheme to successfully develop the potential of MoS2.[14] However, the DP model is constructed based on the structures from the MD trajectory at different temperatures, making it difficult to overcome the high reaction barriers at low temperatures and causing the system gets trapped in high‐entropy structures at high temperatures. To solve this problem, Liu et al. developed the stochastic surface walking (SSW) algorithm for the structure prediction and pathway searching.[15] It is an unbiased global optimization method which can simultaneously study the minimum and saddle point on potential energy surface (PES). Based on SSW, they proposed a potential training scheme for potential called stochastic surface walking global exploration with global neural network potential (SSW‐NN), which has been successfully applied to SiO2[16] and ZnO.[17] The potentials constructed using this approach are also called global neural network (G‐NN) potentials. Compared with other optimization methods, SSW algorithm can cover more energy surface and explore unbiasedly both minima and saddle points on PES.For GaN, Shimizu et al.[18] proposed a NN potential to calculate the phonon modes and thermodynamics properties of large‐scale defective structures, in which the atoms are assumed to move around the equilibrium positions. However, it is not designed for study the irradiation damage processes, and the atoms are potentially knocked much further away from their equilibrium position. Therefore, in this work, we use the SSW‐NN method to construct the Ga‐N binary G‐NN potential that can be applied to simulate the radiation damage process. By comparing the MD simulations using our G‐NN potential with ab initio calculations, our potential is benchmarked to produce similar accuracy in multiple dimensions. The energies and formation energies of different defect configurations are in good agreement with the ab initio results, not only for the regular intrinsic point defects but also for the defect‐pairs. Such potential is also applied to simulate the knock‐on process of the irradiation damage, and the trajectory and the displacement threshold energy also in good agreement with the ab initio calculations. The global neural network (G‐NN) potential constructed in this work can elevate the accuracy and the efficiency of the MD simulations of the irradiation damage of GaN.Results and DiscussionEvaluation of the G‐NN PotentialThe data set initially contains the elemental phase of Ga and N2, as well as the GaN configurations in which Ga:N ratio is kept as 1:1 (GanNn). After 100 000 steps of training based on this data set, a potential with energy root mean square error (RMS) of 1.056 meV per atom and force RMS of 0.036 eV Å−1 is created. We selected 10 GaN structures, 5 Ga elemental structures, and 5 N2 molecular structures to benchmark this potential, resulting an average energy difference of 0.5 meV per atom from density functional theory (DFT) results. To improve the accuracy of our potential under different chemical environment, we increase the data set configurations with different stoichiometry, varying the Ga:N ratio from 0.06 to 0.94. The configurations are also generated using SSW‐NN approach, and the final training data set consists of 61 414 structures (briefly summarized in Table S1 of the Supporting Information). The RMS errors of the energy and the force of our final G‐NN potential are 5.373 meV per atom and 0.090 eV Å−1, respectively, indicating that our potential can accurately describe the defect properties as well as the irradiation process.Benchmark of Defect Properties Using the G‐NN PotentialIn the binary compound GaN, there exist 6 types of intrinsic point defects, including two vacancies (VGa, VN), two interstitials (Gai, Ni) and two antisites (GaN, NGa). To evaluate the accuracy of our potential on defect properties of GaN, we first calculate the defect formation energies for these 6 intrinsic point defects using our G‐NN potential, and compare them with the DFT calculations, as shown in Figure 1. To illustrate the improvement of our potential by adding the nonstoichiometric structures into the training data set, we also plot the formation energies calculated using the potential purely trained from the initial data set that only contains the stoichiometric configurations. Although the potential trained from the stoichiometric configurations can accurately describe the bulk structure and energy of GaN crystal, it shows significant error on the defect formation energies as compared to the DFT results, because the point defects have nonstoichiometric chemical environment. On the contrary, our final G‐NN potential provides formation energies with errors lower than 0.1 eV for five point defects, except that of Ga vacancy is calculated to be 0.28 eV higher than the DFT result. The agreement between our results using G‐NN potential and DFT result indicates that our potential accurately describe not only the bulk structure of GaN but also the point defects. In addition, we also use the potential to optimized the atomic structures of the point defects, and the atomic positions also matches the DFT calculations.1FigureCalculated single defect formation energies based on the NN potential and DFT calculations (under Ga‐rich condition).In order to demonstrate the robustness of the potential, our G‐NN potential is further used to predict the formation energies of different defect‐pairs in GaN. Usually, the formation energies of these defect‐pairs are high, and their concentrations in GaN samples prepared by conventional methods are relatively low. However, in the irradiation damage, these defect‐pairs may exist in large amount. In fact, a series of defects are generated along the collision cascade trajectory during the irradiation damage process, leading to more complicated chemical environment than the simple intrinsic point defects under thermal equilibrium. Therefore, to prove the robustness of our G‐NN potential, it is necessary to evaluate its performance on more complicated situations, such as the defect‐pairs. We utilize the G‐NN potential to predict the formation energies of 20 different defect‐pairs with the optimized structures, and compare them with the DFT results (Figure 2a). The formation energies predicted by G‐NN potential are generally in good agreement with the DFT results with errors lower than 0.5 eV.2Figurea) Comparison of the formation energies of 20 defect‐pairs calculated using the G‐NN and the DFT (under Ga‐rich condition). Their structures are all optimized. b) Comparison of the formation energies of 118 defect‐pairs calculated with the G‐NN and the DFT. c) The error of 118 defect‐pairs formation energies calculated from the difference between NN and DFT (in Ga‐rich GaN). We use 6 colors and shapes to mark different configurations of one type of defect‐pair and their errors increase in the order from 1 to 6.To further evaluate the accuracy, we also change the relative distance between the defects in a pair by 5 times for each pair, resulting a total 100 more defect‐pair configurations. A summary of the difference between G‐NN potential results and DFT results for the total 120 defect‐pair configurations is shown in Figure 2b,c. Most defect‐pairs have errors below 0.5 eV, and only three configurations have errors above 1 eV, indicating the robustness of our G‐NN potential on handling the defect‐pairs. The defect‐pairs with the largest error are VGa−VGa and VGa−NGa, which can be attributed to two possible reasons. On one hand, the G‐NN potential prediction of VGa already has relatively more significant error than the others, so it is reasonable that the prediction of VGa−VGa and VGa−NGa pairs suffer bigger error. On the other hand, we observe that the error generally grows with the values of formation energies, and the formation energies of these two defect‐pairs are both higher than the others. This behavior is because most of the areas explored by SSW are still on the lower potential energy surface, leading to more significant error when describing the structures with higher energy. The energy and formation energies of these structures with poor prediction are listed in Table S2 (Supporting Information). Except these structures, our potential can accurately describe the properties of various point defects and defect‐pairs in GaN, enabling us to use the potential to optimize different defect configurations and simulate the defect evolution process in a larger system.In addition, we also examine the accuracy of our G‐NN potential for a larger supercell size. Based on a 360 atom supercell, we calculate the energy of the supercell and the formation energy of VGa, and the results are in good agreement with the DFT calculations, meaning that no extra error is introduced with the growth of the model supercell. More details of the supercell structure energies and defect formation energies are listed in Table S3 (Supporting Information). Due to higher computational cost of DFT calculations, we do not compare the defect formation energies obtained using G‐NN potential with DFT results for larger supercells, however, the presented results prove that our G‐NN potential can be used to treat the defect properties of much larger systems almost as accurately as DFT.It is worth discussing that we do not consider the charged defect states in our potential, which is different from the NN potential proposed by Shimizu et al.[18] As mentioned, the potential from Shimizu et al. is designed for calculating the phonon modes and thermodynamics properties, where the atoms are assumed to move around the equilibrium positions under thermodynamic equilibrium, and the charged states are relatively important. On the contrary, the G‐NN potential in this study is used for simulating the irradiation damage processes, in which the primary knock‐on atom (PKA) is given high kinetic energy causing cascade collisions. As a result, the atoms are knocked much further away from their equilibrium position. In such cases, the difference in the collision simulations caused by the different charged state should not be significant. One evidence is that the atomic structures of neutral and charged defects are very similar, compared to the atomic displacements in cascade collisions. However, it does not mean the charged defects are not important, and it is definitely worth studying if the LASP package can support the potential fitting for charged systems in the future.Simulation of Irradiation Damage Process Using the G‐NN PotentialOur G‐NN potential is subsequently used to simulate the irradiation damage process. We first calculate the displacement threshold energy Ed of GaN for several main crystallographic directions using G‐NN potential, and they are compared to the Ed obtained using AIMD as well as the results from Xiao et al. using DFT (Table S4, Supporting Information). As shown in Table S4 of the Supporting Information, the MD simulations using G‐NN potential can provide predictions of Ed as accurate as the ab initio MD simulations. The DFT‐calculated Ed is about 39 eV for Ga and 17 eV for N along b direction, which is the minimal Ed among all crystallographic directions, and the G‐NN potential predicts the Ed as about 40 eV for Ga and 17 eV for N, agreeing with the DFT calculations. As comparison, the classical potentials fail to provide accurate simulations of the irradiation damage processes. Specifically, instead of collision, the PKA atoms directly pass through other atoms for initial kinetic energy higher than 120 eV, if the Tersoff potential and Tersoff‐ZBL potential are used.In order to validate that the generated G‐NN potential can describe the irradiation damage processes with AIMD accuracy, we further simulate the knock‐on displacement processes of Ga/N primary knock‐on atoms. We give the initial kinetic energy of Ed in the a and b directions to the Ga/N PKAs, respectively. Figure 3 shows the potential energy evolutions of the knock‐on processes. For all processes, the potential energies along the predicted using our G‐NN potential agree well with the AIMD trajectories of the potential energies, indicating that our potential is capable of simulating the knock‐on process with AIMD‐level accuracy. In addition, using the atomic configurations along the trajectories of the AIMD simulations, we also use our G‐NN potential to calculate the potential energies, and our G‐NN potential almost identically reproduces the potential energies of the AIMD results.3FigureThe evolutions of the potential energy during the displacements of a) 19 eV and b) 17 eV N PKA, as well as c) 73 eV and d) 40 eV Ga PKA.Furthermore, we also examine the atomic positions during the simulations of cascading collisions using G‐NN potential and AIMD. As shown in Figure 4, we plot the atomic configurations around the PKA atom, which is a Ga PKA with initial 40 eV kinetic energy along the [1¯010]$\bar{1}010]$ direction (b direction), at 50, 150, and 250 fs. A VGa and an octahedral Gai are formed in both 150 and 250 fs simulations, and the atoms are in relatively similar positions, indicating that the G‐NN potential can be used to simulate the cascade collisions during the irradiation damage process with ab initio accuracy.4FigureThe trajectories of a Ga 40 eV PKA along [1¯010$\bar{1}010$] direction simulated by a–d) DFT and e–h) NN. The red and blue balls represent Ga PKA and Ga interstitial, respectively, and the green arrow represents the direction of Ga PKA. The two axes a and b represent [1¯21¯0]$\bar{1}2\bar{1}0]$ and [1¯010$\bar{1}010$] direction, respectively.Last but not the least, we apply the G‐NN potential to carry out a relatively large‐scale simulations that are conventionally beyond the capability of AIMD. A 2 3040 atom supercell is used to simulate the knock‐on displacement processes of 500 eV Ga PKA along 11 different directions, and each calculation is carried out for 1000 fs. By monitoring the defect types and their populations, we observe that the number of induced defects first increases before it reaches a peak value, and decreases after due to the recovery of defects. The Wigner–Seitz cell method[19] is employed to identify vacancies, interstitials and antisite defects. A cartoon of the defect evolution is shown in Figure S1 of the Supporting Information, and only displaced atoms or defects are plotted. The peak damage is found to happen at 0.24 ps, and the number of defects keep decreasing afterward. Based on the statistics of the defect types and their populations, the relative ratios of all six intrinsic defects are determined by averaging the numbers of defects for all directions. The defect ratios for the peak damage and the final state are both listed in Table S6 of the Supporting Information. For all directions, more Ga and N vacancies, as well as N interstitials are generated than the Ga interstitials and antisites. This finding contrasts with the Stopping and Range of Ions in Matter (SRIM) and technology computer‐aided design (TCAD) simulations of GaN devices in which only Ga vacancies are considered.[20–22] To confirm our results, a larger supercell of 184 320 atoms is also used to repeat the knock‐on displacement process simulations, and the results agree with those from the 23 040 atom supercell. It is worth noting that our G‐NN potential not only enables the MD simulations of large‐scale systems with accuracy that is comparable to ab initio simulations, but also provide parameters for device simulations. The displacement threshold energy and the information of the relative defect populations determined in MD simulations can be utilized as the input parameter of the TCAD simulations of the proton irradiation damage of GaN‐based devices, such as AlGaN/GaN HEMTs.ConclusionsIn summary, we utilize machine learning algorithms and global PES exploration to develop a G‐NN potential for GaN. Based on a dataset of 61 414 structures, the potential is trained to reduce the RMS errors of energy and force as 5.373 meV per atom and 0.090 eV Å−1. The formation energies of the point defects and defect‐pairs calculated using G‐NN potential are in good agreement with the results of DFT, indicating that the potential is capable of accurately describing the defect properties. In addition, using the G‐NN potential, we also investigate the threshold displacement energies of PKA atoms along different directions, the trajectories and the potential energies of several knock‐on processes, all of which are in good agreement with AIMD calculations. Finally, based on the simulations of the knock‐on processes of a Ga PKA atom, the relative ratios of the induced intrinsic point defects are obtained. Ga and N vacancies, as well N interstitials are found to have higher amounts than the Ga interstitials and the antisites. The G‐NN potential extends the scales of time and size of the atomic‐level simulations of the irradiation damage processes while keeping similar accuracy as ab initio approach. Therefore, it is not only able to provide the fundamental understanding of the irradiation damage processes, but also to refine the parameters such as the threshold displacement energy in the device simulations.Experimental SectionSSW‐NN MethodThe constructions of the G‐NN potential were conducted using LASP code[23] (Large‐scale Atomic Simulation with neural network Potential) developed by Liu et al. group. The random stochastic surface walking (SSW) method[15] was utilized for the global optimization, which can automatically explore the multidimensional PES and identify important structures and reaction paths on the potential energy surface. The G‐NN potential follows the atom‐centered NN potential architecture proposed by Behler and Parrinello,[24] where the input layer utilized the power‐type structure descriptor (PTSD).[25] The G‐NN potential is sampled by SSW global potential energy surface, fitted and iteratively generated by data obtained from high‐precision DFT calculation, until the G‐NN potential has sufficient robustness to quantitatively describe global PES, or namely the SSW‐NN.[26]The training process of the potential divided into three main stages as illustrated in Figure 5. The first stage was to construct an initial dataset containing the most common atomic environment of the system as a training set for building the initial global PES. In this stage, the dataset only includes structures with Ga:N ratio as 1:1, and the elemental phase of Ga and N2. The random stochastic surface walking (SSW) was used to automatically explore and obtain the multidimensional potential energy surface. In total, 247 246 structures were generated, and 6030 of them were randomly selected into the dataset and simulated using first‐principles method. These structures collected and screened from the SSW trajectory were calculated using first‐principles method and added into the dataset. Eventually a dataset of over 6000 configurations is built.5FigureThe schematic plot of the construction process of the G‐NN potential.In the second stage, based on the dataset obtained in the first stage, LASP was used to carry out the initial potential training. In order to ensure the accuracy of the potential, a large number of PTSDs, including 56 two‐body, 60 three‐body, and 8 four‐body descriptors was used. The general forms of the PTSD are described in the following equations1fcrij=0.5×tanh31−rijrc,forrij≤rc0forrij>rc$$\begin{eqnarray}{f}_{\rm{c}}\ \left( {{r}_{{\rm{ij}}}} \right) = \left\{ { \def\eqcellsep{&}\begin{array}{@{}*{1}{c}@{}} {{\rm{0}}{\rm{.5}} \times \tan {{\rm{h}}}^3\left( {1 - \frac{{{r}_{{\rm{ij}}}}}{{{r}_{\rm{c}}}}} \right),for\ {r}_{{\rm{ij}}} \le {r}_{\rm{c}}}\\[6pt] {\ 0\ for\ {r}_{{\rm{ij}}} > {r}_{\rm{c}}} \end{array} } \right.\ \end{eqnarray}$$2Si1=∑j≠iRnrij$$\begin{eqnarray}S_{\rm{i}}^1 = \mathop \sum \limits_{j \ne i} {R}^n\ \left( {{r}_{{\rm{ij}}}} \right)\end{eqnarray}$$3Si2=∑m=−LL∑j≠iRnrijYLmrij212$$\begin{eqnarray}S_{\rm{i}}^2 = {\left[ {\mathop \sum \limits_{{\rm{m}} = - L}^L {{\left| {\mathop \sum \limits_{j \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right){Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right)} \right|}}^2} \right]}^{\frac{1}{2}}\ \end{eqnarray}$$4Si3=21−ζ∑j1+λcosθijkζ·Rnrij·Rmrik·Rprik$$\begin{eqnarray}S_{\rm{i}}^3 = {2}^{1 - \zeta }\ \mathop \sum \limits_j {\left( {1 + \lambda \cos {\theta }_{{\rm{ijk}}}} \right)}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{ik}}}} \right) \cdot {R}^p\left( {{r}_{{\rm{ik}}}} \right)\end{eqnarray}$$5Si4=21−ζ∑j,k,l≠i[1+λcosθijkζ·Rnrij·Rmrjk$$\begin{eqnarray}S_i^4 = {2}^{1 - \zeta }\ \mathop \sum \limits_{j,k,l \ne i} [{\left( {1 + \lambda \cos {\theta }_{{\rm{ijk}}}} \right)}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{jk}}}} \right)\end{eqnarray}$$6Si5=∑m=−LL∑j,k≠iRnrij·Rmrik·Rprjk·YLmrij+YLmrik212$$\begin{eqnarray}S_{\rm{i}}^5 = {\left. {{{\left\{ {\mathop \sum \limits_{m = - {\rm{L}}}^{\rm{L}} \left| {\mathop \sum \limits_{j,k \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{ik}}}} \right) \cdot {R}^p\left( {{r}_{{\rm{jk}}}} \right) \cdot \left[ {{Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right) + {Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ik}}}} \right)} \right]} \right|} \right.}}^2} \right\}}^{\frac{1}{2}}\ \end{eqnarray}$$7Si6=21−ζ∑j,k,l≠i1+λcosδijklζ·RnrijRmrikRpril$$\begin{eqnarray}S_{\rm{i}}^6 = {2}^{1 - \zeta }\ \mathop \sum \limits_{j,k,l \ne i} \left[ {{{\left( {1 + \lambda \cos {\delta }_{{\rm{ijkl}}}} \right)}}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right){R}^m\left( {{r}_{{\rm{ik}}}} \right){R}^p\left( {{r}_{{\rm{il}}}} \right)} \right]\end{eqnarray}$$8Si7=∑j≠iRnrij·1+tanhα·rij−rm$$\begin{eqnarray}S_{\rm{i}}^7 = \mathop \sum \limits_{j \ne i} \left\{ {{R}^n} \right.\ \left. {\left( {{r}_{{\rm{ij}}}} \right) \cdot \left[ {1 + \tan {\rm{h}}\alpha \cdot \left( {{r}_{{\rm{ij}}} - {r}_{\rm{m}}} \right)} \right]} \right\}\end{eqnarray}$$9Si8=∑m=−LL∑j≠iRnrijYLmrij·1+tanhα·rij−rm212$$\begin{eqnarray}S_{\rm{i}}^8 = {\left[ {\mathop \sum \limits_{{\rm{m}} = - {\rm{L}}}^{\rm{L}} {{\left| {\mathop \sum \limits_{j \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right){Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right) \cdot \left[ {1 + \tan {\rm{h}}\left[ {\alpha \cdot \left( {{r}_{{\rm{ij}}} - {r}_{\rm{m}}} \right)} \right]} \right]} \right|}}^2} \right]}^{\frac{1}{2}}\ \end{eqnarray}$$whererijis the internuclear distance between atomiandj,θijkis the angle centered at iatom withjand kbeing neighbors (i,jand kare atom indices). The key factors of PTSD are the cut‐off functionfc, which is decayed to zero beyondrc, power‐type radial function, trigonometric angular functions, and spherical harmonic function. In PTSD, theS1andS2are two‐body functions, theS3,S4, andS5are three‐body functions, and theS6is a four‐body function.The G‐NN potential has a five‐layer (128‐80‐80‐48‐1) feed‐forward NN structure, in total 40 834 fitting parameters. The min‐max scaling method[27] was used to normalize the training dataset. Hyperbolic tangent function was used to activate the hidden layer, while linear transformation was performed on the output layer of all networks. The quasi‐Newton Limited‐memory Broyden–Fletcher–Goldfarb–Shanno (L‐BFGS)[28] method was used to minimize the loss function to match the energy, force and stress from the DFT calculations. In the process of potential training, the number of training steps set to 50 000, and the training weights of energy, force and stress set to 200, 15, and 1, respectively.In the third stage, the global dataset is expanded iteratively to improve the predictive power of the neural network potential. The SSW‐NN method was still adopted as in the second stage, however, the structures found by SSW‐NN are calculated using the constructed potential and classical MD, which is different from the DFT calculations used in the first stage. Structures collected and filtered from the SSW trajectories were used to expand the dataset. The added configurations included 45 categories of nonstoichiometric structures with Ga/N ratio ranging from 1:15 to 15:1, and the numbers of atoms in a unit cell were in the range of 12–128 atoms. The purpose of the potential is focused on the point defect, which is localized in a small area of the whole crystal. As a result, small nonstoichiometric unit cells are good enough to describe the chemical environment of point defects. Adding a few large supercells is necessary to utilize all structure descriptors in the potential training process adequately. Such rule of constructing the data set not only enables more accurate defect descriptions but also accelerates the training process of potential. Then the new dataset was input to the global dataset to start a new neural network training cycle. The G‐NN potential is obtained until an additional training cycle does not obviously improve the accuracy potential.DFT CalculationsAll density functional theory calculations were performed by using VASP[29] (Vienna Ab‐initio Simulation Package) . The projector‐augmented wave method[30–32] and a 520 eV energy cutoff were used for all calculations. The generalized gradient approximation (GGA) exchange‐correlation functional in the form of Perdew–Burkes–Ernzerhof (PBE)[33] was used. The Monkhorst‐Pack k‐point grid was adopted, and the density of the grid was dependent on the size of the first Brillouin zone. The energy convergence criterion was set as 1 × 10−5 eV, and the lattice was optimized until the force on each atom is less than 0.01 eV Å−1.For the calculation of GaN defects, we only consider defects at neutral charge state. A quasicubic 128 atom supercell is selected to calculate the defect properties. The defect formation energy is calculated using[34,35]10ΔHfa,q=Ea,q−Ebulk+∑niμi+qEF$$\begin{eqnarray}\Delta {H}_{\rm{f}}\ \left( {a,q} \right) = \ E\left( {a,q} \right) - E\left( {{\rm{bulk}}} \right) + \sum {n}_{\rm{i}}{\mu }_{\rm{i}} + q{E}_{\rm{F}}\end{eqnarray}$$where E(a, q) and E(bulk) represent the total energies of supercells with and without defect, respectively. ni is the number of atoms that the supercell with the defect differs from the supercell without defect. µi represents the chemical potential of the corresponding element, and this chemical potential is changed with reference to the total energy Ei of the elemental phase of the element. EF denotes the Fermi level. For a certain chemical potential condition, the defect formation energy of electric neutral is a constant for Fermi level in the whole range of the bandgap. For the convenience of later comparison, the formation energy calculated later is obtained under the condition of Ga‐rich.MD SimulationsAll MD simulations were carried out using the LAMMPS (Large‐scale Atomic/Molecular Massively Parallel Simulator) package,[36] which has been widely used to simulate the irradiation damage process of materials including metals and semiconductors. The atomic visualizations were created using the Open Visualization Tool (OVITO).[19] All simulations were based on a 2880 atom supercell with lattice constants of 32, 33, and 31 Å, and periodic boundary conditions were imposed in all directions. The NVT ensemble was used to relax the crystal for 1000 timesteps at 100 K with time step of 1 ps. Once the system reached equilibrium, a Ga atom or N atom was selected as the PKA and given initial kinetic energy to simulate cascade collisions based on NVE ensemble. The whole process lasted 5000 timesteps, and the time step adaptively changed to avoid instability of simulation results. If a stable Frankel defect‐pair was not generated at the end of the cascade event, a higher initial kinetic energy with an energy was given. The initial kinetic energy that ensured the formation of a stable Frankel pair was the threshold displacement energy of the PKA atom.AcknowledgementsC.S. and L.J. contributed equally to this work. This work was supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 12174060 and 12204170, Science Technology Commission of Shanghai Municipality (Explorer project, 21TS1401000), National Key Research and Development Program of China (Nos. 2019YFE0118100 and 2022YFA1404603), Program for Professor of Special Appointment (Eastern Scholar TP2019019), and Project of MOE Innovation Platform.Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.I. Gorczyca, A. Svane, N. E. Christensen, Phys. Rev. B 1999, 60, 8147.T. Roy, E. X. Zhang, Y. S. Puzyrev, D. M. Fleetwood, R. D. Schrimpf, B. K. Choi, A. B. Hmelo, S. T. Pantelides, IEEE Trans. Nucl. Sci. 2010, 57, 3060.L. Lv, J. G. Ma, Y. R. Cao, J. C. Zhang, W. Zhang, L. Li, S. R. Xu, X. H. Ma, X. T. Ren, Y. Hao, Microelectron. Reliab. 2011, 51, 2168.J. Chen, Y. S. Puzyrev, C. X. Zhang, E. X. Zhang, M. W. McCurdy, D. M. Fleetwood, R. D. Schrimpf, S. T. Pantelides, S. W. Kaun, E. C. H. Kyle, J. S. Speck, IEEE Trans. Nucl. Sci. 2013, 60, 4080.M. P. Khanal, S. Uprety, V. Mirkhani, S. Wang, K. Yapabandara, E. Hassani, T. Isaacs‐Smith, A. C. Ahyi, M. J. Bozack, T.‐S. Oh, M. Park, J. Appl. Phys. 2018, 124, 215702.D. Jin, J. A. del Alamo, IEEE Trans. Electron Devices 2013, 60, 3190.A. P. Karmarkar, B. D. White, D. Buttari, D. M. Fleetwood, R. D. Schrimpf, R. A. Weller, L. J. Brillson, U. K. Mishra, IEEE Trans. Nucl. Sci. 2005, 52, 2239.A. Ionascut‐Nedelcescu, C. Carlone, A. Houdayer, H. J. von Bardeleben, J. L. Cantin, S. Raymond, IEEE Trans. Nucl. Sci. 2002, 49, 2733.D. C. Look, D. C. Reynolds, J. W. Hemsky, J. R. Sizelove, R. L. Jones, R. J. Molnar, Phys. Rev. Lett. 1997, 79, 2273.J. Nord, K. Nordlund, J. Keinonen, K. Albe, Nucl. Instrum. Methods Phys. Res., Sect. B 2003, 202, 93.H. He, C. He, J. Zhang, W. Liao, H. Zang, Y. Li, W. Liu, Nucl. Eng. Technol. 2020, 52, 1537.H. Y. Xiao, F. Gao, X. T. Zu, W. J. Weber, J. Appl. Phys. 2009, 105, 123527.H. Wang, X. Guo, L. Zhang, H. Wang, J. Xue, Appl. Phys. Lett. 2019, 114, 244101.H. Wang, G. Xun, J. X. Xue, Deep‐learning interatomic potential for irradiation damage simulations in MoS2 with ab initial accuracy, 2020. https://doi.org/10.48550/arXiv.2010.09547C. Shang, Z. P. Liu, J. Chem. Theory Comput. 2013, 9, 1838.Y.‐F. Li, Z.‐P. Liu, Phys. Rev. Lett. 2022, 128, 226102.S. Chen, S. Ma, Z. P. Liu, J. Phys. Chem. Lett. 2021, 12, 3328.K. Shimizu, Y. Dou, E. F. Arguelles, T. Moriya, E. Minamitani, S. Watanabe, Phys. Rev. B 2022, 106, 054108.A. Stukowski, Mater. Sci. Eng., R 2010, 18, 015012.S. Ahn, C. Dong, W. Zhu, B. J. Kim, Y. H. Hwang, F. Ren, S. J. Pearton, G. Yang, J. Kim, I. Kravchenko, ECS Trans. 2015, 69, 129.E. Patrick, M. E. Law, L. Liu, C. Velez Cuervo, Y. Xi, F. Ren, S. J. Pearton, IEEE Trans. Nucl. Sci. 2013, 60, 4103.H.‐Y. Kim, J. Kim, L. Liu, C.‐F. Lo, F. Ren, S. J. Pearton, J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom. 2012, 30, 012202.S. D. Huang, C. Shang, P. L. Kang, X. J. Zhang, Z. P. Liu, WIREs Comput. Mol. Sci. 2019, 9, e1415.J. Behler, M. Parrinello, Phys. Rev. Lett. 2007, 98, 146401.S. D. Huang, C. Shang, P. L. Kang, Z. P. Liu, Chem. Sci. 2018, 9, 8644.B. Jiang, H. Guo, J. Chem. Phys. 2019, 150, 180901.S. Ioffe, C. Szegedy, Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift, Proceedings of the 32nd International Conference on International Conference on Machine Learning, 2015, 37, 448.D. C. Liu, J. Nocedal, Math. Program. 1989, 45, 503.G. Kresse, J. Furthmuller, Phys. Rev. B: Condens. Matter. 1996, 54, 11169.W. Kohn, L. J. Sham, Phys. Rev. 1965, 140, A1133.G. Kresse, D. Joubert, Phys. Rev. B 1999, 59, 1758.P. E. Blochl, Phys. Rev. B: Condens. Matter 1994, 50, 17953.J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865.C. G. Van de Walle, J. Neugebauer, J. Appl. Phys. 2004, 95, 3851.M. Huang, Z. Zheng, Z. Dai, X. Guo, S. Wang, L. Jiang, J. Wei, S. Chen, J. Semiconduct. 2022, 43, 042101.S. Plimpton, J. Comput. Phys. 1995, 117, 1. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advanced Electronic Materials Wiley

Neural‐Network Potential Simulation of Defect Formation Induced by Knock‐On Irradiation Damage in GaN

Download PDF
8 pages

Loading next page...
 
/lp/wiley/neural-network-potential-simulation-of-defect-formation-induced-by-UjQ1uA0KSP

References (17)

Publisher
Wiley
Copyright
© 2023 Wiley‐VCH GmbH
eISSN
2199-160X
DOI
10.1002/aelm.202300158
Publisher site
See Article on Publisher Site

Abstract

IntroductionGallium nitride is widely used in optoelectronic and electronic devices, such as high electron mobility transistors , light emitting diodes (LEDs), laser diodes, and photodetectors, mainly because it has a wide direct bandgap (3.4 eV) and high electron mobility (2000 cm2 V−1 s−1).[1–4] It also exhibits excellent resistance against irradiations, making it widely used in aerospace and military fields for irradiation‐resistant devices. Its irradiation resistance comes from the strong ionic‐covalent bonding, great thermal stability, and high displacement energy.[5] For both radiation‐resistant and electronic/optoelectronic applications, irradiation damage not only happens in the radiative environment, but also commonly exits in the ion implantation process during device fabrication.[6] The defects caused by irradiation damage greatly impact the device performance, degradation, and reliability. Therefore, understanding the microscopic mechanism of the irradiation damage of GaN is of great importance.Due to the complexity and the picosecond‐scale of the irradiation damage process, it is rather difficult to directly locate the defects and monitor the evolution of the defects using experimental techniques.[7] As a result, the threshold displacement energy (denoted as Ed), which means the least kinetic energy required for an atom to make displacement happen, can only be measured indirectly. For example, by observing the luminescence of GaN LEDs under irradiation, Ionascut‐Nedelcescu et al.[8] determined the Ed of the Ga atoms was determined to be 19±2 eV based on the assumption that the change of luminescence is mainly caused by the displacement of Ga atoms. On the contrary, the investigation of the microscopic mechanism of the irradiation damage mostly relies on the theoretical simulations. For example, Look et al.[9] analyzed the transport characteristics of GaN thin films irradiated by electron beam and obtained the Ed of GaN through fitting the relativistic cross section of atomic displacement as the function of electron energy (Ed (Ga) = 20.5 eV, Ed (N) = 10.8 eV). A more common theoretical approach to utilize molecular dynamics (MD), including the classical MD and the ab initio MD (AIMD), to simulate the atomic displacement process and determine the threshold displacement energy. Nord et al.[10] proposed an empirical Ziegler–Biersack–Litmark (ZBL) shielding potential of GaN and study the irradiation damage, however, the minimal and the average Ed of Ga atoms are determined to be lower than N atoms. He et al.[11] improved the potential by including the short‐range interaction, and the average Ed of Ga is determined to be 45 eV. Interestingly, using ab initio MD, Xiao et al.[12] also evaluated the average Ed of Ga and N as 73.2 and 32.4 eV, indicating the discrepancies between classical MD and AIMD.Conventionally, AIMD simulations are regarded to be more accurate and reliable than classical MD, because the empirical potential cannot accurately describe the saddle states. However, first‐principles calculations are slower and typically limited to hundreds of atoms, preventing long or large‐scale simulations as in classical MD simulations. Therefore, it is of great significance if one can construct a classical MD potential that is able to provides similar accuracy as AIMD. In recent years, significant progress has been achieved in the construction of the potentials with the help from artificial intelligence and machine learning, such as the potentials for Al, MoS2, SiO2, ZnO, etc. Wang et al.[13] constructed a potential for Al based on a hybrid scheme that interpolates the ZBL screened nuclear repulsion potential with a deep learning potential (DP) energy model. Using this potential, they not only obtained the accurate threshold displacement energy of Al, but also simulated the high‐energy radiation damage processes using a supercell containing 600 000 atoms. They also used this scheme to successfully develop the potential of MoS2.[14] However, the DP model is constructed based on the structures from the MD trajectory at different temperatures, making it difficult to overcome the high reaction barriers at low temperatures and causing the system gets trapped in high‐entropy structures at high temperatures. To solve this problem, Liu et al. developed the stochastic surface walking (SSW) algorithm for the structure prediction and pathway searching.[15] It is an unbiased global optimization method which can simultaneously study the minimum and saddle point on potential energy surface (PES). Based on SSW, they proposed a potential training scheme for potential called stochastic surface walking global exploration with global neural network potential (SSW‐NN), which has been successfully applied to SiO2[16] and ZnO.[17] The potentials constructed using this approach are also called global neural network (G‐NN) potentials. Compared with other optimization methods, SSW algorithm can cover more energy surface and explore unbiasedly both minima and saddle points on PES.For GaN, Shimizu et al.[18] proposed a NN potential to calculate the phonon modes and thermodynamics properties of large‐scale defective structures, in which the atoms are assumed to move around the equilibrium positions. However, it is not designed for study the irradiation damage processes, and the atoms are potentially knocked much further away from their equilibrium position. Therefore, in this work, we use the SSW‐NN method to construct the Ga‐N binary G‐NN potential that can be applied to simulate the radiation damage process. By comparing the MD simulations using our G‐NN potential with ab initio calculations, our potential is benchmarked to produce similar accuracy in multiple dimensions. The energies and formation energies of different defect configurations are in good agreement with the ab initio results, not only for the regular intrinsic point defects but also for the defect‐pairs. Such potential is also applied to simulate the knock‐on process of the irradiation damage, and the trajectory and the displacement threshold energy also in good agreement with the ab initio calculations. The global neural network (G‐NN) potential constructed in this work can elevate the accuracy and the efficiency of the MD simulations of the irradiation damage of GaN.Results and DiscussionEvaluation of the G‐NN PotentialThe data set initially contains the elemental phase of Ga and N2, as well as the GaN configurations in which Ga:N ratio is kept as 1:1 (GanNn). After 100 000 steps of training based on this data set, a potential with energy root mean square error (RMS) of 1.056 meV per atom and force RMS of 0.036 eV Å−1 is created. We selected 10 GaN structures, 5 Ga elemental structures, and 5 N2 molecular structures to benchmark this potential, resulting an average energy difference of 0.5 meV per atom from density functional theory (DFT) results. To improve the accuracy of our potential under different chemical environment, we increase the data set configurations with different stoichiometry, varying the Ga:N ratio from 0.06 to 0.94. The configurations are also generated using SSW‐NN approach, and the final training data set consists of 61 414 structures (briefly summarized in Table S1 of the Supporting Information). The RMS errors of the energy and the force of our final G‐NN potential are 5.373 meV per atom and 0.090 eV Å−1, respectively, indicating that our potential can accurately describe the defect properties as well as the irradiation process.Benchmark of Defect Properties Using the G‐NN PotentialIn the binary compound GaN, there exist 6 types of intrinsic point defects, including two vacancies (VGa, VN), two interstitials (Gai, Ni) and two antisites (GaN, NGa). To evaluate the accuracy of our potential on defect properties of GaN, we first calculate the defect formation energies for these 6 intrinsic point defects using our G‐NN potential, and compare them with the DFT calculations, as shown in Figure 1. To illustrate the improvement of our potential by adding the nonstoichiometric structures into the training data set, we also plot the formation energies calculated using the potential purely trained from the initial data set that only contains the stoichiometric configurations. Although the potential trained from the stoichiometric configurations can accurately describe the bulk structure and energy of GaN crystal, it shows significant error on the defect formation energies as compared to the DFT results, because the point defects have nonstoichiometric chemical environment. On the contrary, our final G‐NN potential provides formation energies with errors lower than 0.1 eV for five point defects, except that of Ga vacancy is calculated to be 0.28 eV higher than the DFT result. The agreement between our results using G‐NN potential and DFT result indicates that our potential accurately describe not only the bulk structure of GaN but also the point defects. In addition, we also use the potential to optimized the atomic structures of the point defects, and the atomic positions also matches the DFT calculations.1FigureCalculated single defect formation energies based on the NN potential and DFT calculations (under Ga‐rich condition).In order to demonstrate the robustness of the potential, our G‐NN potential is further used to predict the formation energies of different defect‐pairs in GaN. Usually, the formation energies of these defect‐pairs are high, and their concentrations in GaN samples prepared by conventional methods are relatively low. However, in the irradiation damage, these defect‐pairs may exist in large amount. In fact, a series of defects are generated along the collision cascade trajectory during the irradiation damage process, leading to more complicated chemical environment than the simple intrinsic point defects under thermal equilibrium. Therefore, to prove the robustness of our G‐NN potential, it is necessary to evaluate its performance on more complicated situations, such as the defect‐pairs. We utilize the G‐NN potential to predict the formation energies of 20 different defect‐pairs with the optimized structures, and compare them with the DFT results (Figure 2a). The formation energies predicted by G‐NN potential are generally in good agreement with the DFT results with errors lower than 0.5 eV.2Figurea) Comparison of the formation energies of 20 defect‐pairs calculated using the G‐NN and the DFT (under Ga‐rich condition). Their structures are all optimized. b) Comparison of the formation energies of 118 defect‐pairs calculated with the G‐NN and the DFT. c) The error of 118 defect‐pairs formation energies calculated from the difference between NN and DFT (in Ga‐rich GaN). We use 6 colors and shapes to mark different configurations of one type of defect‐pair and their errors increase in the order from 1 to 6.To further evaluate the accuracy, we also change the relative distance between the defects in a pair by 5 times for each pair, resulting a total 100 more defect‐pair configurations. A summary of the difference between G‐NN potential results and DFT results for the total 120 defect‐pair configurations is shown in Figure 2b,c. Most defect‐pairs have errors below 0.5 eV, and only three configurations have errors above 1 eV, indicating the robustness of our G‐NN potential on handling the defect‐pairs. The defect‐pairs with the largest error are VGa−VGa and VGa−NGa, which can be attributed to two possible reasons. On one hand, the G‐NN potential prediction of VGa already has relatively more significant error than the others, so it is reasonable that the prediction of VGa−VGa and VGa−NGa pairs suffer bigger error. On the other hand, we observe that the error generally grows with the values of formation energies, and the formation energies of these two defect‐pairs are both higher than the others. This behavior is because most of the areas explored by SSW are still on the lower potential energy surface, leading to more significant error when describing the structures with higher energy. The energy and formation energies of these structures with poor prediction are listed in Table S2 (Supporting Information). Except these structures, our potential can accurately describe the properties of various point defects and defect‐pairs in GaN, enabling us to use the potential to optimize different defect configurations and simulate the defect evolution process in a larger system.In addition, we also examine the accuracy of our G‐NN potential for a larger supercell size. Based on a 360 atom supercell, we calculate the energy of the supercell and the formation energy of VGa, and the results are in good agreement with the DFT calculations, meaning that no extra error is introduced with the growth of the model supercell. More details of the supercell structure energies and defect formation energies are listed in Table S3 (Supporting Information). Due to higher computational cost of DFT calculations, we do not compare the defect formation energies obtained using G‐NN potential with DFT results for larger supercells, however, the presented results prove that our G‐NN potential can be used to treat the defect properties of much larger systems almost as accurately as DFT.It is worth discussing that we do not consider the charged defect states in our potential, which is different from the NN potential proposed by Shimizu et al.[18] As mentioned, the potential from Shimizu et al. is designed for calculating the phonon modes and thermodynamics properties, where the atoms are assumed to move around the equilibrium positions under thermodynamic equilibrium, and the charged states are relatively important. On the contrary, the G‐NN potential in this study is used for simulating the irradiation damage processes, in which the primary knock‐on atom (PKA) is given high kinetic energy causing cascade collisions. As a result, the atoms are knocked much further away from their equilibrium position. In such cases, the difference in the collision simulations caused by the different charged state should not be significant. One evidence is that the atomic structures of neutral and charged defects are very similar, compared to the atomic displacements in cascade collisions. However, it does not mean the charged defects are not important, and it is definitely worth studying if the LASP package can support the potential fitting for charged systems in the future.Simulation of Irradiation Damage Process Using the G‐NN PotentialOur G‐NN potential is subsequently used to simulate the irradiation damage process. We first calculate the displacement threshold energy Ed of GaN for several main crystallographic directions using G‐NN potential, and they are compared to the Ed obtained using AIMD as well as the results from Xiao et al. using DFT (Table S4, Supporting Information). As shown in Table S4 of the Supporting Information, the MD simulations using G‐NN potential can provide predictions of Ed as accurate as the ab initio MD simulations. The DFT‐calculated Ed is about 39 eV for Ga and 17 eV for N along b direction, which is the minimal Ed among all crystallographic directions, and the G‐NN potential predicts the Ed as about 40 eV for Ga and 17 eV for N, agreeing with the DFT calculations. As comparison, the classical potentials fail to provide accurate simulations of the irradiation damage processes. Specifically, instead of collision, the PKA atoms directly pass through other atoms for initial kinetic energy higher than 120 eV, if the Tersoff potential and Tersoff‐ZBL potential are used.In order to validate that the generated G‐NN potential can describe the irradiation damage processes with AIMD accuracy, we further simulate the knock‐on displacement processes of Ga/N primary knock‐on atoms. We give the initial kinetic energy of Ed in the a and b directions to the Ga/N PKAs, respectively. Figure 3 shows the potential energy evolutions of the knock‐on processes. For all processes, the potential energies along the predicted using our G‐NN potential agree well with the AIMD trajectories of the potential energies, indicating that our potential is capable of simulating the knock‐on process with AIMD‐level accuracy. In addition, using the atomic configurations along the trajectories of the AIMD simulations, we also use our G‐NN potential to calculate the potential energies, and our G‐NN potential almost identically reproduces the potential energies of the AIMD results.3FigureThe evolutions of the potential energy during the displacements of a) 19 eV and b) 17 eV N PKA, as well as c) 73 eV and d) 40 eV Ga PKA.Furthermore, we also examine the atomic positions during the simulations of cascading collisions using G‐NN potential and AIMD. As shown in Figure 4, we plot the atomic configurations around the PKA atom, which is a Ga PKA with initial 40 eV kinetic energy along the [1¯010]$\bar{1}010]$ direction (b direction), at 50, 150, and 250 fs. A VGa and an octahedral Gai are formed in both 150 and 250 fs simulations, and the atoms are in relatively similar positions, indicating that the G‐NN potential can be used to simulate the cascade collisions during the irradiation damage process with ab initio accuracy.4FigureThe trajectories of a Ga 40 eV PKA along [1¯010$\bar{1}010$] direction simulated by a–d) DFT and e–h) NN. The red and blue balls represent Ga PKA and Ga interstitial, respectively, and the green arrow represents the direction of Ga PKA. The two axes a and b represent [1¯21¯0]$\bar{1}2\bar{1}0]$ and [1¯010$\bar{1}010$] direction, respectively.Last but not the least, we apply the G‐NN potential to carry out a relatively large‐scale simulations that are conventionally beyond the capability of AIMD. A 2 3040 atom supercell is used to simulate the knock‐on displacement processes of 500 eV Ga PKA along 11 different directions, and each calculation is carried out for 1000 fs. By monitoring the defect types and their populations, we observe that the number of induced defects first increases before it reaches a peak value, and decreases after due to the recovery of defects. The Wigner–Seitz cell method[19] is employed to identify vacancies, interstitials and antisite defects. A cartoon of the defect evolution is shown in Figure S1 of the Supporting Information, and only displaced atoms or defects are plotted. The peak damage is found to happen at 0.24 ps, and the number of defects keep decreasing afterward. Based on the statistics of the defect types and their populations, the relative ratios of all six intrinsic defects are determined by averaging the numbers of defects for all directions. The defect ratios for the peak damage and the final state are both listed in Table S6 of the Supporting Information. For all directions, more Ga and N vacancies, as well as N interstitials are generated than the Ga interstitials and antisites. This finding contrasts with the Stopping and Range of Ions in Matter (SRIM) and technology computer‐aided design (TCAD) simulations of GaN devices in which only Ga vacancies are considered.[20–22] To confirm our results, a larger supercell of 184 320 atoms is also used to repeat the knock‐on displacement process simulations, and the results agree with those from the 23 040 atom supercell. It is worth noting that our G‐NN potential not only enables the MD simulations of large‐scale systems with accuracy that is comparable to ab initio simulations, but also provide parameters for device simulations. The displacement threshold energy and the information of the relative defect populations determined in MD simulations can be utilized as the input parameter of the TCAD simulations of the proton irradiation damage of GaN‐based devices, such as AlGaN/GaN HEMTs.ConclusionsIn summary, we utilize machine learning algorithms and global PES exploration to develop a G‐NN potential for GaN. Based on a dataset of 61 414 structures, the potential is trained to reduce the RMS errors of energy and force as 5.373 meV per atom and 0.090 eV Å−1. The formation energies of the point defects and defect‐pairs calculated using G‐NN potential are in good agreement with the results of DFT, indicating that the potential is capable of accurately describing the defect properties. In addition, using the G‐NN potential, we also investigate the threshold displacement energies of PKA atoms along different directions, the trajectories and the potential energies of several knock‐on processes, all of which are in good agreement with AIMD calculations. Finally, based on the simulations of the knock‐on processes of a Ga PKA atom, the relative ratios of the induced intrinsic point defects are obtained. Ga and N vacancies, as well N interstitials are found to have higher amounts than the Ga interstitials and the antisites. The G‐NN potential extends the scales of time and size of the atomic‐level simulations of the irradiation damage processes while keeping similar accuracy as ab initio approach. Therefore, it is not only able to provide the fundamental understanding of the irradiation damage processes, but also to refine the parameters such as the threshold displacement energy in the device simulations.Experimental SectionSSW‐NN MethodThe constructions of the G‐NN potential were conducted using LASP code[23] (Large‐scale Atomic Simulation with neural network Potential) developed by Liu et al. group. The random stochastic surface walking (SSW) method[15] was utilized for the global optimization, which can automatically explore the multidimensional PES and identify important structures and reaction paths on the potential energy surface. The G‐NN potential follows the atom‐centered NN potential architecture proposed by Behler and Parrinello,[24] where the input layer utilized the power‐type structure descriptor (PTSD).[25] The G‐NN potential is sampled by SSW global potential energy surface, fitted and iteratively generated by data obtained from high‐precision DFT calculation, until the G‐NN potential has sufficient robustness to quantitatively describe global PES, or namely the SSW‐NN.[26]The training process of the potential divided into three main stages as illustrated in Figure 5. The first stage was to construct an initial dataset containing the most common atomic environment of the system as a training set for building the initial global PES. In this stage, the dataset only includes structures with Ga:N ratio as 1:1, and the elemental phase of Ga and N2. The random stochastic surface walking (SSW) was used to automatically explore and obtain the multidimensional potential energy surface. In total, 247 246 structures were generated, and 6030 of them were randomly selected into the dataset and simulated using first‐principles method. These structures collected and screened from the SSW trajectory were calculated using first‐principles method and added into the dataset. Eventually a dataset of over 6000 configurations is built.5FigureThe schematic plot of the construction process of the G‐NN potential.In the second stage, based on the dataset obtained in the first stage, LASP was used to carry out the initial potential training. In order to ensure the accuracy of the potential, a large number of PTSDs, including 56 two‐body, 60 three‐body, and 8 four‐body descriptors was used. The general forms of the PTSD are described in the following equations1fcrij=0.5×tanh31−rijrc,forrij≤rc0forrij>rc$$\begin{eqnarray}{f}_{\rm{c}}\ \left( {{r}_{{\rm{ij}}}} \right) = \left\{ { \def\eqcellsep{&}\begin{array}{@{}*{1}{c}@{}} {{\rm{0}}{\rm{.5}} \times \tan {{\rm{h}}}^3\left( {1 - \frac{{{r}_{{\rm{ij}}}}}{{{r}_{\rm{c}}}}} \right),for\ {r}_{{\rm{ij}}} \le {r}_{\rm{c}}}\\[6pt] {\ 0\ for\ {r}_{{\rm{ij}}} > {r}_{\rm{c}}} \end{array} } \right.\ \end{eqnarray}$$2Si1=∑j≠iRnrij$$\begin{eqnarray}S_{\rm{i}}^1 = \mathop \sum \limits_{j \ne i} {R}^n\ \left( {{r}_{{\rm{ij}}}} \right)\end{eqnarray}$$3Si2=∑m=−LL∑j≠iRnrijYLmrij212$$\begin{eqnarray}S_{\rm{i}}^2 = {\left[ {\mathop \sum \limits_{{\rm{m}} = - L}^L {{\left| {\mathop \sum \limits_{j \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right){Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right)} \right|}}^2} \right]}^{\frac{1}{2}}\ \end{eqnarray}$$4Si3=21−ζ∑j1+λcosθijkζ·Rnrij·Rmrik·Rprik$$\begin{eqnarray}S_{\rm{i}}^3 = {2}^{1 - \zeta }\ \mathop \sum \limits_j {\left( {1 + \lambda \cos {\theta }_{{\rm{ijk}}}} \right)}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{ik}}}} \right) \cdot {R}^p\left( {{r}_{{\rm{ik}}}} \right)\end{eqnarray}$$5Si4=21−ζ∑j,k,l≠i[1+λcosθijkζ·Rnrij·Rmrjk$$\begin{eqnarray}S_i^4 = {2}^{1 - \zeta }\ \mathop \sum \limits_{j,k,l \ne i} [{\left( {1 + \lambda \cos {\theta }_{{\rm{ijk}}}} \right)}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{jk}}}} \right)\end{eqnarray}$$6Si5=∑m=−LL∑j,k≠iRnrij·Rmrik·Rprjk·YLmrij+YLmrik212$$\begin{eqnarray}S_{\rm{i}}^5 = {\left. {{{\left\{ {\mathop \sum \limits_{m = - {\rm{L}}}^{\rm{L}} \left| {\mathop \sum \limits_{j,k \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right) \cdot {R}^m\left( {{r}_{{\rm{ik}}}} \right) \cdot {R}^p\left( {{r}_{{\rm{jk}}}} \right) \cdot \left[ {{Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right) + {Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ik}}}} \right)} \right]} \right|} \right.}}^2} \right\}}^{\frac{1}{2}}\ \end{eqnarray}$$7Si6=21−ζ∑j,k,l≠i1+λcosδijklζ·RnrijRmrikRpril$$\begin{eqnarray}S_{\rm{i}}^6 = {2}^{1 - \zeta }\ \mathop \sum \limits_{j,k,l \ne i} \left[ {{{\left( {1 + \lambda \cos {\delta }_{{\rm{ijkl}}}} \right)}}^\zeta \cdot {R}^n\left( {{r}_{{\rm{ij}}}} \right){R}^m\left( {{r}_{{\rm{ik}}}} \right){R}^p\left( {{r}_{{\rm{il}}}} \right)} \right]\end{eqnarray}$$8Si7=∑j≠iRnrij·1+tanhα·rij−rm$$\begin{eqnarray}S_{\rm{i}}^7 = \mathop \sum \limits_{j \ne i} \left\{ {{R}^n} \right.\ \left. {\left( {{r}_{{\rm{ij}}}} \right) \cdot \left[ {1 + \tan {\rm{h}}\alpha \cdot \left( {{r}_{{\rm{ij}}} - {r}_{\rm{m}}} \right)} \right]} \right\}\end{eqnarray}$$9Si8=∑m=−LL∑j≠iRnrijYLmrij·1+tanhα·rij−rm212$$\begin{eqnarray}S_{\rm{i}}^8 = {\left[ {\mathop \sum \limits_{{\rm{m}} = - {\rm{L}}}^{\rm{L}} {{\left| {\mathop \sum \limits_{j \ne i} {R}^n\left( {{r}_{{\rm{ij}}}} \right){Y}_{{\rm{Lm}}}\left( {{r}_{{\rm{ij}}}} \right) \cdot \left[ {1 + \tan {\rm{h}}\left[ {\alpha \cdot \left( {{r}_{{\rm{ij}}} - {r}_{\rm{m}}} \right)} \right]} \right]} \right|}}^2} \right]}^{\frac{1}{2}}\ \end{eqnarray}$$whererijis the internuclear distance between atomiandj,θijkis the angle centered at iatom withjand kbeing neighbors (i,jand kare atom indices). The key factors of PTSD are the cut‐off functionfc, which is decayed to zero beyondrc, power‐type radial function, trigonometric angular functions, and spherical harmonic function. In PTSD, theS1andS2are two‐body functions, theS3,S4, andS5are three‐body functions, and theS6is a four‐body function.The G‐NN potential has a five‐layer (128‐80‐80‐48‐1) feed‐forward NN structure, in total 40 834 fitting parameters. The min‐max scaling method[27] was used to normalize the training dataset. Hyperbolic tangent function was used to activate the hidden layer, while linear transformation was performed on the output layer of all networks. The quasi‐Newton Limited‐memory Broyden–Fletcher–Goldfarb–Shanno (L‐BFGS)[28] method was used to minimize the loss function to match the energy, force and stress from the DFT calculations. In the process of potential training, the number of training steps set to 50 000, and the training weights of energy, force and stress set to 200, 15, and 1, respectively.In the third stage, the global dataset is expanded iteratively to improve the predictive power of the neural network potential. The SSW‐NN method was still adopted as in the second stage, however, the structures found by SSW‐NN are calculated using the constructed potential and classical MD, which is different from the DFT calculations used in the first stage. Structures collected and filtered from the SSW trajectories were used to expand the dataset. The added configurations included 45 categories of nonstoichiometric structures with Ga/N ratio ranging from 1:15 to 15:1, and the numbers of atoms in a unit cell were in the range of 12–128 atoms. The purpose of the potential is focused on the point defect, which is localized in a small area of the whole crystal. As a result, small nonstoichiometric unit cells are good enough to describe the chemical environment of point defects. Adding a few large supercells is necessary to utilize all structure descriptors in the potential training process adequately. Such rule of constructing the data set not only enables more accurate defect descriptions but also accelerates the training process of potential. Then the new dataset was input to the global dataset to start a new neural network training cycle. The G‐NN potential is obtained until an additional training cycle does not obviously improve the accuracy potential.DFT CalculationsAll density functional theory calculations were performed by using VASP[29] (Vienna Ab‐initio Simulation Package) . The projector‐augmented wave method[30–32] and a 520 eV energy cutoff were used for all calculations. The generalized gradient approximation (GGA) exchange‐correlation functional in the form of Perdew–Burkes–Ernzerhof (PBE)[33] was used. The Monkhorst‐Pack k‐point grid was adopted, and the density of the grid was dependent on the size of the first Brillouin zone. The energy convergence criterion was set as 1 × 10−5 eV, and the lattice was optimized until the force on each atom is less than 0.01 eV Å−1.For the calculation of GaN defects, we only consider defects at neutral charge state. A quasicubic 128 atom supercell is selected to calculate the defect properties. The defect formation energy is calculated using[34,35]10ΔHfa,q=Ea,q−Ebulk+∑niμi+qEF$$\begin{eqnarray}\Delta {H}_{\rm{f}}\ \left( {a,q} \right) = \ E\left( {a,q} \right) - E\left( {{\rm{bulk}}} \right) + \sum {n}_{\rm{i}}{\mu }_{\rm{i}} + q{E}_{\rm{F}}\end{eqnarray}$$where E(a, q) and E(bulk) represent the total energies of supercells with and without defect, respectively. ni is the number of atoms that the supercell with the defect differs from the supercell without defect. µi represents the chemical potential of the corresponding element, and this chemical potential is changed with reference to the total energy Ei of the elemental phase of the element. EF denotes the Fermi level. For a certain chemical potential condition, the defect formation energy of electric neutral is a constant for Fermi level in the whole range of the bandgap. For the convenience of later comparison, the formation energy calculated later is obtained under the condition of Ga‐rich.MD SimulationsAll MD simulations were carried out using the LAMMPS (Large‐scale Atomic/Molecular Massively Parallel Simulator) package,[36] which has been widely used to simulate the irradiation damage process of materials including metals and semiconductors. The atomic visualizations were created using the Open Visualization Tool (OVITO).[19] All simulations were based on a 2880 atom supercell with lattice constants of 32, 33, and 31 Å, and periodic boundary conditions were imposed in all directions. The NVT ensemble was used to relax the crystal for 1000 timesteps at 100 K with time step of 1 ps. Once the system reached equilibrium, a Ga atom or N atom was selected as the PKA and given initial kinetic energy to simulate cascade collisions based on NVE ensemble. The whole process lasted 5000 timesteps, and the time step adaptively changed to avoid instability of simulation results. If a stable Frankel defect‐pair was not generated at the end of the cascade event, a higher initial kinetic energy with an energy was given. The initial kinetic energy that ensured the formation of a stable Frankel pair was the threshold displacement energy of the PKA atom.AcknowledgementsC.S. and L.J. contributed equally to this work. This work was supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 12174060 and 12204170, Science Technology Commission of Shanghai Municipality (Explorer project, 21TS1401000), National Key Research and Development Program of China (Nos. 2019YFE0118100 and 2022YFA1404603), Program for Professor of Special Appointment (Eastern Scholar TP2019019), and Project of MOE Innovation Platform.Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.I. Gorczyca, A. Svane, N. E. Christensen, Phys. Rev. B 1999, 60, 8147.T. Roy, E. X. Zhang, Y. S. Puzyrev, D. M. Fleetwood, R. D. Schrimpf, B. K. Choi, A. B. Hmelo, S. T. Pantelides, IEEE Trans. Nucl. Sci. 2010, 57, 3060.L. Lv, J. G. Ma, Y. R. Cao, J. C. Zhang, W. Zhang, L. Li, S. R. Xu, X. H. Ma, X. T. Ren, Y. Hao, Microelectron. Reliab. 2011, 51, 2168.J. Chen, Y. S. Puzyrev, C. X. Zhang, E. X. Zhang, M. W. McCurdy, D. M. Fleetwood, R. D. Schrimpf, S. T. Pantelides, S. W. Kaun, E. C. H. Kyle, J. S. Speck, IEEE Trans. Nucl. Sci. 2013, 60, 4080.M. P. Khanal, S. Uprety, V. Mirkhani, S. Wang, K. Yapabandara, E. Hassani, T. Isaacs‐Smith, A. C. Ahyi, M. J. Bozack, T.‐S. Oh, M. Park, J. Appl. Phys. 2018, 124, 215702.D. Jin, J. A. del Alamo, IEEE Trans. Electron Devices 2013, 60, 3190.A. P. Karmarkar, B. D. White, D. Buttari, D. M. Fleetwood, R. D. Schrimpf, R. A. Weller, L. J. Brillson, U. K. Mishra, IEEE Trans. Nucl. Sci. 2005, 52, 2239.A. Ionascut‐Nedelcescu, C. Carlone, A. Houdayer, H. J. von Bardeleben, J. L. Cantin, S. Raymond, IEEE Trans. Nucl. Sci. 2002, 49, 2733.D. C. Look, D. C. Reynolds, J. W. Hemsky, J. R. Sizelove, R. L. Jones, R. J. Molnar, Phys. Rev. Lett. 1997, 79, 2273.J. Nord, K. Nordlund, J. Keinonen, K. Albe, Nucl. Instrum. Methods Phys. Res., Sect. B 2003, 202, 93.H. He, C. He, J. Zhang, W. Liao, H. Zang, Y. Li, W. Liu, Nucl. Eng. Technol. 2020, 52, 1537.H. Y. Xiao, F. Gao, X. T. Zu, W. J. Weber, J. Appl. Phys. 2009, 105, 123527.H. Wang, X. Guo, L. Zhang, H. Wang, J. Xue, Appl. Phys. Lett. 2019, 114, 244101.H. Wang, G. Xun, J. X. Xue, Deep‐learning interatomic potential for irradiation damage simulations in MoS2 with ab initial accuracy, 2020. https://doi.org/10.48550/arXiv.2010.09547C. Shang, Z. P. Liu, J. Chem. Theory Comput. 2013, 9, 1838.Y.‐F. Li, Z.‐P. Liu, Phys. Rev. Lett. 2022, 128, 226102.S. Chen, S. Ma, Z. P. Liu, J. Phys. Chem. Lett. 2021, 12, 3328.K. Shimizu, Y. Dou, E. F. Arguelles, T. Moriya, E. Minamitani, S. Watanabe, Phys. Rev. B 2022, 106, 054108.A. Stukowski, Mater. Sci. Eng., R 2010, 18, 015012.S. Ahn, C. Dong, W. Zhu, B. J. Kim, Y. H. Hwang, F. Ren, S. J. Pearton, G. Yang, J. Kim, I. Kravchenko, ECS Trans. 2015, 69, 129.E. Patrick, M. E. Law, L. Liu, C. Velez Cuervo, Y. Xi, F. Ren, S. J. Pearton, IEEE Trans. Nucl. Sci. 2013, 60, 4103.H.‐Y. Kim, J. Kim, L. Liu, C.‐F. Lo, F. Ren, S. J. Pearton, J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom. 2012, 30, 012202.S. D. Huang, C. Shang, P. L. Kang, X. J. Zhang, Z. P. Liu, WIREs Comput. Mol. Sci. 2019, 9, e1415.J. Behler, M. Parrinello, Phys. Rev. Lett. 2007, 98, 146401.S. D. Huang, C. Shang, P. L. Kang, Z. P. Liu, Chem. Sci. 2018, 9, 8644.B. Jiang, H. Guo, J. Chem. Phys. 2019, 150, 180901.S. Ioffe, C. Szegedy, Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift, Proceedings of the 32nd International Conference on International Conference on Machine Learning, 2015, 37, 448.D. C. Liu, J. Nocedal, Math. Program. 1989, 45, 503.G. Kresse, J. Furthmuller, Phys. Rev. B: Condens. Matter. 1996, 54, 11169.W. Kohn, L. J. Sham, Phys. Rev. 1965, 140, A1133.G. Kresse, D. Joubert, Phys. Rev. B 1999, 59, 1758.P. E. Blochl, Phys. Rev. B: Condens. Matter 1994, 50, 17953.J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865.C. G. Van de Walle, J. Neugebauer, J. Appl. Phys. 2004, 95, 3851.M. Huang, Z. Zheng, Z. Dai, X. Guo, S. Wang, L. Jiang, J. Wei, S. Chen, J. Semiconduct. 2022, 43, 042101.S. Plimpton, J. Comput. Phys. 1995, 117, 1.

Journal

Advanced Electronic Materials – Wiley

Published: Aug 1, 2023

Keywords: GaN; defect‐pairs; global neural network potential; irradiation damage

There are no references for this article.