Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
I have the $(x, y, z)$ coordinates for four points in space (atoms). How do I calculate the dihedral angle ($\phi$ in the below pictures from WP)?
I'm somewhat hazy on my math, and am trying to implement some function for this in Python, so a lower programmer's level would be appreciated.
It seems a single complete answer has not yet been given. Here is what I think is the most straightforward and numerically stable solution, expressed in a way that is easy to implement. It overlaps a lot with Vhailor's and Jyrki's answers, so I've marked it community wiki.
Given the coordinates of the four points, obtain the vectors $b_1$, $b_2$, and $b_3$ by vector subtraction.
Let me use the nonstandard notation $\langle v \rangle$ to denote $v/\lVert v \rVert$, the unit vector in the direction of the vector $v$. Compute $n_1 = \langle b_1 \times b_2 \rangle$ and $n_2 = \langle b_2 \times b_3 \rangle$, the normal vectors to the planes containing $b_1$ and $b_2$, and $b_2$ and $b_3$ respectively. The angle we seek is the same as the angle between $n_1$ and $n_2$.
The three vectors $n_1$, $\langle b_2 \rangle$, and $m_1 := n_1 \times \langle b_2 \rangle$ form an orthonormal frame. Compute the coordinates of $n_2$ in this frame: $x = n_1 \cdot n_2$ and $y = m_1 \cdot n_2$. (You don't need to compute $\langle b_2 \rangle \cdot n_2$ as it should always be zero.)
The dihedral angle, with the correct sign, is $\operatorname{atan2}(y,x)$.
(The reason I recommend the two-argument $\operatorname{atan2}$ function to the traditional $\cos^{-1}$ in this case is both because it naturally produces an angle over a range of $2\pi$, and because $\cos^{-1}$ is poorly conditioned when the angle is close to $0$ or $\pm\pi$.)
I think what you are looking for (if I understood the picture correctly) is the angle between the planes spanned by b1,b2 and b2,b3.
One way to do this is to calculate the vector angle between $b_1 \times b_2$ and $b_2 \times b_3$, where $\times$ is the cross product of vectors (you should be able to find the formula on wikipedia). You can calculate this angle with various methods, like the one mentionned by Ross. Depending on whether b1 and b3 lie on the plane of the second drawing, my answer and Ross's might give the same result.
Please correct me, if I interpret your question and/or picture incorrectly. Unfortunately that is a live possibility. Judging from the bottom figure I take it that you want the following angle. First you orthogonally project vectors $-\vec{b}_1$ and $\vec{b}_3$ to the plane that is orthogonal to the vector $\vec{b}_2$. The angle $\phi$ is then the angle between these two projections. Note the minus sign. This is necessary, because you seem to want to have $\phi=\pi$, when the two projections point in the same direction.
If this is the correct interpretation, then you can proceed as follows. The projection mapping is given by the formula $$ p(\vec{x})=\vec{x}-\frac{\vec{x}\cdot\vec{b}_2}{\vec{b}_2\cdot\vec{b}_2}\vec{b}_2. $$ Then you get the cosine of the angle $\phi$ from the usual formula $$ -p(\vec{b}_1)\cdot p(\vec{b}_3)=|p(\vec{b}_1)|\, |p(\vec{b}_3)|\cos\phi. $$ This gives you $\cos\phi$ and leaves you with a choice of two angles - one in the interval $[0,\pi]$ the other in $[-\pi,0]$. We need to know the sign of $\sin\phi$ to decide between the two. I read your bottom picture in the way the vector $\vec{b}_2$ is pointing away from the camera. This determines the orientation. We take advantage of this by computing the cross product $$ \vec{u}=p(\vec{b}_1)\times p(\vec{b}_3). $$ Then the inner product $\vec{u}\cdot \vec{b}_2$ has the same sign as $\sin\phi$ --- unless I made a systematic error :-). Let's check. If the angle in your picture is in the interval $[0,\pi]$, then I can twist the thumb of my right hand along $p(-\vec{b}_1)$, the index finger along $p(\vec{b}_3)$, so the middle finger will point in the direction opposite to $\vec{b}_2$. Therefore I need to use $p(\vec{b}_1)$ as opposed to its negative.
To clear up an ambiguity to the otherwise great response by Rahul, using the nonstandard notation for his unit vectors could give two possible answers.
For the orthogonal vector $n_2 =⟨b_2×b_3⟩$ this is not the the cross product normalized by it's components... $$n_2 = b_2×b_3$$ Such that... $$⟨n_2⟩=\frac{n_2}{∥n_2∥}$$ Rather, it is the cross product of the normalized vectors $b_2$ and $b_3$... $$n_2 = ⟨b_2⟩×⟨b_3⟩$$
You can use the dot product: $\vec{b_1}\cdot\vec{b3}=|b_1||b_3|\cos \phi$, so $\phi=\arccos \frac{\vec{b_1}\cdot\vec{b3}}{|b_1||b_3|}$
Solution I like the best is with atan2. Here is my c code for this:
//x1, y1, z1 coordinates of point 1 //x2, y2, z2 coordinates of point 2 and so on //b1, b2, b3 vectors described in text
//b1
b_a[0] = -(x1 - x2);
b_a[1] = -(y1 - y2);
b_a[2] = -(z1 - z2);
//b2
b_c[0] = x2 - x3;
b_c[1] = y2 - y3;
b_c[2] = z2 - z3;
//b3
c_d[0] = x4 - x3;
c_d[1] = y4 - y3;
c_d[2] = z4 - z3;
double n1[3];
double n2[3];
double m[3];
double x, y;
VectorNormalisation(b_c);
VectorNormalisation(b_a);
VectorNormalisation(c_d);
CrossProduct(b_a, b_c, n1);
CrossProduct(b_c, c_d, n2);
CrossProduct(n1, b_c, m);
x = DotProduct(n1, n2);
y = DotProduct(m, n2);
angle = 180.0 / PI * atan2(y, x);
//////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// //Functions used in the code:
double sqr(double x){ return x*x; }
void VectorNormalisation(double *v) { double lenght = sqrt(sqr(v[0]) + sqr(v[1] + sqr(v[2]))); v[0] /= lenght; v[1] /= lenght; v[2] /= lenght; }
double DotProduct(double *v, double *w) { return (v[0] * w[0] + v[1] * w[1] + v[2] * w[2]); }
void CrossProduct(double *v, double *w, double *cross) { //
cross[0] = w[1] * v[2] - w[2] * v[1];
cross[1] = w[2] * v[0] - w[0] * v[2];
cross[2] = w[0] * v[1] - w[1] * v[0];
}
