GROMACS version: 2021.4
GROMACS modification: No
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=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!