The energetic stability of FePO4 polymorphs (berlinite, heterosite, monoclinic, and high-pressure forms) is investigated combining first-principles and experimental methods. Calculations at the density functional theory (DFT and DFT+U) level performed using the local density approximation (LDA) and the generalized gradient approximation (GGA and GGA+U) yield different relative energetic stability for those forms of iron phosphate. To discern the appropriate computational methodology, we have measured the transport and magnetic properties of high pressure (HP)-FePO4; we found that it is an insulating compound with a room-temperature resistivity of 2.107 Ohm.cm, and Fe3+ ions in high spin configuration (t2g2eg3). Both LDA and GGA methods fail to reproduce the physical properties of HP-FePO4, which are well-predicted within the GGA+U framework. This method predicts that berlinite-FePO4 is the most stable form at ambient pressure, whereas HP-FePO4 is the stable form at pressures above 2 GPa. Heterosite and monoclinic FePO4 are predicted as metastable phases in the whole pressure range. Accordingly, we found that both monoclinic and berlinite transform to HP-FePO4 under pressure (pressure range 2-6 GPa, 900 oC). Differential thermal analysis (DTA) and temperature-XRD measurements reveal that HP-FePO4 reverts exothermically and irreversibly to the berlinite form at about 700 °C. Discrepancies of our combined experimental-computational results with a previous high-temperature oxide melt solution calorimetry investigation are discussed.