Hi, everyone! I am a newbies for GROMACS and begin use it in this week. For test, I compare energy calculated from Charmm and Gromacs with CHARMM36-jul2021 force field parameters. And I found there is difference in bond enrgy term between Charmm and Gromacs, desipte of other terms. (for a terminated LEU residue, the difference in bond energy term is about 0.19 KJ/mol. I think it is a bit large).
After tracing bond energy calculation by print KrA, KrB and rA in bonds.cpp file, I found Kr’s non-zero value in dicimal place is missing. For example, terminal NH3 and HC bond type, its Kr in ffbonded.itp and processed topology file(generate by gmx grompp -pp top-file) is 337230.40 (" NH3 HC 1 0.10400000 337230.40 " ), while becomes 337230 in bond energy calculation. The reason is that its value becomes 3.37230e+05 in tpr due to tpr precision (parameter is printed by gmx dump -s tpr-file and the output look like this: “functype[121]=BONDS, b0A= 1.04000e-01, cbA= 3.37230e+05, b0B= 1.04000e-01, cbB= 3.37230e+05” ).
Is right my test? I think this is a general problem, if I am right. Are there some ways to fix this issue, like overriding force field parameter in tpr by reading processed topology file? It seems double-precision version Gromacs can’t fix this issue. Did somebody have a detailed comparision in energy between Charmm and Gromacs?
Another question about energy is that in calculation with charmm force field, is angle energy term merged/summed inot UB energy term??
looking forward to discussing with all of you!
Bests!
hi, hess! I agree with you！
I just stoped here and didn’t go further to check related code which writing tpr file and modify it to try to fix this issue.
wish someone farmiliar with gromacs code and interests in this issue will do this. hahaha~
bty, I think this fix which includes double precision compliation flag for grompp will very helpful for verifying the correctness of force field parameters file tranformed from other formats, like CHARMM or AMBER.
This is not a “fix”. This is a configuration/compilation option of GROMACS.
When we compare GROMACS to other codes, we always use double precision.
Note that the difference you see here is likely a worst case example. Mostly the CHARMM force field will be used with H-bonds constrained, so this contribution is gone. Also the difference in the parameter is smaller than what you write, as the print in gmx dump does not print the full float precision (but the total difference in the potential is what it is).
As you suggested, I perfome the benchmark again using double precision. Now it is ok. GROMACS can reproduce Charmm results very well, no matter bond energy or total energy.
Actually I used double precision in my tests. It improved than single float precision, but not close enough to charmm’s vaule. In my latest test, I do it more carefully, including zero the last three decimals and not centering the system. And also print each bond’s parameters and contribution to total bond energy.