On the Role of Amides and Imides for Understanding GaN Syntheses from Ammonia Solution: Molecular Mechanics Models of Ammonia, Amide and Imide Interactions with Gallium Nitride.

A key requisite to characterizing GaN precipitation from ammonia solution from molecular simulations is the availability of reliable molecular mechanics models for the interactions of gallium ions with NH3 , NH2 - , and NH2- species, respectively. Here, we present a tailor-made force field which is fully compatible to an earlier developed GaN model, thus bridging the analyses of Ga3+ ions in ammonia solution with the aggregation of [Gax (NH)y (NH2 )z ]+3x-2y-z precursors and the modelling of GaN crystals. For this, quantum mechanical characterization of a series of Ga-coordination clusters is used for parameterization and benchmarking the generalized amber force field (GAFF2) and tailor-made refinements needed to achieve good agreement of both structural features and formation energy, respectively. The perspectives of our models for larger scale molecular dynamics simulations are demonstrated by the analyses of amide and imide defects arrangement during the growth of GaN crystal faces.

Introduction

For the study of ion solvation and ion‐ion interactions in solution, molecular simulations have proven as an indispensable tool of atomic‐scale in‐situ investigation complementing experimental characterization. While the major body of these studies is dedicated to ion solvation in water[ , ] there is an increasing number of molecular simulations aiming at the characterization of less‐common solvents such as ammonia. Similarly, also extreme conditions including autoclave setups are becoming assessable to molecular dynamics simulations.

For the specific task of understanding the syntheses of GaN from ammonia (e. g., via ammonothermal routes ), or more generally from gallium amide and imide precursors, molecular simulations encounter a two‐fold challenge. On the one side, best accuracy in treating the atomic interactions would be obtained from quantum chemical approaches. On the other hand, statistical sampling of extended solvation structures and reliable models of the bulk liquid/super‐critical phase embedding the solvated species calls for larger system sizes and longer time scales than currently accessible to first‐principles techniques.

So far, molecular mechanics models proved successful for the analyses of bulk ammonia and the solvation of mono‐ and divalent ions at quite reasonable accuracy – both in terms of energy and statistical sampling.[ , , ] For trivalent ions, namely for gallium, combined QM/MM studies showed that solvation in ammonia implies proton transfer reactions such that the prevailing solvated species is the [Ga(NH2)4] complex unless (extremely) acidic conditions are imposed.[ , , ]

While gallium amide complexes are hence readily available in ammonia solution, the formation of imides and nitrides requires aggregation processes which are hindered by kinetic barriers. From the ongoing efforts in rationalizing GaN syntheses a series of gallium imide/amide complexes was suggested for the stepwise condensation of [Ga(NH2)4] complexes upon aggregation. Moreover, even a bulk Ga2(NH)3 gallium imide compound was suggested as a metastable precursor material for gallium nitride.

In what follows, we revisit the molecular mechanics modelling of Ga⋅⋅N interactions for Ga solution complexes in ammonia, thus considering NH3, amide and imide ions, respectively. On this basis, improved GAFF parameters[ , ] are obtained and combined with our earlier force‐field for bulk gallium nitride to enable ns scale molecular dynamics simulation of (NH2)− and (NH)2− interactions in bulk GaN. Finally, we demonstrate the unprejudiced analyses of amide and imide defects in GaN crystal faces from simulated annealing runs.

Methods and Models

Ab‐initio Calculations

The quantum calculations were performed using the Gaussian 09 package using the LANL2DZ basis set with effective core potential for Ga, whilst all other atoms are described by the 6–31+G(d,p) and 6–311++G(d,p) basis sets. For the geometry optimization runs, we first employed the inexpensive 6–31+G(d,p) basis and the B3LYP functional.[ , ] Next, the resulting structures are subjected to single‐point calculations using 6–311++G(d,p) to provide more accurate assessment of energies, vibrational frequencies and partial charges. For such single‐point calculations of geometry‐optimized structures, electron correlation is treated at the Moller‐Plesset perturbation theory (MP2) level. The choice of B3LYP and MP2 was based on quantum studies of GaN nanostructures and GaN growth from the gas phase reported in Refs. [18-20]. The partial charges were assigned by the Antechamber module of the Amber tool to mimic the electrostatic potential by RESP charges.

Molecular Mechanics and Molecular Dynamics Simulations

All molecular mechanics calculations including the molecular dynamics (MD) simulations based on the force fields were performed using the Large‐scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) program. For the gallium nitride interactions we used a recently developed force‐field model that is based on full +3 and −3 partial charges, respectively. This enables consistent modelling of the Ga3+ interactions with imide (‐2 charged), amide (‐1) and ammonia species, respectively, as discussed in the results section. To assess the electrostatics, we use explicit Coulomb summation for isolated complexes, whilst the damped‐shifted potential with a damping parameter of 0.05 Å−1 was chosen for MD simulations of periodic systems. A time‐step of 0.5 fs was found adequate for the MD simulations to provide numerically stable treatment of the fast intramolecular vibrations of the (NH2)− and (NH)2− ions at a temperature of 300 K. To provide model robustness at higher temperature, during the artificial melting (up to 8000 K) and simulated annealing runs we employed the shake algorithm to restrict the bond lengths and angles of the (NH2)− and (NH)2− ions to their equilibrium values. On this basis, the MD time step could be increased to 1 fs. However, once our models are annealed, we switch back to the flexible force‐field for the final steps of geometry optimization.

Results

Intramolecular Mechanics of (NH2)− and (NH)2− Ions

To develop flexible models of the (NH2)− and (NH)2− species, we first performed geometry optimizations and identified the harmonic frequencies from quantum calculations. Next, simple harmonic potentials of the form:

were fitted to reproduce the equilibrium geometries and the vibrations as denoted in Tables S1 and S2 of the supporting information. From this, a quite satisfactory agreement of the harmonic frequencies is achieved. Moreover, the comparison of ammonia, amide and imide in terms of the N−H bond distances and the force constants of the corresponding bonding potentials nicely reflects a softening of the N−H bonds stemming from the space demand of the electron lone pairs of the nitrogen atom in NH3, (NH2)− and (NH)2−, respectively.

Intermolecular Interaction Models

The common GAFF2 force‐field is well suited for modelling ammonia and ammonia solutions of halides, mono‐ and divalent metal ions.[ , ] Moreover, we recently extended the GAFF2 force‐field for modeling the interactions of ammonium and amide ions in ammonia solution. On this basis, an intuitive estimate of (NH)2− interactions with NH3 and (NH2)− is given by adopting the Lennard‐Jones parameters from the amide model to describe the van‐der‐Waals interactions of imides. Along this line, we employ the Lorentz‐Berthelot combination rules for mixing the van‐der‐Waals parameters:

A full overview of the molecular mechanics models used for the N−N, N−H and H−H interactions is given in Table 1.

Table 1

Lennard‐Jones parameters employed for ammonia, amide and imide molecules. The partial charges q of the amide/imide species were derived from quantum calculations, respectively. For the (NH)2− species, our initial approach was to adopt the Lennard‐Jones parameters of (NH2)−. However, along with the creation of the tailor‐made force‐field, also the imide model was refined.

	qi /e	ϵii /eV	σii /Å	reference/note
NH3
N	−1.09539	4.118×10−4	4.0447	GAFF2 [11]
H	+0.36513	4.335×10−4	1.1065	GAFF2 [11]
(NH2)−
N	−1.59398	3.034×10−3	3.4102	modified GAFF2 [25]
H	+0.29699	4.335×10−5	5.8598	modified GAFF2 [25]
(NH)2−
N−N	−2.24261	1.847×10−6	8.0825	this study
H−H	+0.24261	0	0	this study

q/e

GAFF2

GAFF2

modified GAFF2

modified GAFF2

In turn, the modeling of trivalent metal ions in ammonia solutions is less straight‐forward because of the strong interactions of the metal ion and the lone pairs of the coordinating nitrogen atoms. For the Ga3+ species, an approximate approach to modeling Ga−N interactions is i) to employ molecular mechanics models originally designed for bulk gallium nitride. On the other hand, ii) the formulation of tailor‐made Ga−N interaction potentials with individual parameters for ammonia, amide, imide and nitride could provide more detailed account of the different coordination of Ga3+. See Table 2 for an overview of the different molecular mechanics models used for describing the Ga−N interactions.

Table 2

a Ga−N interaction model based on directly adopting the van‐der‐Waals term from the Ga3+−N3− potential developed for bulk GaN. b Specific Lennard‐Jones parameters optimized for describing the interactions between Ga3+ and NH3, NH2 −, and NH2−, respectively. The models were optimized to best reproduce the structure and formation energy (see Table 3) of the Ga‐based complexes shown in Figure 1.

MMGaN(bulk) Ga3+−NH3/NH2 −/NH2−/N3−
Ga−N	A/eV 608.54	ρ/Å 0.435	Born‐Meyer potential, Ref. [13]
Ga−H	0	0	set to zero
new:	ϵ/eV		σ/Å
Ga−NH3
Ga3+−Nammonia	1.913		2.046
Ga3+−Hammonia	0		0
Ga−NH2 −
Ga3+−Namide	0.936		1.929
Ga3+−Hamide	0		0
Ga−NH2−
Ga3+−Nimide	1.615×10−3		3.132
Ga3+−Himide	0		0

Table 3

a Formation energies ΔE of the complexes shown in Figure 1 as obtained from QM calculations, the naïve Ga−N interaction model MMGaN(bulk) and our new force‐field designed to best reproduce the QM reference. The assessment of the MM‐based formation energies is shown for both, the QM‐optimized structures and individually relaxed structures using the MM models (numbers in square brackets), respectively. b Formation energies of the benchmark complexes shown in Figure 2. In analogy to Table 3a, the MM‐based formation energies are shown for both, the QM‐optimized structures and individually relaxed structures using the MM model (numbers in square brackets), respectively.

complex	ΔEQM /eV	ΔEMM/eV		percent error
		MMGaN(bulk)	new	MMGaN(bulk)	new
a	−37.81	0.56 [−10.52]	−34.13 [−35.87]	101.47 % [72.17 %]	9.73 % [5.13 %]
b	−67.20	−39.17 [−46.83]	−68.08 [−71.14]	41.71 % [30.31 %]	1.30 % [5.86 %]
c	−137.62	−90.60 [−103.64]	−139.36 [−146.62]	34.17 % [24.70 %]	1.27 % [6.54 %]
d	−142.06	−103.96 [−112.78]	−145.06 [−150.91]	26.82 % [20.62 %]	2.11 % [6.22 %]
Benchmark complex	ΔEQM/eV	ΔEMM, new/eV		percent error
1	−210.87	−218.73 [−225.93]		3.73 % [7.14 %]
2	−280.44	−277.37 [−288.05]		1.09 % [2.71 %]
3	−293.47	−293.77 [−300.95]		0.10 % [2.55 %]
4	−289.98	−286.57 [−295.01]		1.17 % [1.73 %]
5	−283.57	−276.45 [−286.30]		2.51 % [0.96 %]
6	−302.79	−306.26 [−311.73]		1.15 % [2.95 %]
7	−302.75	−301.47 [−310.59]		0.42 % [2.59 %]
8	−1099.37	−1115.99 [−1127.82]		1.51 % [2.59 %]

Figure 1

Relaxed structures of a) [Ga(NH3)6]3+, b) [Ga(NH2)4]−, c) [(H2N)3Ga(μ‐NH)Ga(NH2)3]2−, and d) [(H2N)2Ga(μ‐NH)2Ga(NH2)2]2− complexes optimized with QM and MM methods. Atomic representation is chosen as colored ball‐and‐stick, grey shade and empty circles to indicate the QM reference, the approximate MMGaN(bulk) model and our new force‐field, respectively. Colors‐ white: H, blue: N(ammonia), orange: N(amide), red: N(imide), brown: Ga .

To elucidate the former approach (i), in what follows denoted by the suffix “GaN(bulk)” to indicate the approximate nature of the Ga−N interaction model (Table 2a), we prepared a small series of Ga‐ammonia/amide and imide complexes as illustrated in Figure 1. Our choice of complexes was motivated by focusing as specific as possible on different aspects of gallium ion interaction with (a) ammonia, (b) amide and (c+d) imide ions. All of these complexes were identified as precursors during ammonothermal syntheses by Niewa and coworkers.[ , ] While complex a) is only observed at ammono‐acid conditions, we find it particularly useful for directly assessing gallium‐ammonia interactions. In turn, the aggregates b)–d) are adopted from experimental characterization of neutral to basic ammonia solutions. In addition to this, the complexes a) and b) have also been observed as structural motifs of precipitates from ammonothermal syntheses.[ , , ] For this series of aggregates we find the MMGaN(bulk) model to perform somewhat reasonable for the imide species. However, we observed increasingly large deviation of both structure and formation energy of aggregates comprising amide and ammonia molecules, respectively (Table 3, Figure1).

Figure 2

Relaxed structures of benchmark complexes for comparison of the new MM model with QM references. In analogy to Figure 1, atomic representation is chosen as colored ball‐and‐stick and empty circles to indicate the QM reference and our new force‐field, respectively. Colors‐ white: H, blue: N(ammonia), orange: N(amide), red: N(imide), cyan: N .

As a consequence, we decided to create a tailor‐made molecular mechanics model with individual Ga−N parameters for nitrogen atoms of ammonia, amide and imide, respectively. For this, conventional Lennard‐Jones type potential energy functions were used to best reproduce the QM references of the complexes a–d as shown in Figure 1. In a step‐by‐step manner, we first created Ga−N interaction models for ammonia and amide, by means of the a) [Ga(NH3)6]3+ and b) [Ga(NH2)4]− complexes, respectively. The corresponding parameters were optimized with respect to the force fitting procedure in GULP using all Ga−N forces at equal weight. This strategy was found appropriate to nicely reproduce both the equilibrium structures (Figure1) and formation energies (Table 3) of the reference set, respectively. Next, the imide‐amide complexes c) and d) were used to identify the Ga‐imide interactions. At this stage, only the imide parameters were optimized, whilst keeping the Ga‐amide force‐field as derived from the fitting of complex b). We tested different weights of the two imide‐amide complexes for the parameter fitting based on the Ga‐imide forces, and eventually found a 1 : 25 ratio to best reproduce both structure and formation energy of c) and d), respectively.

However, to demonstrate the transferability of our models, a wider set of benchmark complexes 1–8 as illustrated in Figure 2 was studied. To cover the whole variety of Ga−N and N−N interactions, this benchmark also includes nitride‐containing complexes. On the basis of the 4+8=12 structures investigated for parameter fitting and benchmarking, respectively, we find quite satisfactory agreement of both equilibrium structures and formation energies (Table 3). Indeed, formation energies are within <10 % error margins and the root‐mean‐square deviation (rmsd) of the MM models as compared to the QM reference are ∼0.1 Å for small (up to 3 Ga ions) and <0.5 Å for the larger (4 Ga ions) complexes, respectively. Most of the rmsd stems from orientation mismatch of the NH3, (NH2)− and (NH)2− species, in decreasing order of magnitude, respectively. Indeed, for all complexes investigated, the nearest‐neighbor Ga−N distances are nicely reproduced at average error margins less than 0.1 Å (Figures 1 and 2).

Modeling of (NH2)− and (NH)2− Defects in GaN Crystals and (0001) Faces

Ammonothermal syntheses of gallium nitride imply a cascade of amide and imide based intermediates, and hydrogen defects are thus of crucial importance for understanding nucleation and growth. Moreover, the properties of GaN semiconductors clearly depend on the occurrence and spatial distribution of defects.[ , , ] This particularly applies to amide/imide incorporation, as overall charge neutrality calls for Ga3+ vacancies, . To explore the interplay of Ga‐vacancies and H+ interstitials from theory, a common approach is to start from periodic models of a few unit cells and intuitively prepare a series of defect structures. After structural relaxation, typically employing DFT calculations, energy rankings help to short‐list favorable defect constellations. On this basis, Lyons et al. used model systems of a few tens of atoms to elucidate isolated Ga‐vacancies and combined H‐ complexes, including their interplay with the band structure of the semiconductor material.[ , ]

In turn, our MM model allows the application of ns‐scale MD simulations for exploring systems comprising 1000 s of atoms at quite feasible computational costs. Focusing on the structure and energy of H‐defects in GaN, this offers to study defect formation from computational experiments, mimicking crystallization from the melt as an unprejudiced approach to identifying favorable arrangements. For this, the assessment of large model systems is of crucial importance to simultaneously explore combined H‐ complexes and stand‐alone Ga vacancies and H‐interstitials, respectively.

As a demonstrator system, we choose for a (0001) crystal face of wurtzite GaN. This was constructed as substrate‐melt sandwich using a 2D‐periodic (0001) wurtzite slab with frozen Ga and N atomic positions as the substrate layer. The later was cut as a rectangular model by means of 8×14×3 replications of the hexagonal unit cell along the , and directions, respectively. The resulting slab was relaxed at 0 K, leading to model dimensions of 46.11, 46.58, and 15.76 Å3. After fixing the atom positions of the substrate layer, an analogously prepared GaN adlayer of 8×14×4 replications of the hexagonal unit cell was added. To achieve amorphization, the adlayer was temporarily compressed by 10 % and then propagated without restraints at high temperature. The latter was chosen as much as 8000 K, in order to obtain a highly dynamic melt in which good mixing is ensured even within the ns time scales of our MD runs. Indeed, such extreme conditions are needed to avoid any ordering of the melt/substrate interface. In turn, a repulsive wall was implemented 7.5 nm above the melt to avoid persistent evaporation of ions at such high temperature.

Using a small series of parallel simulation runs, we explored the annealing of a) the pure GaN system as described above, b) replacing 24 N3− ions by 24 (NH2)− in the melt, c) replacing 48 N3− by 48 (NH)2− and d) replacing 36 N3− by 12 (NH2)− and 24 (NH)2−, respectively. Moreover, in each of the systems b), c) and d) we removed 16 Ga3+ ions to impose charge neutrality. All ion removal/replacements were implemented for randomly picked atoms from the liquid film of (a).

After 5 ns propagation at 8000 K, the pure GaN substrate/melt interface (a) was subjected to annealing runs starting from 10 consecutive snapshots of the melt taken in intervals of 1 ns. For each of these starting points, cooling by 4, 2 and 1 K/ps to zero temperature was investigated. This provides 3×10 annealing runs which final configurations were analyzed in terms of structural features and overall surface energy. Illustrations of the pure GaN system (a) as taken before and after the annealing process are shown in Figure 3. Driven by the fixed substrate, we find the melt to crystallize in the wurtzite structure as expected for bulk GaN. However, at the outer surface of the GaN growth front we observed local deformations that lead to a ‘flattening’ of the hexagonal layers near the (0001) surface (Figure 3, right). This strongly resembles surface effects observed for the (0001) faces of ZnO nanorods which were attributed to the reduction of polarity at both, the anionic and cationic surface layers. To confirm this finding from a computational viewpoint, we compared the potential energy of the final structures as obtained from the 3×10 series of annealing runs. Figure 4 shows the surface energy Es as a function of the cooling rate, based on the difference of the potential energy of the crystal growth front and a 3D‐periodic bulk GaN crystal with the same number of ions, respectively. From this, cooling at a rate of 1 K/ps is suggested as a robust route to reasonably relaxed structures. Indeed, all 10 annealing runs performed at this rate lead to practically the same surface energy of 0.596±0.001 eV/Å2 and equivalent final structures as the one illustrated in Figure 3, right.

Figure 3

Representative snapshots of the pure GaN substrate‐melt interface before (left) and after (right) the annealing run. Upon cooling to zero Kelvin, the GaN melt adopts the wurtzite structure from the substrate layer. However, next to the interface with vacuum the puckered hexagons normal to [0001] prefer a flattened arrangement which leads to a smoothened surface layer. As a consequence, the outermost 2–3 (0001) layers show reduced dipoles as compared to the polarity of the wurtzite bulk. Colors‐ grey/green: Ga (fixed substrate), brown/cyan: Ga (mobile).

Figure 4

Surface energy of the (0001) face per area as a function of the cooling rate applied to the pure GaN model (a). The solid curve Es(ave) refers to averages taken over sets of 10 annealing runs each, whilst the error margins indicate the standard deviation, respectively. In turn, the dotted curve shows the surface energy Es(min) of the final configuration of lowest energy. Note that both curves converge at a cooling rate of 1 K/ps, thus demonstrating that all 10 annealing runs lead to equivalent structures.

Consequently, cooling by 1 K/ps was also applied to the annealing runs of the defect systems (b‐d) in order to find favorable arrangements of amide and imide species within the GaN growth front. While the starting temperature is chosen generously large to provide good mixing of the N3−, (NH)2− and (NH2)− ions, upon cooling we observe the segregation of the less charged amide and imide species to the surface of the melt. After solidification, we therefore find all of the (NH2)− ions and most of the (NH)2− at the surface of the crystal slab models. A statistical analysis of the different defect arrangements is provided in Table 4, whereas representative snapshots are shown in Figure 5, respectively.

Table 4

Defect statistics as obtained from averages over the 10 annealing runs performed for the models featuring (b) amide, (c) imide and (d) both amide/imide defects in (0001) crystal faces, respectively. While most of the (NH2)− ions were segregated to the GaN surface during crystallization of the melt, the (NH)2− ions remain in the crystal face and bulk phase. The numbers in square bracket denote the standard deviation. See also Figure 5 for illustrations of the annealed crystal faces.

	(b) – (NH2)−	(c) – (NH)2−	(d) – (NH2)−/(NH)2−
Total count: (NH2)−/(NH)2−/N3−	24/0/1992	0/48/1968	12/24/1980
(NH2)−/(NH)2− – in bulk	0 %/–	–/37 %	0 %/43 %
(NH2)−/(NH)2− – crystal face	17 %/–	–/63 %	27 %/57 %
(NH2)−/(NH)2− – segregated on surface	83 %/–	–/0 %	73 %/0 %
Defect arrangements in the bulk phase (total counts):
δ_N^3-→(NH_2)​^- – isolated	–	–	–
δ_N^3-→(NH)^2- – isolated	–	17 [±4]	10 [±2]
δ_N^3-→(NH_2)​^- ⋯□_Ga	–	–	–
δ_N^3-→(NH)^2- ⋯□_Ga	–	0 [±1]	1 [±1]
□_Ga – isolated	2 [±2]	1 [±1]	1 [±1]

Figure 5

Representative defect arrangements of the (0001) GaN crystal faces after re‐crystallization of (b) amide, (c) imide and (d) amide/imide containing melts. The least charged (NH2)− species is not incorporated in the wurtzite bulk, and may only be found in the outermost crystal layer or (preferred) in agglomerates located above the crystalline part of the simulation models. In turn, the imide species tends to incorporate into the outermost crystal layer, but also into the adjacent layer below. Only a small fraction of the imide ions show incorporation into the GaN bulk. See Table 4 for the statistics of all defect arrangements. Colors‐ white: H, orange: N(amide), red: N(imide), grey/green: Ga (fixed substrate), brown/cyan: Ga (mobile).

– isolated

–

–

–

– isolated

–

–

–

–

–

– isolated

While the solidification of the pure GaN model (a) leads to an ideally planar (0001) surface that reproduces the preparation of the 8×14×4 replications of the hexagonal unit cell (Figure 3), the models (b)–(d) deviate from this stoichiometry and hence exhibit rough/stepped crystal faces (Figure 5). As a consequence, these models allow insights into the interplay of amide/imide arrangements and the overall structure of the crystal face. However, in lack of a rigorous analysis of crystal growth as a function of ions deposited, our small selection of defect systems mainly offers qualitative insights and the quantitative information concluded from the ratio of different defect arrangements should be interpreted in terms of trends.

The observed preferences for the defects suggest that amide species are mainly found at the boundaries of small clusters (including those shown in Figures 1 and 2) and at the extremities of amorphous deposits on GaN surfaces. Likewise, we find the (NH2)− species at the corners of surface steps, whilst the (NH)2− ions tend to stabilize both surface steps and planar sketches of the surface. Moreover, we identified a small fraction of imide ions in the second layer beneath the surface. In turn, defects also exist in the GaN bulk, however mainly as isolated point defects, leaving charge neutralization to dispersed Ga vacancies.

This is somewhat surprising, as the +1 and +2 charges of the and defects intuitively lead to the expectation that amide and imide incorporation could be linked with nearby Ga3+ vacancies, . We thus examined the annealed structures for missing Ga species at positions suggested by superimposing the ideal GaN wurtzite structure. To characterize complexes of and defects, a distance delimiter of 2 Å was found appropriate to identify (at least one of) the hydrogen atoms of imide or amide species next to Ga vacancies, respectively. The few amide ions that were incorporated in the ordered arrangement of a step on the crystal surface are indeed all associated to a neighboring vacancy. However, only 10 % of the imide defects on the surface were found next to missing Ga atom, whereas many of our annealing runs showed not even a single complex in the bulk (Table 4).

Conclusions

The in‐depth understanding of GaN formation from ammonia solution poses an ongoing challenge to both experiment and theory. In the present work, we elaborated molecular mechanics models that will help to escape the limited time and length scales of quantum calculations whilst maintaining much of their accuracy. Indeed, we argue that our tailor‐made interaction model should provide reasonable assessment of [Gai Nj (NH)k (NH2)l (NH3)m]+3i−3j−2k−l aggregates, hence spanning the vast range of amide/imide containing intermediates from solvated Ga3+ ions in liquid ammonia to bulk GaN.

Benefitting from the ns time scales accessible to our models, we demonstrated simulated annealing studies of crystal formation. Triggered by a template slab, this allowed to identify aspects of (0001) GaN crystal faces, including preferences of amide and imide association/incorporation to rough/stepped surfaces, respectively. Incorporation of imide species ( defects) shows only minor effects of the GaN lattice, both in the bulk and at the (0001) surface. In turn, the association of amide species implies drastic reorganization of GaN crystal faces, and embedding in the bulk appears to be strongly disfavored. This suggests that the balance of nitride/imide and amide species – a parameter assessable to choosing pH conditions – strongly affects the expected growth rate and quality of GaN crystals. Rigorous confirmation of this calls for explicit consideration of proton transfer reactions as a function of pH. Indeed, as a further perspective of the present Ga‐nitride/imide/amide and ammonia models, we suggest the combination with QM/MM assessment of proton transfer events. This shall enable the modelling of GaN nucleation and crystal growth much in analogy to previous studies of ZnO nucleation and growth.[ , ]

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.