A two-phase molecular dynamics simulation of coexisting solid and liquid has been carried out to investigate the melting point of wurtzite-type GaN crystals. The melting point is determined by examining the movement of the interface between the solid and liquid during the simulation. The potential is a two-body interatomic one composed of the long-range Coulomb interaction, the Gilbert-type short-range repulsion, the covalent bonding and covalent repulsion of the modified Morse type, and the van der Waals interaction. The melting point and the interface morphology depend on the crystallization direction. The melting point Tm(K ) increases with pressure P(GPa ), but there appears a discontinuity in the vicinity of 8-9GPa. This is due to the solid-electrolyte-like behavior of Ga atoms with a partial charge in the high-pressure region. The discontinuity has not yet been confirmed by experiment. The least-squares fitted result is Tm=2538+177P -4.62P2 at pressures lower than 8GPa and Tm=2825+210P-5P2 at pressures higher than 9GPa. The Clausius-Clapeyron relation is confirmed using calculated thermodynamic data.