Hi Cerda,
I see that you have used a FF calculator, which is fine, but the geometry should be optimized to reach the nearest minimum in the potential energy surface before performing the dynamical matrix calculation. Without geometry optimization, your script is generating a lot of -ve normal mode frequencies. Given the size of the configuration, this could take thousands of optimization steps before the forces go below 0.01 eV/A. A good practice before computing the free energy is to check the vibrational frequencies using the 'VibrationalMode' analysis object. For nonlinear molecules, all eigenfrequencies except the first 6 should be positive. If not, then either the FF is not good enough or the geometry is not a local minimum to perform harmonic analysis.
Please try geometry optimization until the vibrational frequencies are as mentioned above before performing IdealGasThermoChemistry - negative frequencies result in nan during vibrational entropy calculation under ideal gas approximation. When you use PhononDOS, then we only consider the positive frequencies for calculation, that is why it does not give you nan. But as I mentioned earlier, PhononDOS uses Quasi Harmonic Approximation meant for solids and does not include contributions for rotation and translation. Also, it is a constant volume approximation - meaning that it only gives you the Helmholtz free energy and not the Gibbs free energy.
Let me know if there are further questions.