Hi @andre. Please also note that the normal modes used to compute the Duschinsky with the duschinsky function need to be mass-weighted. I am not sure if Firefly prints out mass-weighted data but you can use:
Li = Li[6:].T * m[:, np.newaxis]**0.5
Lf = Lf[6:].T * m[:, np.newaxis]**0.5
and then call the duschinsky function. For more information please check the documentation here.
I also highly recommend using DFT or a post-HF method with a basis set with proper size to get the input vibrational data. Please let us know if you face any issues.