Interaction potentials for modelling GaN precipitation and solid state polymorphism

We outline a molecular mechanics model for the interaction of gallium and nitride ions ranging from small complexes to nanoparticles and bulk crystals. While the current GaN force fields allow the modelling of either bulk crystals or single ions dispersed in solution, our model covers both and hence paves the way to describing aggregate formation and crystal growth processes from molecular simulations. The key to this is the use of formal  +3 and  −3 charges on the gallium and nitride ions, whilst accounting for the charge transfer in GaN crystals by means of additional potential energy terms. The latter are fitted against experimental data of GaN in the wurtzite structure and benchmarked for the zinc-blende and rock-salt polymorphs. Comparison to quantum chemical references and experiment shows reasonable agreement of structures and formation energy of [GaN]n aggregates, elastic properties of the bulk crystal, the transition pressure of the wurtzite to rock-salt transformation and intrinsic point defects. Furthermore, we demonstrate force field transferability towards the modelling of GaN nanoparticles from simulated annealing runs.


Introduction
Gallium nitride is a direct band-gap semiconductor with appealing properties thanks to its wide band-gap (~3.2 eV) [1], high thermal conductivity (~2.3 W (cm K) −1 ) [2] and favorable electron mobility (~1500 cm 2 V −1 s −1 ) [2]. Crystals of GaN permit broad applications, ranging from (opto-) electronic devices to gas sensors. This is however hampered by the ongoing challenge of producing high-quality single crystalswhich inspired considerable research efforts from both experiment and theory.
GaN was first synthesized by Maruska and Tietjen in 1969 using a hydride vapor phase epitaxy (HVPE) method [3]. However, its potential use became more apparent when Nakamura developed a metalorganic chemical vapor deposition (MOCVD) route for GaN growth in 1991 [4,5]. Using the buffer-layer technique and developing the idea of double heterostructures, Akasaki, Amano and Nakamura [5,6] demonstrated the first blue light-emitting diode using GaN and its ternary alloy with indium, In x Ga 1−x N [7]. This important contrib ution led to the Nobel Prize in Physics in 2014. Commonly used substrate materials include Si, SiC, or sapphire. Due to inevitable lattice mismatch, the growth on foreign surfaces leads to high dislocation density (typically higher than 10 8 cm −2 ), which results in low performance of the final device [8]. An ideal solution to this problem would be to grow GaN on GaN substrate, hence its native crystal surface. However, the syntheses of sufficiently large GaN single crystals remain an immense challenge. Indeed GaN has high melting temper ature and chemical resistance (it is stable up S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
to ~1200 K at an associated nitrogen pressure of 1 atm or to ~4000 K at 12 GPa) [9,10]-which obscures conventional growth (such as Czochralski and Bridgman) methods.
At present, among the two main routes to GaN crystals (vapor phase and flux growth methods [8]), nucleation from solution and in particular ammonothermal syntheses is proposed as greatly promising for GaN growth with high crystalline quality [11,12]. The underlying process of ammonothermal synthesis is based on dissolving metallic Ga, its amide or polycrystalline GaN in liquid ammonia at supercritical condition in an autoclave. Acidic or basic mineralizers are generally added to the system to increase the solubility of Ga. By means of implementing local cooling, GaN is then gradually grown according to a temperature gradient in autoclave [13]. To date, the understanding of GaN growth by ammonothermal synthesis is still limited. However, a number of studies were subjected to identifying precursors to GaN nucleation. Niewa and coworkers, for example, investigated dissolved intermediates in the ammonobasic GaN crystal growth and detected tetraamidogallate ions [Ga(NH 2 ) 4 ] − as the dominating gallium-containing species towards gallium nitride formation [14].
The need for molecular scale insights has motivated a variety of modelling and simulation studies. Early theoretical studies mainly focused on the gas-phase reaction mechanism of GaN growth by means of ab-initio calculations [15][16][17][18][19]. Goddard and coworkers, for example, explored the reaction pathways of the GaN growth on the GaN (0 0 0 1) surface from various gas precursors that can act as gallium or nitride sources [19]. While quantum mechanical (QM) calculations are well-suited for such comparably simple model systems, the description of solution-based processes calls for large models and extended statistical sampling. For this, molecular mechanics (MM) and combined QM/MM approaches are suggested as a compromise in accurate treating interatomic forces and sufficient sampling to ensure statistical significance. As a key requisite, a small series of MM potentials for GaN crystal have been developed. The underlying force fields were based on either the Stillinger-Weber and Tersoff-Brenner formalisms [20,21], or, more recently, consider mixed covalent-ionic character. For the latter, the potentials include long-range Coulombic terms, based on residual charges on the Ga and N ions. This reflects a considerable advantage in terms of force field flexibility as covalent bonds no longer need to be pre-defined (and hence permanent) throughout the molecular simulation. Along those lines, Zapol et al devised a GaN force field using partial charges of +2 and −2 for Ga and N, respectively, whilst considering electronic polarization by combining Buckingham potentials with a shell model applied to the nitrogen ions [22]. Similarly, Gao et al modified the Buckingham potential to mimic the polarization effects without the shell model [23].
When focusing on bulk GaN, the residual charges may be fitted as free parameters as a part of the force field creation procedure, hence allowing deviation from the formal charges. This, however, significantly limits transferability to small clusters or single ions, for which re-parameterization of the residual charges for each ion complex would be needed. On the other side, force field models suited for ions dispersed in solution must reflect the proper +3 charge of Ga 3+ (along with −1, −2 and −3 for NH − 2 , NH 2− and N 3− , respectively) [24,25]. To pave the way towards the holistic molecular modelling of GaN nucleation and growth, we desire a force field based on formal charges, hence allowing universal modelling of single ions, ionic complexes, nano-scale aggregates and the bulk solid. This is the aim of the present study. In what follows, we develop a hybrid Lennard-Jones/Born-Mayer-type GaN potential for this purpose. Our model is then benchmarked against experimental observables and compared to quantum calculations of [GaN] n nanoclusters. Moreover, to demonstrate the versatility and efficiency of this molecular mechanics approach, crystallization of GaN droplets comprising ~1500 formulae units is illustrated as a function of cooling time.

Methods and models
The molecular mechanics models aim at bulk gallium nitride and [Ga x N y ] aggregates with x ≈ y . Considering the covalent/ ionic interactions, our models were based on a combination of Coulomb, Lennard-Jones (LJ) and Born-Mayer (BM) potentials as given by: where r ij is the distance between atoms i and j . The cutoff radius for short-range non-Coulombic interactions was set as 12 Å. The partial charges q i , q j were fixed to the formal charges of +3 (Ga) and −3 (N) to ensure force field transferability from bulk crystals to clusters and single ions. In turn, all electronic polarization effects are mimicked by the LJ and BM terms in equation (1), by parameterization of ɛ, σ, A, and ρ respectively. We first attempted to parameterize only a pure Buckingham-type force field ( A · exp(−r ij /ρ) − C/r 6 ij ) for describing non-Coulombic interactions in wurtzite GaN. From this, a  figure 3 for benchmarking against all known polymorphs and Ga n N n aggregates, respectively.

Lennard-Jones potential parameters
parameter set that can reasonably reproduce the structural parameters and elastic constants of all known GaN polymorphs (see supporting information (stacks.iop.org/ JPhysCM/32/205401/mmedia)) can be obtained. However, when this model is used for molecular dynamics simulations of GaN nanoparticle and bulk GaN at high temperature (>~1500 K), we found the Ga-N pairs may cross the repulsive part of the BM potential barrier. This is leading to what is sometimes called the Buckingham catastrophe, i.e. the dominance of Coulomb attraction (and also −C/r 6 ij ) in r ij → 0 and thus the artificial fusion of atoms. To avoid this problem, we decided to rely on LJ type potentials for better description of the Ga-Ga and N-N short-range repulsion, whilst BM type terms are used for accurately modelling the Ga-N interactions. The BM potential is only applied to cation-anion interactions, as the Coulombic repulsion   [34]. b Ueno et al [26]. c Polian et al [35]. d From first-principle calculations [36]. e Uehara et al [37]. of the Ga-Ga and N-N interactions readily prevents short distances of cation-cation and anion-anion pairs, respectively. In other words, the contrib ution of short-range terms to Ga-Ga and N-N interaction is model led by the LJ approach, offering the use of the Lorentz-Berthelot mixing rules for the combination with other LJ potentials when studying solvent-GaN or other systems involving further molecular species. The parameter optimization is based on least-square fitting of the weighted sum of deviation from experimental data: Were f denotes the lattice parameters, atomic positions and the elastic constants of bulk GaN in the wurtzite structure, hence leaving the zinc-blende and rock-salt polymorphs for model validation [26]. The weight factors w m were systematically adjusted to achieve a balanced fit of reasonable accuracy for all reference data [27]. The parameterization was done at zero temperature and pressure using the relax fit and genetic algorithm implemented in the GULP package [27]. The density functional theory calculations were performed using the B3LYP functional [28,29] and the LanL2DZ basis set [30,31]. This level of theory, along with the series of Ga n N n aggregates with n = 1..20 was adopted from the recent study of Brena and Ojamäe [32]. Using the Gaussian 09 package [33], all clusters were subjected to geometry optimization to get DFT-based reference structures and formation energies. Likewise, our force field model is used to perform MM-based geometry optimization, allowing for force field validation in terms of both aggregate structure and formation energy.

Results and discussion
Our molecular mechanics model of GaN aims at wide transferability to cover not only the bulk phases, but also nanocrystals and ion aggregates. We therefore decided to only use a limited set of input data for parameterization, whilst reserving other properties and structures for validation. Along those lines, the wurtzite crystal structure, including the elastic constants as known from the experiment, were taken as inputs for force field optimization. The resulting force field parameters are denoted in table 1.
The wurtzite, zinc-blende and rock-salt structures of GaN are illustrated in figure 1. The lattice constants a, c of wurtzite along with the internal coordinate u describing the respective position of Ga and N in the unit cell are all reproduced within error margins below 3.4%. The assessment of the primitive unit cell volume benefits from error compensation and shows good agreement between experimental and calculated volumes (22.87 versus 23.03 Å 3 /f.u.). For the elastic constants, which were also used as inputs to force field optimization, the agreement with the experimental values was driven to the limits of the experimental error margins [23]. For instance, the calculated bulk modulus for wurtzite phase was found as 235.3 GPa, whilst experimental data ranges from 188 GPa [34] to 210 GPa [35] and 237 GPa [26], respectively.
First benchmarks to assessing the quality of our GaN model were derived from analyzing the zinc-blende and the rock-salt structures. The calculated lattice parameters of these two solid polymorphs are in good agreement with the experimental values (table 2). For the bulk modulus of the zinc-blende phase, there is no experimental value to compare with our calculated value of 235.2 GPa. Force field potentials developed by Zapol et al [22] and Gao et al [23] however predicted the bulk modulus as 227 GPa and 222.3 GPa, respectively. For the rock-salt phase, Xia et al [34] fitted the Birch-Murnaghan equation of state to the experimental data to assess the bulk modulus as 248 GPa. Later on, Uehara et al [37] derived the bulk modulus of rock-salt phase as 323 ± 27 GPa. Considering the error margins of high-pressure experiments, this is still in reasonable agreement with our prediction of 320.7 GPa. Overall, these assessments show the validity of our interatomic interaction potential in describing GaN in all relevant crystalline forms, i.e. wurtzite, zinc-blende and rock-salt phases.
At ambient conditions, the wurtzite structure represents the stable form of GaN, whilst zinc-blende and rock-salt are metastable in decreasing order of stability [26]. Much in analogy to previous force field developments [22,23], this is also reproduced by our model. However, at high pressure the more compact rock-salt structure (nearest neighbor coordination number of 6, as compared to 4 for wurtzite and zinc-blende, see also figure 1) becomes increasingly stable. High-pressure anvil cell experiments indicate the wurtzite to rock-salt transformation to occur at a transition pressure between 37 GPa [34] and 52 GPa [26], depending on the experimental setup. To explore the pressure-dependence of solid phase stability, we used our force field to screen pressure from 0 to 70 GPa whilst allowing full structural relaxation. From this, enthalpy profiles as functions of pressure are obtained. As can be seen from figure 2, the wurtzite phase was found as more stable than the zinc-blende structure for the entire pressure range explored (at 0 K). Based on our force field, the transition pressure for the wurtzite to rock-salt transformation is estimated from the intersection of the corresponding enthalpy profiles as 36.5 GPa-which is in excellent agreement with the experiment (37 GPa) from Xia et al [34].
Intrinsic point defects in GaN crystal were investigated by employing the Mott-Littleton approach [38,39]. In this approach, the lattice for energy minimization is divided into two regions: (i) the inner region, which contains a defect or an ensemble of many defects and (ii) the outer region, which extends to infinity. A radius of 10 Å for the inner region (comprising of ~380 ions) and 16 Å for the outer region (comprising of ~1100 ions) were found to achieve a desired convergence. Schottky and Frenkel defects are investigated in the wurtzite and zinc-blende GaN crystal structures. While a previous study considered a defect at infinite dilution [22], we modelled the Schottky defects as nearest-neighbor sites such that the binding energy due to defects aggregation is taken into account. For Frenkel defects, a single ion (Ga 3+ or N 3-) is moved and placed at close proximity (second-nearest interstitial site), thus creating a pair of vacancy and interstitial defects. To avoid immediate recombination, the interstitial was held fixed during the defect calculation.
The Schottky defect formation energies in the wurtzite and zinc-blende phases were found as 4.1 and 3.8 eV/defect, respectively. Based on macroscopic analyses, Van Vechten estimated the formation enthalpy for Schottky defect in GaN real crystal as 3.9 eV/defect [40]. Due to the small difference between the Schottky defect formation energies of wurtzite and zinc-blende phases and for simplicity, in what follows we focus on the zinc-blende phase for investigating Ga and N Frenkel defects, respectively. In the zinc-blende structure, Frenkel defects can be classified based on the interstitial occupying a hexagonal or tetrahedral site (T). Moreover, the T site can be further categorized into a site that is surrounded by four N atoms (denoted as T N ) and T Ga site by four Ga atoms. Among all these possible sites, we found the interstitial at T N site leads to the lowest formation energy for Ga Frenkel defect (6.0 eV/ defect) and T Ga site for N Frenkel defect (5.4 eV/defect). The formation of N Frenkel defect is thus more favorable than that of Ga Frenkel defects. These results are in good agreement with previous theoretical studies by Zapol et al [22] who predicted the formation energies of the Schottky, N Frenkel, Ga Frenkel defects as 4.8, 6.1, and 6.9 eV/defect and by Chisholm et al [41] as 4.74, 6.66, and 7.42 eV/defect, respectively. Moreover, the simple ion vacancy defect energies in both wurtzite and zinc-blende phases were calculated as 48.9 and 48.4 eV for removing a Ga 3+ and N 3− ion, respectively.
A particularly challenging benchmark to force fields derived from fitting to bulk phases is given by the comparison of finite structures and their formation energy. Inspired by the analyses of Wu et al [42] which was focused to Al n N n aggregates (n = 1-20), we performed DFT-based structure optimizations of the homologous series of Ga n N n aggregates. This was compared to Ga n N n aggregate relaxation from our molecular mechanics model. In terms of structure, figure 3 highlights reasonable agreement with the exception of n = 1-3, whereas the matching of force field optimized and DFT reference aggregates improves with increasing size. This is attributed to the (deliberate) lack of coordination-dependent electronic polarizability in the present force field-as has also been observed for ZnS clusters in a recent study by Woodley et al [43].
The trend towards increasing force field accuracy for increasingly large (GaN) n aggregates is also reflected by the cohesive energy, as shown normalized per formula unit (n = 1) in figure 4. For cohesive energy calculation, we use E coh (n) = E(n)/n − E(1); where E(n) is the total energy of the (GaN) n aggregate and n = 1 refers to the Ga-N dimer. The cohesive energy gradually decreases as the cluster size increases (figure 4), in line with the DFT study of Brena and Ojamäe [32] who explored the range n = 6 to 48. Within the explored range of aggregates (n = 1..20) we find significant deviation of the cohesive energy calculated by our force field as compared to the DFT reference for n < 6, whereas larger aggregates are nicely reproduced both in terms of structure and energy.
Non-stoichiometric clusters were also investigated. We selected (GaN) 14 and (GaN) 15 (shown in figure 3(b)) in parallel relaxation runs removing a single ion (Ga 3+ or N 3− ) in one cluster whilst inserting it into the other. Generally, the non-stoichiometric clusters based on QM and MM optimization are in reasonable agreement, except for DFT-optimized [Ga 14 N 13 ] 3+ for which we found that two four-member rings were changed to two six-member rings and vice versa ( figure  5(a)). In terms of energy, the equation for cohesive energy calculation cannot be applied to non-stoichiometric clusters since it is based on one formula unit. To make the comparison with the QM reference, we hence calculated the energy costs related to transferring ions between the clusters as shown in figure 5. We found that the MM-based reaction energies are in reasonable agreement with the DFT-based reference for both the transfer of a N 3− ion from (GaN) 14 to [Ga 15 N 14 ] 3+ (figure . This is in agreement with our findings related to single point defects in crystal that the formation energy of N 3− vacancy is lower than that of Ga 3+ vacancy by −0.5 eV for both the wurtzite and the zincblende structures, respectively. To illustrate the perspective of the molecular mechanics model, we decided to investigate the self-organization of larger (GaN) n (n = 1572) nanoparticles from molecular dynamics simulations of droplet annealing. The gallium and nitride ions were first randomly arranged (whilst respecting the van-der-Waals distances to avoid overlaps) into a spherical droplet. This randomized structure was modelled in vacuum and propagated at 5000 K for 5 ps to obtain a droplet model. The starting temperature of 5000 K is far above the reported melting points of GaN bulk crystals [10] (e.g. about 2493 K at pressure higher than 6 GPa [45]), and was chosen to ensure the annealing is subjected to liquid droplets. Using the molten GaN droplet as a starting point, annealing was then explored as a function of the cooling time. As seen in figure 6(a1), at shortest cooling time (0.005 ns) the energy of the resulting configuration is high and the final structure is less ordered, reflecting the solidification of amorphous nanoparticles. However, when we increased the cooling time, the configuration energy significantly decreases and finally converges. These longer quenching runs allow extended droplet reorganization and hence nucleation and growth of a crystalline interior embedded by compact nanoparticle surfaces. Fitting a mono-exponential function to the configuration energy as a function of cooling time results in a characteristic relaxation time of τ = 0.22 ns. From this fitting, the energy change from amorphous (instant quenching; zero cooling time) to ordered structures (infinite cooling time) is obtained as ( E ∞ − E 0 ) = −0.39 eV per formula unit ( figure 6).
Whilst the atoms at the surface form the six-membered ring layer with some bubble-like structures, a close inspection into the crystalline interior of the nanoparticle reveals that the structure is neither wurtzite nor zinc-blende. In fact, the structure in the interior is best described by four-membered rings and dominant six-membered rings, conforming to the body centred tetragonal (BCT) zeolite structure (figure 6(a2)). Xiao et al [46] also found this BCT-like structure in their study on surface transformation of GaN nanorods using molecular dynamics simulation with the force field developed by Zapol et al [22]. In addition, the BCT-like structures on Figure 6. (a1) Configuration energy calculated as a function of cooling time from a series of droplet annealing runs. Ga and N atoms are shown in brown and blue color, respectively. After quenching, further 1 ns relaxation at 300 K was performed for structural and energy sampling. The configuration energy was sampled from the last 0.5 ns thereof and a snapshot in the inset is rendered from the last frame of the MD simulation. A mono-exponential function, other wurtzite-type crystals have been predicted by ab-initio and classical molecular mechanics simulations, for example ZnO [47], ZnS [48], CdSe [49] and BN [50], and experimental evidences have also been verified for ZnO [48] and ZnS [51]. This atomic arrangement is expected to facilitate the transition of atomic arrangement between the interior and the surface during the quench runs. To underpin this phenomenon, we equilibrated the bulk GaN in the wurtzite-type phase at 300 K and then cut it into a spherical (GaN) 1572 cluster (figure 6(b1). After relaxing this model at 300 K for 1 ns, we found the interior is still arranged according to wurtzite structure, but the surface gets significantly disordered (figures 6(b2) and (b3)). Although the interior of the structure is arranged in the wurtzite structure, which is the most stable phase of GaN, this is compromised by the poorly ordered atomic arrangement at the surface. The energy sampled from the last 0.5 ns of the relaxation run is found as −84.67 eV/f.u., which is 0.14 eV/f.u. larger than the cohesive energy identified for the droplet annealing runs.
To further discriminate the role of bulk energy and surface stress on the nanoparticles, we modelled a periodic BCT type crystal of GaN by relaxing the unit cell as a function of pressure (figure 2). From this, we find the wurtzite structure clearly favored over BCT in the bulk crystal. The observed BCT type ordering of the nanoparticles hence suggests that for small nanoparticles the BCT structure offers a better compromise of bulk energy and a surface stress than the wurtzite analog. Accordingly, GaN is expected to follow a two-step nucleation mechanism with BCT as an intermediate structure.
The structural sequence is suggested as triggered by crystallite size R. Adopting classical nucleation theory here, the relation of the bulk energy term and surface energy changes as a function of crystallite size R-as the former scales by ~R 3 whilst the latter is proportional to ~R 2 .

Conclusion
The present study represents more than yet another GaN force field. Instead, we elaborated a fully transferrable molecular mechanics model that bridges from small (GaN) n aggregates to nano-crystals and all relevant solid phases. This was achieved from imposing formal +3 and −3 charges on the gallium and nitride ions, respectively, and by fitting Lennard-Jones and Born-Mayer-type potentials to mimic the covalent aspect of Ga-N bonds. The new force field received rigorous benchmarking with respect to solid phase polymorph structures, elastic properties, lattice energy, and intrinsic point defects. Moreover, comparison to reference DFT calculations shows the assessment of nanoparticles and (GaN) n aggregates with excellent accuracy for n > 6. This paves the way to molecular simulations of GaN nucleation and growth-either from solution (which involves amide and imide precursors [24]), or via deposition from the vapor. We note that systematic analyses of GaN nanoparticle morphogenesis and phase stability exceed the scope of the present work and remain to be addressed in future studies. Nevertheless, we already demonstrate the preference of BCT structure domains in small crystallites, whilst the bulk material is dominated by the wurtzite structure.