1.21 How do I calculate the elastic constants as a function of temperature?

There are two ways. Both of them are a two step calculation. The fastest, but less accurate, method uses the ``quasi-static'' approximation. First use what='elastic_constants_geo'. This option computes the elastic constants at T = 0 K for all the geometries used by what='mur_lc_t' and saves them in the directory elastic_constants. In the second step, using the same input, run again THERMO_PW with what='mur_lc_t' after copying in your working directory the directory elastic_constants obtained in the previous step. Elastic constants are read only when lmurn=.FALSE.. The ``quasi-static'' approximation means that the code interpolates the T = 0 K elastic constants found in the first step at the geometry that minimizes the Helmholtz (or Gibbs at finite pressure) free energy at temperature T .
The second method uses the ``quasi-harmonic'' approximation and requires many phonon calculations at many geometries. First use what= 'elastic_constants_geo' and set use_free_energy= .TRUE.. This option computes the temperature dependent elastic constants taking as unperturbed geometries those used by what='mur_lc_t' and saves them in the directory anhar_files. In the second step, using the same input, run again THERMO_PW with what= 'mur_lc_t' after copying in your working directory the directory anhar_files obtained in the previous step. Elastic constants are read only when lmurn= .FALSE.. The ``quasi-harmonic'' approximation means that the code interpolates the temperature dependent elastic constants found in the first step as second derivatives of the Helmholtz free energy (corrected to give the stress-strain elastic constants), at the geometry that minimizes the Helmholtz free energy.
The plot of the elastic constants will be found in a file whose default name is output_anhar.el_cons.ps. Note that when the elastic constant are available in a file, with the option lmurn=.FALSE. the thermodynamic properties are calculated also for non-cubic solids. The main approximation of the present implementation is that the atomic coordinates are relaxed only at T = 0 K minimizing the energy and the free energy is not calculated as a function of the atomic positions. Please be careful to use the same tmin, tmax, and deltat in these calculations. Only a tmax larger in the what='elastic_constants_geo' calculation than in what='mur_lc_t' is allowed.