There is no way to calculate the free energy within quasi harmonic debye model by PHONOPY. One should use GIBBS code to do that.

REF to GIBBS code: http://cpc.cs.qub.ac.uk/summaries/adsy

Your attempt should first be to obtain acceptable positive phonon frequencies. Negative frequencies have no meaning. These are actually indicating imaginary frequencies as square root of a negative quantity. If you have a solid, stabilize your structure, include anharmonic terms if necessary but frequencies must be positive. Once you have phonon frequencies, it is a trivial matter to claculate the Free energy even in harmonic approximation - just requires a summation of a sinh function of phonon frequencies over the Brillouin zone. No need to go to any professional coding.

Yes, but Quasi harmonic Debye model (QDM) does not need the phonons to calculate the free energy right? And that's why I am using that model for my calculations as I have unstable structures. Making an unstable structure stable is computationally too expensive if someone use SCAILD or Ab-initio MD which deals with phonon and QDM is a kind of analytic solution which only needs volume-energy data and Poisson's ratio as input. In this respect one can use GIBBS code as there are plenty of papers published on this code.

I cant advise you working with unstable structures. To obatin minimum configuration is basic and so also getting positive phonon frequencies- otherwise you while working with unstable structures do not know for what are you calculating and the results keep on changing.

Yes you're absolutely right, usually people don't work with unstable structures. But we are interested to capture the phase transition temperature between the competing phases of some element.

Now, Ti, Hf, Zr are the metals which have a different phase at elevated temperature which you can't see from DFT at 0K as they have different phase at 0K. So, what people does- calculate their free energy (of both the phases) and plot free energy difference as a function of temperature and at some point you will see that those will cross each other leading to the phase transition temperature.

N.B: Unstable structures at 0K (e.g. Ti, Zr and Hf beta phase) show imaginary frequency and hence you can't apply phonon calculations (e.g. QHA) but there are several ways to handle these imaginary frequencies to calculate the free energy. Such methods involve - 1. Self-consistency Ab-Initio  lattice Dynamics (SCAILD), fast free energy calculations, AIMD where these methods work with phonon DOS and make the low temperature structure stable at some elevated temperature via entropy stabilization, anharmonic effect, large phonon displacements etc. followed by the calculation of the free energy but QDM is somewhat analytic, fast and produces reasonable results with experiments.

If you still have some query, please let me know Sir!

Capturing phase transition temperature based on free energy is fine and that is often done by finding out the temperature at which the two structures give same free energy. but these these structures must be always stabilized for minimum energy and thus phonon frequencies calculated. If you have different solids you will get positive phonons and you have to minimize the equilibrium config carefully. If you can proceed without phonon frequency approach, I have no idea. But I would wonder how would you call claim a structure without minimum energy or stable or positive phonons.

But how would you get minimum energy for two phases of an element at a particular temperature? I mean one of the phase gives minimum energy right? so what abut the second one? then how would you compare their free energy to get transition temperature?

E.g:, Ti has alpha and beta phase and beta phase is only stable after 1155K. Then how would you get free energy of beta phase below 1155K as it has imaginary phonon frequencies?

So, that's why it is stabilized at some elevated temperature by entropy stabilization or by some other means explained earlier. But you can only have minimum energy for only one phase and in that sense the other phase will always be unstable.

It is better to extrapolate the temperature of equivalent free energy of two different structures in their minimum energy configurations. Their slopes will be different.

Although I am using QDM model for the moment, but I don't like analytic equations as it involves less physical significance unlike phonons. For the accurate calculations, we will do AIMD or something.

Similarly extrapolating temperature I think also will be crude approximation (it's just a technical thing). But could you please give me some references where they used extrapolating temperature method for calculating free energy of unstable structures? And does this technique really reproduce experimental data? Would be good to have literatures on this as I have no idea about this.