Questions on dihedral angle calculations

[Annika Viswesh]

Dear all, 

I am trying to use MDAnalysis to calculate the dihedral angles (chi1, chi2, etc.) of my residue at various timesteps and have some questions. 

To better understand the manner in which MDAnalysis calculates Dihedral angles, I looked at the source code, in particular the function _calc_dihedral_angle formula in calc_distances.h. 

I also looked into the formula for calculating dihedral angle in Wikipedia and other sources, and I noticed the formula for calculating the dihedral angle is slightly different than how MDAnalysis does it. Here are the differences and my questions:

1. In most formulas, to calculate the normal vectors, n1 = va x vb and n2 = vb x vc, where va = coord2 - coord1 and vb = coord3 - coord2 and vc = coord4 - coord3. 

However, MDAnalysis uses n1 = -va x vb and n2 = -vb x vc, so can you let me know why MDAnalysis does it this way? 

2. Some versions of calculating the dihedral angle formula state that 

alpha (our dihedral angle) = atan2(norm(cross(n1, n2)), dot(n1, n2))

However, MDAnalysis calculates alpha in the following way (according to my interpretation of the source code)

alpha = atan2(dot(cross(n1, n2), (vb/|vb|)), dot(n1, n2))

After drawing a few diagrams to convince myself, the angle between cross(n1, n2) and (vb/|vb|) must be nonzero, so norm(cross(n1, n2)) should not equal dot(cross(n1, n2), (vb/|vb|)),

So my question is - why does MDAnalysis calculate the numerator as dot(cross(n1, n2), (vb/|vb|)) instead of norm(cross(n1, n2))?

I would greatly appreciate it if you can help answer my questions.

Thank you for your help.

Sincerely,

Annika

[Fiona Naughton]

Hi Annika,

> However, MDAnalysis uses n1 = -va x vb and n2 = -vb x vc, so can you let me know why MDAnalysis does it this way? 

I'm not actually sure why we do it this way, though someone else might know better than me. From a glance through history, seems it was changed to to this some 10 years ago (https://github.com/MDAnalysis/mdanalysis/commit/5bd25e7e1ae874e8080ba6f7700647e0c3253d4b) - perhaps Richard might remember the reasoning...?

 

However, given the final angle calculation involves cross(n1,n2) and dot(n1,n2), as far as I can tell the signs should cancel out and we get the same answer.

> So my question is - why does MDAnalysis calculate the numerator as dot(cross(n1, n2), (vb/|vb|)) instead of norm(cross(n1, n2))?

cross(n1,n2) and vb/|vb| should be parallel (or anti-parallel); the 3D geometry is also doing my head in 😅, but if it helps to think it this way: n1 x n2 will be perpendicular to both n1 and n2, and since n1 = va x vb and n2 = vb x vc, vb is also perpendicular to both n1 and n2 - in the general case where n1 and n2 aren't themselves parallel, this will only be possible if vb and n1 x n2 are parallel.

As for why we use the dot product there rather than the norm: norm(n1 x n2) will always be positive, so the resulting alpha is in the range 0-180 and won't distinguish between dihedrals representing a clockwise or anticlockwise rotation. With the dot product, its sign changes depending the direction of rotation, so we get alpha in the range -180-180 and can distinguish them.

Hope that helped!

 - Fiona

[Annika Viswesh]

Hi Fiona,

Thank you so much for the clarification - yes this helped!

Sincerely,

Annika

--
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/72L7S1dnoI0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/b7cf847f-d372-4efc-b7ae-de7be348fd3dn%40googlegroups.com.