The relative thermodynamic stability of the four BN polytypes at different fixed pressures p and temperature T is evaluated by minimization of the Gibbs free energy G (T, p); (1) where E0 (V) is the total energy at T = 0 directly obtained from the DFT calculations, as described in section II A. Fvib (T, V) is the Helmholtz free energy contributed by the lattice vibrations (phonons), given by; ω(q,ν,V) is the phonon frequency at the wave vector q and the band index ν. ℏ and kB are the reduced Planck constant and the Boltzmann constant, respectively. It is worth noting that, in addition to q and ν, the phonon frequencies ω in Eq. 2 are volume dependent, as the phonon calculations are performed at the quasi-harmonic level using the PHONOPY package for phonon calculations , in which the Parlinski-Li-Kawazoe method with a finite displacement of 0.01 Å is used to calculate the force constants within sufficiently large supercells, i.e., 4 × 4 × 1 primitive hexagonal unit cells for h-BN (64 atoms), 4 × 4 × 1 conventional hexagonal unit cells for r-BN (96 atoms), 2 × 2 × 2 conventional cubic unit cells for c-BN (64 atoms), and 3 × 3 × 3 primitive hexagonal unit cells for w-BN (108 atoms). The supercells are then sampled with the 21 × 21 × 21 Monkhorst-Pack q-point grids in order to secure the convergence of the phonon frequencies and thus Fvib (T, V). To determine the term pV and include it in G (T, p) at fixed temperatures, the sums of the first two terms on the right-hand side of Eq. 1 at different fixed volumes are fitted to the third-order Birch-Murnaghan equation of state (EOS) . The pressure p is thus calculated by;