We present a systematic ab initio study of the temperature- and pressure-dependent thermoelastic properties of hcp beryllium within the quasiharmonic approximation (QHA). The accuracies of the zero static internal stress approximation (ZSISA) and of the volume-constrained ZSISA that are widely applied in ab initio thermodynamic calculations are quantified. Particularly, the effect of ZSISA on the calculation of C11 and C12 is compared with a numerical approach which minimizes the free energy with respect to the atomic positions at each strain. In beryllium, minor deviations are found within ZSISA, which gives elastic constants (ECs) in good agreement with the full free energy minimization. A substantial difference is found between the QHA and the quasistatic approximation (QSA), with the former closer to experiments. Within the QSA, we compare the ECs obtained by interpolating from a set of geometries along the “stress-pressure” isotherm at 0 K (within the constant-volume ZSISA) with a more general interpolation on a two-dimensional grid of crystal parameters, which allows the calculation of the ECs along the 0 kbar isobar. This paper provides a practical approach for the investigation of the thermoelastic properties of hcp metals at extreme conditions.