early stages of GaN(0001) epitaxial growth are modeled using a lattice-based kinetic Monte Carlo simulation. The simulation explicitly includes both the constituent atoms and is performed on the wurtzite crystal structure of GaN. The deposition flux is chosen to mimic conditions in a molecular beam epitaxy chamber. The surface diffusion barriers used for the growth simulation are obtained from ab initio density functional theory based calculations. The evolution of submonolayer islands from small random clusters to ordered triangular islands is captured in these simulations. Further the submonolayer island density statistics is calculated to illustrate standard scaling behavior. The island density as a function of coverageis used to show the different regimes of submonolayer growth and the transition to multi-layer growth. The island density dependence on deposition flux, Ga:N flux ratio, and temperature are also shown to be reasonable and consistent with experiments. We further highlight the importance of N-rotation, which is a diffusion mechanism for nitrogen adatoms under moderate gallium excess. (C) 2017 Elsevier B.V. All rights reserved.