I have derived force constants for MD simulations from vibrational frequency calculations of QM Gaussian calculations. The QM to MD force constant conversion is from my own software. Now that I've implemented those force constants in Gromacs I want to validate my force constants back against those frequency values from Gaussian.
So I am using normal mode analysis in Gromacs i.e.,
Energy minimisation to < 0.001 force threshold via cg or L-BFGS.
Pulled out double precision from the trajectory file.
Integration using nm (normal mode) analysis to generate a force field hessian. Then pull the eigenfreq from the hessian.
The first few frequency values from my QM vib freq are (I've truncated):
18.3, 210.0, 342.9, 415.3, 475.0, 527.2, 653.6, 705.7, 742.8
The first few frequency values from my MD hessian are:
1.1, 1.2, 1.6, 1.7, 2.1, 2.4, 7.5, 46.1, 46.6, 61.2.
As you can see, very different.
I feel as though comparing two sets of frequencies from two energy minima is woefully inaccurate, for a start, I can't see how it would be possible to ensure that both structures are in the same minimum.
I would like to hear anyone's thoughts on how to compare back from MD structural data to QM structural data.
Thanks
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.
- Badges
- Science topic