I tried the two options (a) grad conp prop phon or b)opti free zsisa) you mentioned there, but I am using the pGFN-FF model.
Consequently, I could not run b) calculation, as the free method is not implemented for this model. For a), the calculation runs only if I don’t use the shrink keyword, which I understand it is needed so we can integrate across the Brillouin zone to obtain a converged free energy difference.

I am interested on calculating these free
energy differences for different phases of a framework. Is there any other way I can properly calculate the free energy using pGFN-FF in GULP?

Explicit k points aren’t implemented for pGFN-FF yet, but there are 2 ways you can compute phonon frequencies sampled at different k points to get the free energy:

Provided you’re using version 6.2 onwards then the code will automatically create a supercell to determine the real space force constants and then phase them to get the right k points & so shrink should work.

There is a direct equivalence between k points are supercells & so you can build larger supercells until you converge the free energy per formula unit & so this avoids the need to go away from the gamma point.

Thank you so much for our answer. I was testing both approaches.

a) without the shrink keyword, what caught up my attention was that, considering only one cell of my smaller system (a ~ 13 Å, b ~ 20 Å and c ~ 35 Å, lattice angs = 90º, previously optimized at constant volume and based on a crystallographic structure), only the first 3 frequencies were imaginary modes with small magnitude (-0.04, -0.03 and -0.03). If I try to do the same calculation, but replicating the optimized cell (2x2x2) to check the Free energy convergence, a huge amount of imaginary modes appear (from -891.00 to -5.08).
The same occurred for other systems that I previously optimized at constant pressure.

b) Using the shrink keyword and gulp 6.2 for the same structure above: imaginary frequencies also appear, less than in the case of 2x2x2 cells, but ranging from -414 to -10. I tried different shrinking factors, keeping them proportional to the cell parameters, but these imaginary modes were always present.

Therefore, I can’t know for sure if only one cell is enough in a) case, or how to solve the imaginary frequencies problem. I also tried to re-optimize the systems, exploring other optimizers (rfo) and using the anneal keyword, but could not lose those imaginary freqs.

You can have imaginary frequencies at k points even the frequencies at gamma are fine. This usually indicates that the system wants to distort in a way that is incommensurate with the cell (i.e. it needs a supercell to get rid of the imaginary modes). Usually you look for the point in the BZ where the frequencies are most imaginary and this will tell you the supercell dimensions.
If you can post your input then it’s easier to check whether this is the case or if something else is happening.