Bug #1109

gmxcheck and GROMOS 96 bonds

Added by Reid Van Lehn about 8 years ago. Updated about 8 years ago.

analysis tools
Target version:
Affected version - extra info:
Affected version:


I was using gmxcheck to check a simple trajectory of surfactants, and get many errors of the type "Distance between atoms X and Y is 0.153, should be 0.023" where the actual distance is always approximately correct and the "correct" value is the square of the actual distance. This always occurs for G96 type bonds, where I am defining the bonds using the macros in the ffbonded.ff file for the gromos53a6.ff force field (e.g. gb_27 in that example). A typical line from my topology is then:

[ bonds ]
3 4 2 gb_27

I also noticed the same errors when using gmxcheck on the output of the energy minimization in Justin Lemkul's KALP-15 tutorial, which also uses GROMOS96 bonds. I did not get errors for other steps in that tutorial, presumably because the bonds are constrained using LINCs.

I am using Gromacs 4.5.5 and it appears that the relevant lines (212 - 228) from the source code in gmxcheck.c are:

switch (ftype) {
case F_BONDS:
case F_G96BONDS:
b0 = idef->iparams[type].harmonic.rA;
case F_MORSE:
b0 = idef->iparams[type].morse.b0;
b0 = idef->iparams[type].cubic.b0;
case F_CONSTR:
b0 = idef->iparams[type].constr.dA;

It's not clear to me why b0 would be squared here since the switch statement doesn't look to discriminate between function 1 and function 2 bonds (if I interpret this correctly).

Justin Lemkul suggested that this is a small output bug and should be reported here.

Associated revisions

Revision 852a4ffa (diff)
Added by Justin Lemkul about 8 years ago

Corrected gmxcheck output for Gromos96 bonds.

There was a small output bug since Gromos96 bonds are modified to have
the reference value stored as (b0)^2 rather than just b0.

Fixes #1109.

Change-Id: Ie4516515f7cacb35d262c2059b5e461962340fec


#1 Updated by Justin Lemkul about 8 years ago

The reference bond length is squared in convparm.c, in function assign_param(), due to the functional form of Gromos96 bonds, which uses squared terms (equation 4.36 in the manual). The code in gmxcheck.c simply reads these values, which have been stored in the .tpr file. Can you modify gmxcheck.c to the following, recompile, and see if it works:

    case F_BONDS:
      b0 = idef->iparams[type].harmonic.rA;
    case F_G96BONDS:
      b0 = sqrt(idef->iparams[type].harmonic.rA);

I think the solution is this simple, but I just want to be sure. If it works, I'll implement the fix in the development code.

#2 Updated by Reid Van Lehn about 8 years ago

I implemented the change exactly as Justin suggested above and it works correctly. If only they were all this easy, right?

#3 Updated by Justin Lemkul about 8 years ago

  • Category set to analysis tools
  • Status changed from New to In Progress
  • Assignee set to Justin Lemkul
  • Target version set to 4.6

Great, thanks. I have pushed the change into the development repo; hopefully things will be corrected in the next version.

#4 Updated by Erik Lindahl about 8 years ago

  • Target version changed from 4.6 to 4.6.1

I expect this patch to make 4.6, but just to have our priority list for 4.6 clear I'm bumping the target version since it's not a critical bug.

#5 Updated by Justin Lemkul about 8 years ago

  • Status changed from In Progress to Closed

Also available in: Atom PDF