The first principles total energies E0 at T = 0 K as a function of volume V of the four BN polytypes are derived from the density functional theory (DFT), in which the project augmented wave (PAW) method , as implemented in the Vienna ab initio simulation package (VASP), is used and the generalized gradient approximation (GGA), as proposed by Perdew, Burke, and Ernzerhof , is employed for modeling the exchange-correlation interactions. For the four BN polytypes, the energy cutoff for plane waves, included in the expansion of wave functions, is set to 600 eV, and the Monkhorst-Pack k-point mesh of 15 × 15 × 15 is chosen for the Brillouin zone integration. Since the standard DFT calculations cannot accurately describe the van der Waals interactions, the DFT-D3 correction method with Becke-Jonson damping is used to account for the interactions. For each BN polytype, the optimal equilibrium volume V0 at T = 0 K, minimizing E0, is determined by the minimum point of the total energy curve E0(V), calculated for a set of different fixed V. During the calculations of E0 at each fixed V, the unit-cell shape and all internal atomic coordinates within the cell are fully relaxed so that the total force, acting on each atom, is less than 10 -6 eV/Å. Note that, for this particular case, the calculated E0(V) are all assured to be converged within an accuracy of 1 meV/atom with respect to both the energy cutoff and the number of k-points.