  An accurate and easily extendable method to deal with lattice dynamics of
solids is offered. It is based on first-principles molecular dynamics
simulations and provides a consistent way to extract the best possible harmonic
- or higher order - potential energy surface at finite temperatures. It is
designed to work even for strongly anharmonic systems where the traditional
quasiharmonic approximation fails. The accuracy and convergence of the method
are controlled in a straightforward way. Excellent agreement of the calculated
phonon dispersion relations at finite temperature with experimental results for
bcc Li and bcc Zr is demonstrated.
