Aiming to shed a light in the analysis of the controversial experimental and theoretical results concerning the structural phase transitions of gallium phosphide (GaP), we performed a combination of molecular dynamics simulations (MD) and ab-initio calculations based on density functional theory (DFT) with vibrational corrections within the harmonic approximation, for pressures up to 50 GPa. We show that the transition sequence depends on temperature for both cases. For low temperatures, our MD simulations predict a transition sequence zinc blende ⇒ β-tin ⇒ Immm ⇒ NaCl, while DFT
predicts the sequence zinc blende ⇒ Sc16 ⇒ β-tin ⇒ Immm. Our DFT calculations with vibrational corrections also indicate that, as we increase the temperature, both formalisms feature the same results: the Sc16 phase gradually loses stability, justifying its absence in experimental results despite early theoretical predictions, and β-tin and NaCl become more stable. The kinetic effects are also studied by dynamically changing the pressure in MD simulations up to 200 GPa. We show that, as we increase the pressure, a transition to the β-tin structure should not occur at room temperature due a large energy barrier.  For higher temperatures, we show that the first transition is zinc blende ⇒ NaCl. The differences in the observed transition sequence at low and high temperatures indicate that the vibrational corrections included in the DFT calculations are fundamental to the description of the stability of higher pressure structures, as also observed on other semiconductors. Moreover, the good agreement between DFT and MD shows the suitability of the classical effective many-body interaction potential, used to describe the different high-pressure phases of GaP.