Is any bugs in gmx angle?

User discussions
1

GROMACS version:2020.2
GROMACS modification: No
Dear all gmx users and developers,
I want to use “gmx angle” the dihedral changes as time in my molecule, my order is :
gmx angle -f xxx.pdb -n index.ndx -type dihedral -ov
after I get the angaver.xvg file, I calculate the mean value of the dihedral using python3, but the mean value is greatly different from the given by the Gromacs. And I carefully compare the angles at different time, this is because the gmx calculated the dihedrals uncorrectly, but the angles in angaver.xvg is correct. I don’t find the reasons…
Could anyone help me?

2

What evidence do you have of a bug? How does your output differ from that of gmx angle?

3

Dr. Lemkul (@jalemkul ),

I confirm that there is a bug in gmx_angle.cpp file while calculating average as mentioned by @littlema . The bug mainly related to conversion to the correct range of angles (-180 to 180). I fixed the lines and recomplied, now I am getting the expected value. I checked gromacs-5.1.4 (which uses different algo to calculate average), which is okay and gives the correct value.

Correct: gmx_angle.cpp [line: 425-434] line fixed: 427, 430

425    else
426   { /* Incorrect  for Std. Dev. */
427    real delta, b_aver = correctRadianAngleRange(aver_angle[0]);
428       for (i = 0; (i < nframes); i++)
429        {
430           delta = correctRadianAngleRange(aver_angle[i]) - b_aver;
431           b_aver += delta;
432           aver += b_aver;
433        }
434   }

Buggy:

425    else
426   { /* Incorrect  for Std. Dev. */
427       real delta, b_aver = aver_angle[0];    // bug
428       for (i = 0; (i < nframes); i++)
429        {
430           delta = correctRadianAngleRange(aver_angle[i] - b_aver);  // bug 
431           b_aver += delta;
432           aver += b_aver;
433        }
434   }

We need to fix also the block right above else also (line: 408-424).
Sincerely,
Masrul

4

@Masrul is this still relevant in the version at hand (2020.2)?

5

Dr. @jalemkul Yes, bug is true for 2020.X

6

Please file an issue on https://gitlab.com/gromacs with your proposed fix.

7

yeah， I have checked the code of version 2020.2, it is same as Masrul said.

8

Thank you very much!! My dear friend! And I think Masrul can file the issue on the gitlab. He does the great job! @Masrul @jalemkul

9

Dr. Lemkul,

I thought about this, not sure if it is a bug. Maybe both one is correct.

For example:

dih1 = -160
dih2 = 90

Gromacs-5.1.4: Average = -35 ( Correct, considering algebraic mean )

Gromacs-2020.X: Average = 145 (Correct, considering periodicity of 360).

Obviously both are correct.

Line diagram: (The way, I see the problem)

|<------------------------------360 --------------------------------->|
dih1=-160----------------------------dih2=90-------------dih1[p]=200

avg(dih1,dih2)=-35
avg(dih2,dih1[p])=145

Now, I am wondering, probably, doing average in periodic quantities is not meaningful and rather ambiguous.

Let me know, what is your opinion?

Sincerely,
Masrul

10

Indeed positive and negative angles have important implications, and a “correct” numerical answer can have a totally wrong physical meaning. Imagine an even more extreme case, if you have two frames in which the dihedrals of interest are φ1 = 179° and φ2 = -179° according to gmx angle (or any other program). The value of <φ> is zero! The interpretation is that the dihedral is in the cis configuration while it is decidedly trans.

The correctRadianAngleRange function is to normalize the range of angles and put everything within {-180°,180°} so the user then has to interpret what this means in a consistent way.

1 Like
11

@jalemkul
I think I get it. In order to get the correct physical meaning, the φ2 should =360-179=181°, which is used in Gmx 2020.2 to get the mean value.And in the angaver.xvg file, it is still -179° to get the angle in the range {-180° +180°}.Am I right?

12

Yes, you need to make sure you’re applying a useful convention for the values. Each is equally valid but may be misleading.