20 How do I calculate the Helmholtz (or Gibbs) free energy as a function of temperature keeping into account the thermal expansion?

Use what='mur_lc_t'. The Helmholtz (or Gibbs at finite pressure) free energy is shown as a function of temperature in the file output_anharm.ps. Note that the absolute value of this energy depends on the pseudopotentials as the total energy. You can however compare the free energies for different crystal structures and predict if a phase transition occurs and at which temperature (also as a function of pressure).