Have you tried validating your MD approach (system setup, unit conversion, output interpretation etc.) using a simple model, say, a single water molecule? This seems to be a technical issue of this kind. I wouldn't expect two different minima to give this large differences in frequencies, particularly given that high-frequency modes are mostly local and thus largely independent of global conformational states (do they converge to similar values if you do not truncate the list, by the way?).

Another possibility is erroneous unit conversion when calculating force constants, but this you can probably easily check against standard FF values.

Hi, thanks for the response. It won't be unit conversion, as I've managed to recreate identical results from published work. My concern with FF comparison is that there is no exact details on the QM paramterisation i.e., grid integration size, whether it is tight, ultratight or normal. All of these will make a difference in the returned vibrational frequency. 

I am thinking, as I have already had molecules parameterised using my approach and they reproduce experimental data, without including their parameters (it's a different project), it would be best to correlate my derived Force Constants with something from the energy landscape. I've found a great correlation between the Laplacian Bond Order (using QTAIM's multiwfn software) with the magnitude of that bond's derived force constant. 

What do you mean exactly by QM parametrization - the QM setup of freq calculation? Even with possible contributions from anharmonicity, I would not expect the discrepancy to be so large. For reasonably small molecules, most of the time I was able to  reproduce IR spectra at least qualitatively using very simple QM models.

I'd still do some form of sanity check, e.g. identify localized normal modes such as C=O or C-H stretching and see whether you're getting reasonable values from both approaches. If only the low-frequency modes diverge, then what you say might be reasonable.

I am referring to reproducing the IR spectra of a QM model from an MD model using normal analysis mode in Gromacs. 

Generating the IR spectra is not a problem. It generating the corresponding spectra from MD model in vacuum that are being the problem.  

Ok, I see. Still, I do not know whether it's only the low freq modes that cause problems. Off the top of my head, I can think of at least several simple diagnostics:

By the way, how do you treat dihedrals/impropers in your MD model? You do realize that you cannot infer dihedral force constants from NMA. And if missing, they could indeed make your molecule quite 'loose', although I have no idea how much.

Hi Miłosz, thanks for those suggestions. I'll act on those. With regards to your very last point, you are quite right, which is why I am only calculating bond angle and bond stretch force constants, not dihedral.