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?

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

?

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

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

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

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

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.

@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?

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