How do I calculate a dihedral angle given Cartesian coordinates?

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.

-
Regarding the level the answer should be at, are you familiar with: addition and subtraction of vectors? dot product and cross product? the normal vector to a plane? – Rahul Jun 23 '11 at 4:09
Dot product I recall, cross a little less so, though I went and read the article on it, so a bit better. – Nick T Jun 23 '11 at 4:16

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.

1. Given the coordinates of the four points, obtain the vectors $b_1$, $b_2$, and $b_3$ by vector subtraction.

2. 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$.

3. 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.)

4. 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$.)

-
A note about the direction of $\phi$, according to the diagram (and the definition on wikipedia), it does not use the right-hand-rule along vector $b_2$. It instead uses the left hand rule about $b_2$. – stvn66 Oct 27 at 15:43

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.

-

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|}$

-
The range of arccos is only 0 to $\pi$, how do I account for the other angles? – Nick T Jun 23 '11 at 4:34
That's the angle between the vectors $b_1$ and $b_3$, not the dihedral angle, which as its name suggests is between two planes. – Rahul Jun 23 '11 at 4:45
@Nick : There are always two angles between two planes (just like between two vectors in $R^2$), one of them will be less than $pi$, and by convention it's the one we choose. There isn't really a way to give an "orientation" to angles like in $R^2$ so we can't do better than that and have a canonical way of choosing which angle we should measure, so we always measure the smaller one. – Vhailor Jun 23 '11 at 4:52
@Vhailor: In this case, the direction of $b_2$ does provide a way to choose an orientation. – Rahul Jun 23 '11 at 5:26

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⟩$$

-

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];


}

-