It is a widely used approach to calculate the equilibrium volume of a solid via minimization of its Gibbs free energy G=F+pV with respect to volume V. For an insulator F(V,T)=E(V)+Fph(V,T), where Fph(V,T) is the phonon contribution and E(V) is “static” elastic energy of a solid. E(V) is usually calculated via DFT by taking several values of the unit cell volume, which is obtained by changing the lattice constants. This approach is realized in codes such as Gibbs and QHA. It is straightforward for an isotropic solid (with cubic symmetry) since all three lattice constants change with temperature by the same relative amount i.e. Δa/a=Δb/a=Δc/c. However, for anisotropic solids the lattice constants may change by different amounts, thermal expansion becomes anisotropic and one has to use more general Gruneisen formalism as implemented e.g. here https://doi.org/10.1016/j.commatsci.2018.04.036 . I was not able to find codes implementing general Gruneisen formalism except this one https://github.com/gangliu-github/gruneisen-formula . However, this code just calculates anisotropic thermal expansion from given files of Phonopy calculations for strained supercells, and does not provide the strains that should be applied for a particular crystal structure. Therefore, my questions are:

1) Does somebody know codes implementing Gruneisen formalism for calculation of anisotropic coefficient of thermal expansion “from scratch” ?

2) Is the conventional approach described above (where the relative change of all lattice constants is the same) correct for anisotropic solids at least for calculating the average thermal expansion (and bulk modulus)?