Conventional methods to calculate the thermodynamics of crystals evaluate the harmonic phonon spectra and therefore do not work in frequent and important situations where the crystal structure is unstable in the harmonic approximation, such as the body-centered cubic (bcc) crystal structure when it appears as a high-temperature phase of many metals. A method for calculating temperature dependent phonon spectra self consistently from first principles has been developed to address this issue. The method combines concepts from Born's inter-atomic self-consistent phonon approach with first principles calculations of accurate inter-atomic forces in a super-cell. The method has been tested on the high temperature bcc phase of Ti, Zr and Hf, as representative examples, and is found to reproduce the observed high temperature phonon frequencies with good accuracy.Many elements, alloys, and compounds appear in crystal structures which should not be energetically stable. The inter-atomic interactions places these systems at energy saddle points on the potential surface for atomic positions corresponding to the lattice sites of these structures rather than minima for statically stable structures. The body centered cubic (bcc) structure prevails as the simplest and best known example. Although a stable structure at low temperatures for several elements in the Periodic Table, bcc becomes unstable in the harmonic approximation [1,2,3] for the group IVB elements, for the rare-earth elements, for the actinides, and for several alkaline-earth elements. Nevertheless, at elevated temperatures the bcc structure emerges as the stable crystal structure for all these elements. Zener considered this enigma long ago and proposed a possible explanation: the large vibrational entropy of the bcc crystal structure makes it thermodynamically favourable at finite temperatures [5]. Also, Grimvall et al [4] pointed out the importance of electronic entropy in the stabilization of the bcc crystal structure of the group IVB elements Ti and Zr.So far no satisfactory, quantitative explanation has been presented for this situation. Density functional theory (DFT) [6] forms the basis of contemporary microscopic solid state theory and allows, in principle, to calculate different properties of crystals completely ab initio, without any fitting parameters. In particular, phonon spectra in the harmonic approximation can be efficiently evaluated in this way [7]. However, for the bcc phases mentioned above the phonon spectra in the harmonic approximation reveal imaginary phonon frequencies of e.g. Zr [8,9] for some wave-vectors, which shows that the bcc phase is from a lattice dynamics point of view unstable (hence these elements are energetically unstable, and are referred to as dynamically unstable in the bcc phase). A straightforward calculation using DFT molecular dynamics (MD) [10] should in principle be able to reproduce the stability of the bcc phase for the above discussed elements, since MD implicitly include temperature effects. However, MD suffers...